Adding classical Bose–Chaudhuri–Hocquenghem code to ECC module (#263)

Co-authored-by: Stefan Krastanov <[email protected]>

Author: Fe-r-oz

Project.toml

@@ -53,7 +53,7 @@ LDPCDecoders = "0.3.1"
 LinearAlgebra = "1.9"
 MacroTools = "0.5.9"
 Makie = "0.20"
-Nemo = "0.38, 0.39, 0.40, 0.41, 0.42, 0.43, 0.44, 0.45"
+Nemo = "0.42, 0.43, 0.44, 0.45"
 Plots = "1.38.0"
 PrecompileTools = "1.2"
 PyQDecoders = "0.2.0"

docs/src/references.bib

@@ -341,3 +341,41 @@ @book{djordjevic2021quantum
   year={2021},
   publisher={Academic Press}
 }
+
+@article{hocquenghem1959codes,
+  title={Codes correcteurs d'erreurs},
+  author={Hocquenghem, Alexis},
+  journal={Chiffers},
+  volume={2},
+  pages={147--156},
+  year={1959}
+}
+
+@article{bose1960class,
+  title={On a class of error correcting binary group codes},
+  author={Bose, Raj Chandra and Ray-Chaudhuri, Dwijendra K},
+  journal={Information and control},
+  volume={3},
+  number={1},
+  pages={68--79},
+  year={1960},
+  publisher={Elsevier}
+}
+
+@article{bose1960further,
+  title={Further results on error correcting binary group codes},
+  author={Bose, Raj Chandra and Ray-Chaudhuri, Dwijendra K},
+  journal={Information and Control},
+  volume={3},
+  number={3},
+  pages={279--290},
+  year={1960},
+  publisher={Elsevier}
+}
+
+@book{error2024lin,
+  title={Error Control Coding},
+  author={Lin, Shu and Costello, Daniel},
+  year={2024},
+  publisher={Pearson}
+}

docs/src/references.md

@@ -43,6 +43,10 @@ For classical code construction routines:
 - [raaphorst2003reed](@cite)
 - [abbe2020reed](@cite)
 - [djordjevic2021quantum](@cite)
+- [hocquenghem1959codes](@cite)
+- [bose1960class](@cite)
+- [bose1960further](@cite)
+- [error2024lin](@cite)
 
 # References
 

src/ecc/ECC.jl

@@ -8,7 +8,7 @@ using DocStringExtensions
 using Combinatorics: combinations
 using SparseArrays: sparse
 using Statistics: std
-using Nemo: ZZ, residue_ring, matrix
+using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible 
 
 abstract type AbstractECC end
 
@@ -57,6 +57,13 @@ function iscss(c::AbstractECC)
     return iscss(typeof(c))
 end
 
+"""
+Generator Polynomial `g(x)`
+
+In a [polynomial code](https://en.wikipedia.org/wiki/Polynomial_code), the generator polynomial `g(x)` is a polynomial of the minimal degree over a finite field `F`. The set of valid codewords in the code consists of all polynomials that are divisible by `g(x)` without remainder.
+"""
+function generator_polynomial end
+
 parity_checks(s::Stabilizer) = s
 Stabilizer(c::AbstractECC) = parity_checks(c)
 MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); kwarg...)
@@ -347,4 +354,5 @@ include("codes/toric.jl")
 include("codes/gottesman.jl")
 include("codes/surface.jl")
 include("codes/classical/reedmuller.jl")
+include("codes/classical/bch.jl")
 end #module

src/ecc/codes/classical/bch.jl

@@ -0,0 +1,104 @@
+"""The family of Bose–Chaudhuri–Hocquenghem (BCH) codes, as discovered in 1959 by Alexis Hocquenghem [hocquenghem1959codes](@cite), and independently in 1960 by Raj Chandra Bose and D.K. Ray-Chaudhuri [bose1960class](@cite).
+
+The binary parity check matrix can be obtained from the following matrix over GF(2) field elements:
+
+```
+1		(α¹)¹			(α¹)²			(α¹)³			...		(α¹)ⁿ ⁻ ¹
+1		(α³)¹			(α³)²			(α³)³			...		(α³)ⁿ ⁻ ¹
+1		(α⁵)¹			(α⁵)²			(α⁵)³			...		(α⁵)ⁿ ⁻ ¹
+.		  .			  .			  .			...		   .
+.		  .			  .			  .			...		   .
+.		  .			  .			  .			...		   .
+1		(α²ᵗ ⁻ ¹)¹		(α²ᵗ ⁻ ¹)²		(α²ᵗ ⁻ ¹)³		...		(α²ᵗ ⁻ ¹)ⁿ ⁻ ¹
+```
+
+The entries of the matrix are in GF(2ᵐ). Each element in GF(2ᵐ) can be represented by an `m`-tuple (a binary column vector of length `m`). If each entry of `H` is replaced by its corresponding `m`-tuple, we obtain a binary parity check matrix for the code.
+
+The BCH code is cyclic as its generator polynomial, `g(x)` divides `xⁿ - 1`, so `mod (xⁿ - 1, g(x)) = 0`.
+
+You might be interested in consulting [bose1960further](@cite) and [error2024lin](@cite) as well.
+
+The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/q-ary_bch).
+"""
+
+abstract type AbstractPolynomialCode <: ClassicalCode end
+
+"""
+`BCH(m, t)`
+- `m`: The positive integer defining the degree of the finite (Galois) field, GF(2ᵐ).
+- `t`: The positive integer specifying the number of correctable errors.
+"""
+struct BCH <: AbstractPolynomialCode
+    m::Int 
+    t::Int 
+    function BCH(m, t)
+        if m < 3 || t < 0 || t >= 2 ^ (m - 1)
+            throw(ArgumentError("Invalid parameters: `m` and `t` must be positive. Additionally, ensure `m ≥ 3` and `t < 2ᵐ ⁻ ¹` to obtain a valid code."))
+        end
+        new(m, t)
+    end
+end
+
+"""
+Generator Polynomial of BCH Codes
+
+This function calculates the generator polynomial `g(x)` of a `t`-bit error-correcting BCH code of binary length `n = 2ᵐ - 1`. The binary code is derived from a code over the finite Galois field GF(2).
+
+`generator_polynomial(BCH(m, t))`
+
+- `m`: The positive integer defining the degree of the finite (Galois) field, GF(2ᵐ).
+- `t`: The positive integer specifying the number of correctable errors.
+
+The generator polynomial `g(x)` is the fundamental polynomial used for encoding and decoding BCH codes. It has the following properties:
+
+1. Roots: It has `α`, `α²`, `α³`, ..., `α²ᵗ` as its roots, where `α` is a primitive element of the Galois Field GF(2ᵐ).
+2. Error Correction: A BCH code with generator polynomial `g(x)` can correct up to `t` errors in a codeword of length `2ᵐ - 1`.
+3. Minimal Polynomials: `g(x)` is the least common multiple (LCM) of the minimal polynomials `φᵢ(x)` of `αⁱ` for `i` from `1` to `2ᵗ`.
+
+Useful definitions and background:
+
+Minimal Polynomial: The minimal polynomial of a field element `α` in GF(2ᵐ) is the polynomial of the lowest degree over GF(2) that has `α` as a root.
+
+Least Common Multiple (LCM): The LCM of two or more polynomials `fᵢ(x)` is the polynomial with the lowest degree that is a multiple of all `fᵢ(x)`. It ensures that `g(x)` has all the roots of `φᵢ(x)` for `i = 1` to `2ᵗ`.
+
+Conway polynomial: The finite Galois field `GF(2ᵐ)` can have multiple distinct primitive polynomials of the same degree due to existence of several irreducible polynomials of that degree, each generating the field through different roots. Nemo.jl uses [Conway polynomial](https://en.wikipedia.org/wiki/Conway_polynomial_(finite_fields)), a standard way to represent the primitive polynomial for finite Galois fields `GF(pᵐ)` of degree `m`, where `p` is a prime number.
+"""
+function generator_polynomial(b::BCH)
+    GF2ʳ, a = finite_field(2, b.m, "a")
+    GF2x, x = GF(2)["x"]
+    minimal_poly = FqPolyRingElem[]
+    for i in 1:2 * b.t
+        if i % 2 != 0
+            push!(minimal_poly, minpoly(GF2x, a ^ i))
+        end 
+    end
+    gx = lcm(minimal_poly)
+    return gx
+end
+
+function parity_checks(b::BCH)
+    GF2ʳ, a = finite_field(2, b.m, "a")
+    HField = Matrix{FqFieldElem}(undef, b.t, 2 ^ b.m - 1)
+    for i in 1:b.t
+        for j in 1:2 ^ b.m - 1
+            base = 2 * i - 1  
+            HField[i, j] = (a ^ base) ^ (j - 1)
+        end
+    end
+    H = Matrix{Bool}(undef, b.m * b.t, 2 ^ b.m - 1)
+    for i in 1:b.t
+        row_start = (i - 1) * b.m + 1
+        row_end = row_start + b.m - 1
+        for j in 1:2 ^ b.m - 1
+            t_tuple = Bool[]
+            for k in 0:b.m - 1
+                push!(t_tuple, !is_zero(coeff(HField[i, j], k)))
+            end 
+            H[row_start:row_end, j] .=  vec(t_tuple')
+        end
+    end 
+    return H
+end
+
+code_n(b::BCH) = 2 ^ b.m - 1
+code_k(b::BCH) = 2 ^ b.m - 1 - degree(generator_polynomial(BCH(b.m, b.t)))

test/runtests.jl

@@ -64,6 +64,7 @@ end
 @doset "ecc_encoding"
 @doset "ecc_gottesman"
 @doset "ecc_reedmuller"
+@doset "ecc_bch"
 @doset "ecc_syndromes"
 @doset "ecc_throws"
 @doset "precompile"

test/test_ecc_bch.jl

@@ -0,0 +1,85 @@
+using Test
+using LinearAlgebra
+using QuantumClifford
+using QuantumClifford.ECC
+using QuantumClifford.ECC: AbstractECC, BCH, generator_polynomial
+using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible 
+
+"""
+- To prove that `t`-bit error correcting BCH code indeed has minimum distance of at least `2 * t + 1`, it is shown that no `2 * t` or fewer columns of its binary parity check matrix `H` sum to zero. A formal mathematical proof can be found on pages 168 and 169 of Ch6 of Error Control Coding by Lin, Shu and Costello, Daniel. 
+- The parameter `2 * t + 1` is usually called the designed distance of the `t`-bit error correcting BCH code.
+"""
+function check_designed_distance(matrix, t)
+    n_cols = size(matrix, 2)
+    for num_cols in 1:2 * t
+        for i in 1:n_cols - num_cols + 1
+            combo = matrix[:, i:(i + num_cols - 1)]
+            sum_cols = sum(combo, dims = 2)
+            if all(sum_cols .== 0)
+                return false  # Minimum distance is not greater than `2 * t`.
+            end
+        end
+    end
+    return true  # Minimum distance is at least `2 * t + 1`.
+end
+
+@testset "Testing properties of BCH codes" begin
+    m_cases = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
+    for m in m_cases
+        n = 2 ^ m - 1
+        for t in rand(1:m, 2)
+            H = parity_checks(BCH(m, t))
+            @test check_designed_distance(H, t) == true
+            # n - k == degree of generator polynomial, `g(x)` == rank of binary parity check matrix, `H`.
+            mat = matrix(GF(2), parity_checks(BCH(m, t)))
+            computed_rank = rank(mat)
+            @test computed_rank == degree(generator_polynomial(BCH(m, t)))
+            @test code_k(BCH(m, t)) == n - degree(generator_polynomial(BCH(m, t)))
+            # BCH code is cyclic as its generator polynomial, `g(x)` divides `xⁿ - 1`, so `mod (xⁿ - 1, g(x))` = 0.
+            gx = generator_polynomial(BCH(m, t))
+            GF2x, x = GF(2)["x"] 
+            @test mod(x ^ n - 1, gx) == 0
+        end
+    end
+
+    #example taken from Ch6 of Error Control Coding by Lin, Shu and Costello, Daniel
+    @test parity_checks(BCH(4, 2))  ==    [1  0  0  0  1  0  0  1  1  0  1  0  1  1  1;
+                                           0  1  0  0  1  1  0  1  0  1  1  1  1  0  0;
+                                           0  0  1  0  0  1  1  0  1  0  1  1  1  1  0;
+                                           0  0  0  1  0  0  1  1  0  1  0  1  1  1  1;
+                                           1  0  0  0  1  1  0  0  0  1  1  0  0  0  1;
+                                           0  0  0  1  1  0  0  0  1  1  0  0  0  1  1;
+                                           0  0  1  0  1  0  0  1  0  1  0  0  1  0  1;
+                                           0  1  1  1  1  0  1  1  1  1  0  1  1  1  1]
+    
+    # Examples taken from https://web.ntpu.edu.tw/~yshan/BCH_code.pdf.
+    GF2x, x = GF(2)["x"]
+    GF2⁴, a = finite_field(2, 4, "a")
+    GF2⁶, b = finite_field(2, 6, "b")
+    @test defining_polynomial(GF2x, GF2⁴) ==  x ^ 4 + x + 1
+    @test is_irreducible(defining_polynomial(GF2x, GF2⁴)) == true
+    @test generator_polynomial(BCH(4, 2)) == x ^ 8 + x ^ 7 +  x ^ 6 +  x ^ 4 + 1
+    @test generator_polynomial(BCH(4, 3)) == x ^ 10 + x ^ 8 +  x ^ 5 +  x ^ 4 +  x ^ 2 + x + 1
+ 
+    # Nemo.jl uses [Conway polynomial](https://en.wikipedia.org/wiki/Conway_polynomial_(finite_fields)), a standard way to represent the primitive polynomial for finite Galois fields `GF(pᵐ)` of degree `m`, where `p` is a prime number. 
+    # The `GF(2⁶)`'s Conway polynomial is `p(z) = z⁶ + z⁴ + z³ + z + 1`. In contrast, the polynomial given in https://web.ntpu.edu.tw/~yshan/BCH_code.pdf is `p(z) = z⁶ + z + 1`. Because both polynomials are irreducible, they are also primitive polynomials for `GF(2⁶)`.
+    
+    test_cases = [(6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6), (6, 7), (6, 10), (6, 11), (6, 13), (6, 15)]
+    @test defining_polynomial(GF2x, GF2⁶) == x ^ 6 + x ^ 4 + x ^ 3 + x + 1
+    @test is_irreducible(defining_polynomial(GF2x, GF2⁶)) == true    
+    for i in 1:length(test_cases)
+        m, t = test_cases[i]
+        if t == 1
+            @test generator_polynomial(BCH(m, t)) == defining_polynomial(GF2x, GF2⁶)
+        else
+            prev_t = test_cases[i - 1][2]
+            @test generator_polynomial(BCH(m, t)) == generator_polynomial(BCH(m, prev_t)) * minpoly(GF2x, b ^ (t + prev_t - (t - prev_t - 1)))
+        end
+    end
+    
+    results = [57 51 45 39 36 30 24 18 16 10 7]
+    for (result, (m, t)) in zip(results, test_cases)
+        @test code_k(BCH(m, t)) == result
+        @test check_designed_distance(parity_checks(BCH(m, t)), t) == true
+    end
+end

test/test_ecc_throws.jl

@@ -1,7 +1,9 @@
 using Test
 using QuantumClifford
-using QuantumClifford.ECC
-using QuantumClifford.ECC: ReedMuller
+using QuantumClifford.ECC: ReedMuller, BCH
 
 @test_throws ArgumentError ReedMuller(-1, 3)
 @test_throws ArgumentError ReedMuller(1, 0)
+
+@test_throws ArgumentError BCH(2, 2)
+@test_throws ArgumentError BCH(3, 4)