diff --git a/Project.toml b/Project.toml index 5312a5b60..4296812b8 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/references.bib b/docs/src/references.bib index 09820c204..f67523219 100644 --- a/docs/src/references.bib +++ b/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} +} \ No newline at end of file diff --git a/docs/src/references.md b/docs/src/references.md index eaaf7cb51..f4d3fc0b2 100644 --- a/docs/src/references.md +++ b/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 diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl index fa56c7510..0e6216eb7 100644 --- a/src/ecc/ECC.jl +++ b/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 diff --git a/src/ecc/codes/classical/bch.jl b/src/ecc/codes/classical/bch.jl new file mode 100644 index 000000000..70580a909 --- /dev/null +++ b/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))) diff --git a/test/runtests.jl b/test/runtests.jl index b4846e7cb..1177d673f 100644 --- a/test/runtests.jl +++ b/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" diff --git a/test/test_ecc_bch.jl b/test/test_ecc_bch.jl new file mode 100644 index 000000000..2885cd2ae --- /dev/null +++ b/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 diff --git a/test/test_ecc_throws.jl b/test/test_ecc_throws.jl index 82e02303b..19468e1ad 100644 --- a/test/test_ecc_throws.jl +++ b/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)