hqcr 0.1.0

// GF(2^8) arithmetic for Reed-Solomon.
// Primitive polynomial: p(x) = x^8 + x^4 + x^3 + x^2 + 1 (0x11D).
// Precomputed tables: GF_LOG[256] (discrete log base α) and GF_EXP[512]
// (antilog, doubled to avoid modular reduction in multiply).
// Exports: gf_add, gf_mul, gf_inv. Zero is handled as a special case.

// ── Table generation ──────────────────────────────────────────────────────────
//
// GF(2^8) has 256 elements. The non-zero elements form a cyclic group of
// order 255 generated by α (the root of p(x), i.e. the field element 0x02).
//
// GF_EXP[i] = α^i  for i = 0..254
// GF_LOG[α^i] = i  for i = 0..254  (GF_LOG[0] = 0, sentinel — never read)
//
// The table is doubled (size 512): GF_EXP[i] = GF_EXP[i mod 255].
// This lets gf_mul use GF_EXP[log_a + log_b] without a % 255 reduction,
// since log_a + log_b ≤ 254 + 254 = 508 < 512.

const fn build_tables() -> ([u8; 512], [u8; 256]) {
    let mut exp = [0u8; 512];
    let mut log = [0u8; 256];

    // x tracks the current power of α as a u16 so we can detect the x^8 term
    // without u8 overflow.
    let mut x: u16 = 1;
    let mut i = 0usize;
    while i < 255 {
        exp[i] = x as u8;
        exp[i + 255] = x as u8; // doubled copy
        log[x as usize] = i as u8;

        // Multiply by α = 0x02: left-shift, then reduce if x^8 term appears.
        // Reduction: x^8 ≡ x^4 + x^3 + x^2 + 1 (low 8 bits of 0x11D = 0x1D).
        x <<= 1;
        if x & 0x100 != 0 {
            x ^= 0x11D;
        }

        i += 1;
    }
    // exp[255] was already set above (i=0 wrote exp[0+255]=exp[255]=1).
    // Indices 510 and 511 remain 0 but are never accessed (max sum of two logs
    // is 508).

    (exp, log)
}

const TABLES: ([u8; 512], [u8; 256]) = build_tables();

/// α^i for i = 0..508. GF_EXP[i] == GF_EXP[i mod 255] for i < 510.
pub(crate) const GF_EXP: &[u8; 512] = &TABLES.0;

/// Discrete log base α. GF_LOG[α^i] == i for i in 0..255.
/// GF_LOG[0] is a sentinel (value 0); never pass 0 to functions that use it.
pub(crate) const GF_LOG: &[u8; 256] = &TABLES.1;

// ── Arithmetic ────────────────────────────────────────────────────────────────

/// Addition in GF(2^8) is bitwise XOR.
#[inline(always)]
pub fn gf_add(a: u8, b: u8) -> u8 {
    a ^ b
}

/// Multiplication in GF(2^8) via log/antilog tables.
/// gf_mul(a, b) = α^(log a + log b)  [with special case for zero].
#[inline]
pub fn gf_mul(a: u8, b: u8) -> u8 {
    if a == 0 || b == 0 {
        return 0;
    }
    GF_EXP[GF_LOG[a as usize] as usize + GF_LOG[b as usize] as usize]
}

/// Multiplicative inverse in GF(2^8).
/// a^{-1} = α^{255 - log a}  (since α^255 = 1).
/// Panics on zero in debug builds; undefined behaviour in release.
#[inline]
pub fn gf_inv(a: u8) -> u8 {
    debug_assert!(a != 0, "gf_inv(0) is undefined");
    GF_EXP[255 - GF_LOG[a as usize] as usize]
}

/// Division: a / b = a * b^{-1}.
/// Returns 0 if a == 0. Panics on b == 0 in debug builds.
#[inline]
pub fn gf_div(a: u8, b: u8) -> u8 {
    if a == 0 {
        return 0;
    }
    gf_mul(a, gf_inv(b))
}

/// Evaluate polynomial p (coefficients p[0] + p[1]·x + … ) at field element x
/// using Horner's method.
pub fn gf_poly_eval(p: &[u8], x: u8) -> u8 {
    let mut acc = 0u8;
    for &coeff in p.iter().rev() {
        acc = gf_add(gf_mul(acc, x), coeff);
    }
    acc
}

// ── Tests ─────────────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn exp_log_are_inverses() {
        // GF_EXP[GF_LOG[x]] == x for all non-zero x
        for x in 1u8..=255 {
            assert_eq!(GF_EXP[GF_LOG[x as usize] as usize], x, "x={x}");
        }
    }

    #[test]
    fn generator_order_is_255() {
        // α has multiplicative order 255: α^255 == 1, no smaller power does
        assert_eq!(GF_EXP[255], 1);
        for i in 1..255usize {
            assert_ne!(GF_EXP[i], 1, "α^{i} == 1 but order should be 255");
        }
    }

    #[test]
    fn mul_commutativity() {
        // Sample a grid; full 256×256 is only 65 536 iterations — cheap.
        for a in 0u8..=255 {
            for b in 0u8..=255 {
                assert_eq!(gf_mul(a, b), gf_mul(b, a), "a={a} b={b}");
            }
        }
    }

    #[test]
    fn mul_identity_and_zero() {
        for x in 0u8..=255 {
            assert_eq!(gf_mul(1, x), x);
            assert_eq!(gf_mul(0, x), 0);
            assert_eq!(gf_mul(x, 0), 0);
        }
    }

    #[test]
    fn mul_inv_roundtrip() {
        for x in 1u8..=255 {
            assert_eq!(gf_mul(x, gf_inv(x)), 1, "x={x}");
        }
    }

    #[test]
    fn known_values() {
        // α^1 = 0x02  (the generator is just 2)
        assert_eq!(GF_EXP[1], 0x02);
        // α^8 = x^8 mod p(x) = x^4+x^3+x^2+1 = 0x1D
        assert_eq!(GF_EXP[8], 0x1D);
        // α^255 = 1
        assert_eq!(GF_EXP[255], 0x01);
    }

    #[test]
    fn distributivity() {
        // a*(b+c) == a*b + a*c  for a few triples
        let triples = [(0x53u8, 0xCAu8, 0x7Fu8), (0x01, 0xFF, 0xAB), (0x80, 0x80, 0x01)];
        for (a, b, c) in triples {
            let lhs = gf_mul(a, gf_add(b, c));
            let rhs = gf_add(gf_mul(a, b), gf_mul(a, c));
            assert_eq!(lhs, rhs, "a={a:#04x} b={b:#04x} c={c:#04x}");
        }
    }
}