In C99, a header <complex.h> is provided that implements complex numbers.

As we've seen, a qubit is a two-state configuration with complex probability amplitudes distributed between the two states.
This gives us the following implementation of a single qubit: $$\lvert\psi\rangle = \alpha\lvert 0\rangle + \beta\lvert 1\rangle,\qquad \alpha,\beta\in\mathbb C$$

#include <complex.h>
 
typedef struct {
    double complex a; // alpha
    double complex b; // beta
} qubit_t;

$$\lVert\psi\lVert^2 = |\alpha|^2 + |\beta|^2 = \alpha\bar{\alpha} + \beta\bar{\beta}$$

double qubit_norm2(qubit_t psi)
{
    return creal(psi.a * conj(psi.a) + psi.b * conj(psi.b))
}

$$\lvert\psi\rangle\rightarrow\frac{\lvert\psi\rangle}{\sqrt{\langle\psi\lvert\psi\rangle}}$$

qubit_t qubit_normalize(qubit_t psi)
{
    double norm = sqrt(qubit_norm2(psi));
    return norm == 0.0 ? psi : { psi.a / norm, psi.b / norm };
}

$$\langle\phi\lvert\psi\rangle = \overline{\phi_\alpha}\psi_\alpha + \overline{\phi_\beta}\psi_\beta$$

double complex qubit_inner(qubit_t phi, qubit_t psi)
{
    return conj(phi.a) * psi.a + conj(phi.b) * psi.b;
}

$$P(0) = |\alpha|^2 = \alpha\bar{\alpha}$$

double qubit_prob0(qubit_t psi)
{
    return creal(psi.a * conj(psi.a));
}

$$P(1) = |\beta|^2 = \beta\bar{\beta}$$

double qubit_prob1(qubit_t psi)
{
    return creal(psi.b * conj(psi.b));
}

$$\lvert\psi\rangle\rightarrow\frac{P_i\lvert\psi\rangle}{\sqrt{\langle\psi\lvert P_i\lvert\psi\rangle}}$$ $$P_0 = \lvert 0\rangle\langle 0\lvert$$ $$P_1 = \lvert 1\rangle\langle 1\lvert$$

#define C_ONE  (1.0 + I * 0.0)
#define C_ZERO (0.0 + I * 0.0)
qubit_t qubit_measure(qubit_t psi)
{
    double r = (double)rand() / RAND_MAX;
    if (r < qubit_prob0(psi)) {
        return { C_ONE, C_ZERO };
    } else {
        return { C_ZERO, C_ONE };
    }
}

$$U = \begin{pmatrix}u_{00} & u_{01} \\ u_{10} & u_{11} \end{pmatrix}$$ $$\lvert\psi'\rangle = U\lvert\psi\rangle$$

typedef struct {
    double complex u00, u01;
    double complex u10, u11;
} gate2_t;
 
qubit_t gate2_apply(gate2_t U, qubit_t psi)
{
    return (qubit_t){
        U.u00 * psi.a + U.u01 * psi.b,
        U.u10 * psi.a + U.u11 * psi.b
    };
}

Pauli-X operator swaps amplitudes. It's a $180^\deg$ rotation around the X-axis of the Bloch sphere. $$X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$$

gate2_t gate_X(void)
{
    return (gate2_t){
        0, 1,
        1, 0
    };
}

Pauli-Z operator flips the relative phase. Probabilities remain unchanged, but interference patterns change. $$Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$$

gate2_t gate_Z(void)
{
    return (gate2_t){
        1, 0,
        0, -1
    };
}

Hadamard gate mixes the amplitudes evenly. If a qubit is in a pure state $\lvert\psi\rangle = \lvert 0\rangle$, after Hadamard states are evenly mixed $1/\sqrt{2}\cdot(\lvert 0\rangle + \lvert 1\rangle)$. If you apply Hadamard again, you return back to $\lvert\psi\rangle$, so Hadamard is its own inverse. This is the gate that enables interference-based algorithms. $$H = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$$

gate2_t gate_H(void)
{
    double s = 1.0 / sqrt(2.0);
    return (gate2_t){
        s,  s,
        s, -s
    };
}

This connects the Schrodinger equation with gate logic.

$$R_z(\theta) = e^{-i\theta\sigma_z/2} = \begin{pmatrix} e^{-i\theta / 2} & 0 \\ 0 & e^{i\theta/2} \end{pmatrix}$$

gate2_t gate_Rz(double theta)
{
    double complex e_minus = cexp(-I * theta / 2.0);
    double complex e_plus  = cexp( I * theta / 2.0);
 
    return (gate2_t){
        e_minus, 0,
        0,       e_plus
    };
}

$$\lvert\psi\rangle \sim e^{i\phi}\lvert\psi\rangle$$

qubit_t qubit_global_phase(qubit_t psi, double phi)
{
    double complex phase = cexp(I * phi);
    psi.a *= phase;
    psi.b *= phase;
    return psi;
}

$$\lvert\psi\rangle = \cos(\theta / 2)\cdot\lvert 0\rangle + e^{i\phi}\sin(\theta/2)\cdot\lvert 1 \rangle$$

void qubit_to_bloch(qubit_t psi, double *theta, double *phi)
{
    psi = qubit_normalize(psi);
 
    double a_mag = cabs(psi.a);
    *theta = 2.0 * acos(a_mag);
 
    double arg_a = carg(psi.a);
    double arg_b = carg(psi.b);
    *phi = arg_b - arg_a;
}