comnum-fec 0.3.0

//! Binary BCH codes, shortened: the outer code of DVB-S2 (`comnum_systems::dvbs2`).
//!
//! # Model
//!
//! In the Galois field $\mathrm{GF}(2^{m})$ generated by a primitive element $\alpha$
//! (root of a primitive polynomial $p(x)$), the narrow-sense BCH code of designed
//! distance $2t + 1$ has the generator polynomial
//!
//! $$
//! g(x) = \operatorname{lcm}\left(m_{1}(x), m_{3}(x), \dots, m_{2t-1}(x)\right),
//! $$
//!
//! with $m_{i}$ the minimal polynomial of $\alpha^{i}$ (Lin & Costello, ch. 6). The code
//! of length $2^{m} - 1$ is shortened to $n$ bits by fixing its first information bits
//! to zero. The encoder is systematic: the $k$ message bits $m_{k-1}, \dots, m_{0}$ (the
//! highest degree first, in transmission order) followed by the $n - k$ coefficients of
//! $x^{n-k} m(x) \bmod g(x)$. The decoder computes the syndromes $S_{j} = r(\alpha^{j})$,
//! the error-locator polynomial by the Berlekamp–Massey algorithm and its roots by a
//! Chien search over the $n$ positions; it corrects up to $t$ errors and reports a
//! failure beyond (when the locator has fewer roots than its degree). The field
//! arithmetic is that of [`comnum_math::gf2m`].
//!
//! # Examples
//! ```
//! use comnum_fec::bch::Bch;
//! // the (15, 7) BCH code, t = 2, of GF(16) with p(x) = 1 + x + x⁴
//! let code = Bch::new(4, 0b1_0011, 2, 15)?;
//! assert_eq!(code.k(), 7);
//! assert_eq!(code.generator(), &[1, 0, 0, 0, 1, 0, 1, 1, 1]);   // 1 + x⁴ + x⁶ + x⁷ + x⁸
//! let msg = [1, 0, 1, 1, 0, 0, 1];
//! let mut c = code.encode(&msg)?;
//! c[2] ^= 1;
//! c[11] ^= 1;
//! let (decoded, corrected) = code.decode(&c)?;
//! assert_eq!((decoded, corrected), (msg.to_vec(), Some(2)));
//! # Ok::<(), comnum_core::Error>(())
//! ```
//!
//! # References
//!
//! 1. S. Lin and D. J. Costello, *Error Control Coding*, 2nd ed., Pearson, 2004, ch. 6.
//! 2. ETSI EN 302 307-1 V1.4.1, *Digital Video Broadcasting (DVB); Second generation
//!    framing structure, channel coding and modulation systems (DVB-S2)*, 2014,
//!    § 5.3.1.

use crate::algebraic::{berlekamp_massey, chien_search};
use crate::block::{BlockCode, container_traits, expect_len};
use comnum_core::signal::hard_decisions;
use comnum_core::{Result, ensure};
use comnum_math::gf2::{poly_mul, poly_rem};
use comnum_math::gf2m::Field;

/// A shortened binary BCH code (see the module).
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Bch {
    field: Field,
    t: usize,
    n: usize,
    /// $g(x)$, lowest degree first.
    generator: Vec<u8>,
}

impl Bch {
    /// The code of $\mathrm{GF}(2^{m})$ with the primitive polynomial `poly` (bit $i$ the
    /// coefficient of $x^{i}$), correcting `t` errors, shortened to `n` bits.
    ///
    /// # Errors
    /// [`comnum_core::Error::Config`] for a polynomial that is not primitive, or a length
    /// that leaves no message bit or exceeds $2^{m} - 1$.
    pub fn new(m: u32, poly: u32, t: usize, n: usize) -> Result<Bch> {
        let field = Field::new(m, poly)?;
        ensure(t >= 1, || "BCH: t must be at least 1".into())?;
        let mut generator = vec![1u8];
        let mut used: Vec<Vec<u8>> = Vec::new();
        for i in (1..2 * t).step_by(2) {
            let mi = field.minimal_polynomial(i);
            if !used.contains(&mi) {
                generator = poly_mul(&generator, &mi);
                used.push(mi);
            }
        }
        let r = generator.len() - 1;
        ensure(n <= field.order() && n > r, || {
            format!("BCH: the length {n} must exceed the {r} parity bits and not exceed {}", field.order())
        })?;
        Ok(Bch { field, t, n, generator })
    }

    /// Codeword length $n$.
    pub fn n(&self) -> usize {
        self.n
    }

    /// Message length $k = n - \deg g$.
    pub fn k(&self) -> usize {
        self.n + 1 - self.generator.len()
    }

    /// Errors corrected $t$.
    pub fn t(&self) -> usize {
        self.t
    }

    /// The generator polynomial $g(x)$, lowest degree first.
    pub fn generator(&self) -> &[u8] {
        &self.generator
    }

    /// The codeword of $k$ message bits: the message, then the parity bits (highest
    /// degree first).
    ///
    /// # Errors
    /// [`comnum_core::Error::Shape`] for another number of bits.
    pub fn encode(&self, msg: &[u8]) -> Result<Vec<u8>> {
        expect_len(msg.len(), self.k(), "bits")?;
        // the parity bits: the remainder of m(x) x^r by g(x)
        let mut c: Vec<u8> = msg.iter().map(|b| b & 1).collect();
        c.extend(poly_rem(msg, &self.generator));
        Ok(c)
    }

    /// Decode $n$ hard bits: the message and the number of corrected errors (`None` if
    /// the word has more errors than the code corrects and they were detected; the
    /// message bits are then returned as received).
    ///
    /// # Errors
    /// [`comnum_core::Error::Shape`] for another number of bits.
    pub fn decode(&self, r: &[u8]) -> Result<(Vec<u8>, Option<usize>)> {
        expect_len(r.len(), self.n, "bits")?;
        let f = &self.field;
        let n = self.n;
        let two_t = 2 * self.t;
        // syndromes S_j = r(α^j), position i carrying x^{n-1-i}
        let mut s = vec![0u32; two_t + 1];
        for (i, &b) in r.iter().enumerate() {
            if b & 1 == 1 {
                let e = n - 1 - i;
                for (j, sj) in s.iter_mut().enumerate().skip(1) {
                    *sj ^= f.pow(j * e);
                }
            }
        }
        let k = self.k();
        let msg = r[..k].iter().map(|b| b & 1).collect::<Vec<u8>>();
        if s.iter().all(|&v| v == 0) {
            return Ok((msg, Some(0)));
        }
        // Berlekamp–Massey on S_1 … S_{2t}; at most t errors
        let (lambda, _) = berlekamp_massey(f, &s[1..]);
        let degree = lambda.len() - 1;
        if degree > self.t {
            return Ok((msg, None));
        }
        // Chien search: position i is in error if Λ(α^{-(n-1-i)}) = 0; as many roots as
        // the degree of Λ
        let order = f.order();
        let roots = chien_search(f, &lambda, n, |i| (n - 1 - i) % order);
        if roots.len() != degree {
            return Ok((msg, None));
        }
        let mut out: Vec<u8> = r.iter().map(|b| b & 1).collect();
        for &i in &roots {
            out[i] ^= 1;
        }
        let found = roots.len();
        out.truncate(k);
        Ok((out, Some(found)))
    }
}

impl BlockCode for Bch {
    fn k(&self) -> usize {
        Bch::k(self)
    }

    fn n(&self) -> usize {
        Bch::n(self)
    }

    fn encode_block(&self, info: &[u8]) -> Result<Vec<u8>> {
        self.encode(info)
    }

    /// Hard decisions of the LLRs, then [`Bch::decode`]; a decoding failure returns the
    /// message bits as received (as [`crate::rs::ReedSolomon`]).
    fn decode_block(&self, llr: &[f64]) -> Result<Vec<u8>> {
        expect_len(llr.len(), Bch::n(self), "LLRs")?;
        Ok(self.decode(&hard_decisions(llr))?.0)
    }
}

container_traits!(Bch);

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

    /// The encoder of the BCH code before `poly_rem`: its own register of 64-bit words.
    fn former_encode(code: &Bch, msg: &[u8]) -> Vec<u8> {
        let r = code.generator.len() - 1;
        let words = r.div_ceil(64);
        let mut g = vec![0u64; words];
        for (d, &b) in code.generator[..r].iter().enumerate() {
            if b == 1 {
                g[d / 64] |= 1 << (d % 64);
            }
        }
        let mut reg = vec![0u64; words];
        let top = r - 1;
        for &b in msg {
            let fb = ((reg[top / 64] >> (top % 64)) & 1) as u8 ^ (b & 1);
            for w in (0..words).rev() {
                let carry = if w > 0 { reg[w - 1] >> 63 } else { 0 };
                reg[w] = (reg[w] << 1) | carry;
            }
            if r % 64 != 0 {
                reg[words - 1] &= (1u64 << (r % 64)) - 1;
            }
            if fb == 1 {
                for (x, y) in reg.iter_mut().zip(&g) {
                    *x ^= y;
                }
            }
        }
        let mut c: Vec<u8> = msg.iter().map(|b| b & 1).collect();
        c.extend((0..r).rev().map(|d| ((reg[d / 64] >> (d % 64)) & 1) as u8));
        c
    }

    #[test]
    fn the_codewords_are_those_of_the_former_encoder() {
        let mut rng = Rng::new(17);
        // degrees 10 to 192: one word, several words, a multiple of 64
        for (m, poly, t, n) in
            [(5, 0b10_0101, 2, 31), (16, 0x1_002D, 12, 16_200), (16, 0x1_002D, 4, 1000), (8, 0x11D, 8, 255)]
        {
            let code = Bch::new(m, poly, t, n).unwrap();
            for _ in 0..3 {
                let msg: Vec<u8> = (0..code.k()).map(|_| rng.bit()).collect();
                assert_eq!(code.encode(&msg).unwrap(), former_encode(&code, &msg), "m {m}, t {t}, n {n}");
            }
        }
    }

    #[test]
    fn the_generators_of_dvb_s2_are_the_products_of_the_published_polynomials() {
        // ETSI EN 302 307-1, Table 6a: g1 … g12 of the normal FECFRAME, GF(2^16)
        let table6a: [&[usize]; 12] = [
            &[0, 2, 3, 5, 16],
            &[0, 1, 4, 5, 6, 8, 16],
            &[0, 2, 3, 4, 5, 7, 8, 9, 10, 11, 16],
            &[0, 2, 4, 6, 9, 11, 12, 14, 16],
            &[0, 1, 2, 3, 5, 8, 9, 10, 11, 12, 16],
            &[0, 2, 4, 5, 7, 8, 9, 10, 12, 13, 14, 15, 16],
            &[0, 2, 5, 6, 8, 9, 10, 11, 13, 15, 16],
            &[0, 1, 2, 5, 6, 8, 9, 12, 13, 14, 16],
            &[0, 5, 7, 9, 10, 11, 16],
            &[0, 1, 2, 5, 7, 8, 10, 12, 13, 14, 16],
            &[0, 2, 3, 5, 9, 11, 12, 13, 16],
            &[0, 1, 5, 6, 7, 9, 11, 12, 16],
        ];
        let f = Field::new(16, 0b1_0000_0000_0010_1101).unwrap();
        for (i, g) in table6a.iter().enumerate() {
            let mut p = vec![0u8; 17];
            for &d in *g {
                p[d] = 1;
            }
            assert_eq!(f.minimal_polynomial(2 * i + 1), p, "g{}", i + 1);
        }
    }

    #[test]
    fn up_to_t_errors_are_corrected_and_more_are_mostly_detected() {
        // the short DVB-S2 code of rate 1/2: GF(2^14), t = 12, (7200, 7032)
        let code = Bch::new(14, 0b100_0000_0010_1011, 12, 7200).unwrap();
        assert_eq!(code.k(), 7032);
        let mut rng = Rng::new(4);
        // `errors` distinct positions drawn at random
        let mut positions: Vec<usize> = (0..code.n()).collect();
        let (mut detected, mut wrong) = (0, 0);
        for errors in (0..=12).chain((13..=24).flat_map(|e| [e; 4])) {
            let msg: Vec<u8> = (0..code.k()).map(|_| rng.bit()).collect();
            let mut c = code.encode(&msg).unwrap();
            rng.shuffle(&mut positions);
            for &p in &positions[..errors] {
                c[p] ^= 1;
            }
            let (m, n) = code.decode(&c).unwrap();
            if errors <= 12 {
                // every pattern of at most t errors is corrected and counted
                assert_eq!((&m, n), (&msg, Some(errors)), "{errors} errors");
            } else if n.is_none() {
                // detected: the message bits come back as received
                assert_eq!(m, c[..code.k()], "{errors} errors");
                detected += 1;
            } else {
                // a decoded word within t of the received one cannot be the sent one
                assert_ne!(m, msg, "{errors} errors");
                wrong += 1;
            }
        }
        // a miscorrection needs the t roots of a wrong locator among the 7200 positions of
        // the shortened code out of 2^14 - 1: rare, none in these 48 patterns
        assert_eq!((detected, wrong), (48, 0), "{detected} detected, {wrong} miscorrected");
        // exactly t distinct errors
        let msg: Vec<u8> = (0..code.k()).map(|_| rng.bit()).collect();
        let mut c = code.encode(&msg).unwrap();
        for p in [3, 100, 2000, 5000, 7031, 7032, 7100, 7199, 1, 2, 4000, 6000] {
            c[p] ^= 1;
        }
        assert_eq!(code.decode(&c).unwrap(), (msg, Some(12)));
    }
}