comnum-math 0.3.1

//! The Galois fields $\mathrm{GF}(2^{m})$: the arithmetic on which `comnum_fec` builds the
//! BCH and Reed–Solomon codes.
//!
//! # Model
//!
//! The field is generated by a primitive element $\alpha$, a root of a primitive
//! polynomial $p(x)$ of degree $m$ over $\mathrm{GF}(2)$: its $2^{m} - 1$ non-zero
//! elements are the powers $\alpha^{0}, \dots, \alpha^{2^{m}-2}$. An element is stored as
//! the $m$ bits of its polynomial representation (bit $i$ the coefficient of
//! $\alpha^{i}$), so that the addition is the exclusive or, and a table of the powers
//! and of the logarithms gives the product
//!
//! $$
//! \alpha^{i} \alpha^{j} = \alpha^{(i + j) \bmod (2^{m} - 1)}.
//! $$
//!
//! The minimal polynomial of $\alpha^{i}$ is the product over its conjugates,
//! $m_{i}(x) = \prod_{j} \left(x - \alpha^{i 2^{j}}\right)$, a binary polynomial.
//!
//! [`primitive_polynomial`] gives the conventional primitive polynomial of every degree
//! 2 to 16: the ones of the tables of Lin & Costello (Table 2.7) and of the coding
//! standards (DVB, CD, QR codes use $x^{8} + x^{4} + x^{3} + x^{2} + 1$ for $m = 8$).
//!
//! # Examples
//! ```
//! use comnum_math::gf2m::{Field, primitive_polynomial};
//! // GF(16) with p(x) = 1 + x + x⁴: α⁴ = α + 1 = 0b0011
//! let f = Field::new(4, primitive_polynomial(4).unwrap())?;
//! assert_eq!(f.pow(4), 0b0011);
//! assert_eq!(f.mul(f.pow(7), f.pow(12)), f.pow(4));          // α⁷ α¹² = α¹⁹ = α⁴
//! assert_eq!(f.minimal_polynomial(1), vec![1, 1, 0, 0, 1]);  // p(x) itself
//! # Ok::<(), comnum_core::Error>(())
//! ```
//!
//! # References
//!
//! 1. S. Lin and D. J. Costello, *Error Control Coding*, 2nd ed., Pearson, 2004, ch. 2
//!    (Table 2.7: primitive polynomials, Table 2.8: GF($2^{4}$)).

use comnum_core::{Result, ensure};

/// The conventional primitive polynomial of degree `m` (bit $i$ the coefficient of
/// $x^{i}$) for $2 \le m \le 16$, `None` otherwise.
///
/// Lin & Costello, Table 2.7, except $m = 8$ where the polynomial of the coding
/// standards (ETSI EN 300 744, ISO/IEC 18004, the compact disc) is given,
/// $x^{8} + x^{4} + x^{3} + x^{2} + 1$.
pub fn primitive_polynomial(m: u32) -> Option<u32> {
    const POLYS: [u32; 15] = [
        0x7,     // x² + x + 1
        0xB,     // x³ + x + 1
        0x13,    // x⁴ + x + 1
        0x25,    // x⁵ + x² + 1
        0x43,    // x⁶ + x + 1
        0x89,    // x⁷ + x³ + 1
        0x11D,   // x⁸ + x⁴ + x³ + x² + 1
        0x211,   // x⁹ + x⁴ + 1
        0x409,   // x¹⁰ + x³ + 1
        0x805,   // x¹¹ + x² + 1
        0x1053,  // x¹² + x⁶ + x⁴ + x + 1
        0x201B,  // x¹³ + x⁴ + x³ + x + 1
        0x4443,  // x¹⁴ + x¹⁰ + x⁶ + x + 1
        0x8003,  // x¹⁵ + x + 1
        0x1100B, // x¹⁶ + x¹² + x³ + x + 1
    ];
    (2..=16).contains(&m).then(|| POLYS[m as usize - 2])
}

/// The Galois field $\mathrm{GF}(2^{m})$ by exponential and logarithm tables.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Field {
    m: u32,
    exp: Vec<u32>,
    log: Vec<u32>,
}

impl Field {
    /// The field of degree `m` (2 to 20) generated by the primitive polynomial `poly`
    /// (bit $i$ the coefficient of $x^{i}$).
    ///
    /// # Errors
    /// [`comnum_core::Error::Config`] for a degree outside 2 to 20, a polynomial of another
    /// degree, or one that is not primitive (its powers of $x$ do not run through every
    /// non-zero element).
    pub fn new(m: u32, poly: u32) -> Result<Field> {
        ensure((2..=20).contains(&m), || format!("GF(2^m): the field degree must be 2 to 20, got {m}"))?;
        ensure(poly >> m == 1, || format!("GF(2^m): {poly:#x} is not a polynomial of degree {m}"))?;
        let q = 1usize << m;
        let mut exp = vec![0u32; 2 * q];
        let mut log = vec![0u32; q];
        let mut x = 1u32;
        for (i, e) in exp.iter_mut().enumerate().take(q - 1) {
            *e = x;
            ensure(i == 0 || x != 1, || format!("GF(2^m): {poly:#x} is not a primitive polynomial of degree {m}"))?;
            log[x as usize] = i as u32;
            x <<= 1;
            if x & (1 << m) != 0 {
                x ^= poly;
            }
        }
        ensure(x == 1, || format!("GF(2^m): {poly:#x} is not a primitive polynomial of degree {m}"))?;
        for i in q - 1..2 * q {
            exp[i] = exp[i - (q - 1)];
        }
        Ok(Field { m, exp, log })
    }

    /// The degree $m$.
    pub fn m(&self) -> u32 {
        self.m
    }

    /// The number of non-zero elements $2^{m} - 1$.
    pub fn order(&self) -> usize {
        (1 << self.m) - 1
    }

    /// $\alpha^{e}$ for any non-negative exponent.
    pub fn pow(&self, e: usize) -> u32 {
        self.exp[e % self.order()]
    }

    /// The exponent $e$ of a non-zero element $\alpha^{e}$ (0 for the element 0).
    pub fn log(&self, a: u32) -> usize {
        self.log[a as usize] as usize
    }

    /// The product $ab$.
    pub fn mul(&self, a: u32, b: u32) -> u32 {
        if a == 0 || b == 0 { 0 } else { self.exp[(self.log[a as usize] + self.log[b as usize]) as usize] }
    }

    /// The quotient $a / b$ of a non-zero `b` ($0 / b = 0$).
    ///
    /// # Panics
    /// In debug builds, for `b = 0` (a division by zero is a broken invariant of the
    /// caller: release builds return a meaningless element instead of checking it in the
    /// inner loops of the decoders).
    pub fn div(&self, a: u32, b: u32) -> u32 {
        debug_assert!(b != 0, "division by zero in GF(2^{})", self.m());
        if a == 0 {
            0
        } else {
            let n = self.order() as u32;
            self.exp[((self.log[a as usize] + n - self.log[b as usize]) % n) as usize]
        }
    }

    /// The inverse $a^{-1}$ of a non-zero element.
    ///
    /// # Panics
    /// In debug builds, for `a = 0` (as [`Field::div`]).
    pub fn inv(&self, a: u32) -> u32 {
        self.div(1, a)
    }

    /// The minimal polynomial of $\alpha^{i}$ (binary coefficients, lowest degree first).
    pub fn minimal_polynomial(&self, i: usize) -> Vec<u8> {
        // the conjugates α^{i 2^j}
        let n = self.order();
        let mut conj = vec![i % n];
        let mut e = (2 * i) % n;
        while e != i % n {
            conj.push(e);
            e = (2 * e) % n;
        }
        // Π (x - α^e) over GF(2^m)
        let mut p = vec![1u32];
        for &c in &conj {
            let r = self.pow(c);
            let mut next = vec![0u32; p.len() + 1];
            for (d, &a) in p.iter().enumerate() {
                next[d + 1] ^= a;
                next[d] ^= self.mul(a, r);
            }
            p = next;
        }
        p.iter().map(|&a| u8::from(a == 1)).collect()
    }
}

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

    #[test]
    #[cfg(debug_assertions)]
    #[should_panic(expected = "division by zero")]
    fn the_inverse_of_zero_is_a_broken_invariant_in_debug_builds() {
        let _ = Field::new(4, 0x13).unwrap().inv(0);
    }

    #[test]
    fn a_polynomial_of_the_wrong_degree_is_an_error_not_a_panic() {
        // x⁸ + x⁴ + x³ + x² + 1 for GF(16), and x³ + x + 1 for GF(16)
        assert!(Field::new(4, 0x11D).is_err());
        assert!(Field::new(4, 0xB).is_err());
        assert!(Field::new(4, 0).is_err());
        assert!(Field::new(8, u32::MAX).is_err());
        // x⁴ + x³ + x² + x + 1 has the degree but is not primitive (order 5)
        assert!(Field::new(4, 0x1F).is_err());
        assert!(Field::new(4, 0x13).is_ok());
    }

    #[test]
    fn every_default_polynomial_is_primitive_and_the_arithmetic_is_consistent() {
        for m in 2..=16 {
            let f = Field::new(m, primitive_polynomial(m).unwrap()).unwrap();
            assert_eq!(f.order(), (1 << m) - 1);
            let a = f.pow(5);
            let b = f.pow(f.order() - 2);
            assert_eq!(f.div(f.mul(a, b), b), a, "m = {m}");
            assert_eq!(f.mul(a, f.inv(a)), 1, "m = {m}");
            assert_eq!(f.log(f.pow(5)), 5 % f.order());
        }
        assert!(primitive_polynomial(1).is_none() && primitive_polynomial(17).is_none());
        // x⁴ + x³ + x² + x + 1 has order 5: not primitive
        assert!(Field::new(4, 0b1_1111).is_err());
    }
}