Binary Golay codes: [24, 12, 8] and [23, 12, 7]

CHANGELOG.md

@@ -7,6 +7,7 @@
 
 ## v0.9.4 - 2024-06-28
 
+- Addition of classical Golay codes to ECC module
 - Addition of a constructor for concatenated quantum codes -- `Concat`.
 - Addition of multiple unexported classical code constructors.
 - Failed compactification of gates now only raises a warning instead of throwing an error. Defaults to slower non-compactified gates.

docs/src/references.bib

@@ -386,3 +386,27 @@ @article{knill1996concatenated
   journal={arXiv preprint quant-ph/9608012},
   year={1996}
 }
+
+@article{golay1949notes,
+  title={Notes on digital coding},
+  author={Golay, Marcel JE},
+  journal={Proc. IEEE},
+  volume={37},
+  pages={657},
+  year={1949}
+}
+
+@book{huffman2010fundamentals,
+  title={Fundamentals of error-correcting codes},
+  author={Huffman, W Cary and Pless, Vera},
+  year={2010},
+  publisher={Cambridge university press}
+
+}
+
+@article{bhatia2018mceliece,
+  title={McEliece cryptosystem based on extended Golay code},
+  author={Bhatia, Amandeep Singh and Kumar, Ajay},
+  journal={arXiv preprint arXiv:1811.06246},
+  year={2018}
+}

docs/src/references.md

@@ -48,6 +48,9 @@ For classical code construction routines:
 - [bose1960class](@cite)
 - [bose1960further](@cite)
 - [error2024lin](@cite)
+- [golay1949notes](@cite)
+- [huffman2010fundamentals](@cite)
+- [bhatia2018mceliece](@cite)
 
 # References
 

src/ecc/ECC.jl

@@ -69,6 +69,8 @@ In a [polynomial code](https://en.wikipedia.org/wiki/Polynomial_code), the gener
 """
 function generator_polynomial end
 
+function generator end
+
 parity_checks(s::Stabilizer) = s
 Stabilizer(c::AbstractECC) = parity_checks(c)
 MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); kwarg...)
@@ -361,4 +363,5 @@ include("codes/surface.jl")
 include("codes/concat.jl")
 include("codes/classical/reedmuller.jl")
 include("codes/classical/bch.jl")
+include("codes/classical/golay.jl")
 end #module

src/ecc/codes/classical/golay.jl

@@ -0,0 +1,84 @@
+"""
+The family of classical binary Golay codes were discovered by Edouard Golay in his 1949 paper [golay1949notes](@cite), where he described the binary `[23, 12, 7]` Golay code.
+
+There are two binary Golay codes:
+
+1. Binary `[23, 12, 7]` Golay code: The perfect code with code length `n` 23 and dimension `k` 12. Minimum distance is 7, implying it can detect or correct up to 3 errors. By puncturing in any of the coordinates of parity check Matrix `H` = `[24, 12, 8]`, we obtain a `[23, 12, 7]` Golay code.
+
+2. Extended Binary `[24, 12, 8]` Golay code: Obtained by adding a parity check bit to `[23, 12, 7]`. The bordered reverse circulant matrix `(A)` of `[24, 12, 8]` Golay code is self-dual, i.e., A₂₄ is same as A₂₄'. 
+
+Parity Check Matrix `(H)`: `H` is defined as follows: `H₂₄ = [I₁₂ | A']` where `I₁₂` is the 12 x 12 identity matrix and `A` is a bordered reverse circulant matrix.
+
+Construction method for `A` [huffman2010fundamentals](@cite): The columns of `A` are labeled by ∞, 0, 1, 2, ..., 10. The first row contains 0 in column ∞ and 1 elsewhere. To obtain the second row, a 1 is placed in column ∞ and a 1 is placed in columns 0, 1, 3, 4, 5, and 9; these numbers are precisely the squares of the integers modulo 11. That is, 0² = 0, 1² ≡ 10² ≡ 1 (mod 11), 2² ≡ 2² ≡ 4 (mod 11), etc. The third row of `A` is obtained by putting a 1 in column ∞ and then shifting the components in the second row one place to the left and wrapping the entry in column 0 around to column 10. The fourth row is obtained from the third in the same manner, as are the remaining rows.
+
+Puncturing and then extending any column in​ with an overall parity check `H₂₃` reconstructs the original parity check matrix `H₂₄`. Thus, all punctured codes are equivalent.
+
+The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/golay).
+"""
+struct Golay <: ClassicalCode
+    n::Int 
+
+    function Golay(n)
+        if !(n in (23, 24))
+            throw(ArgumentError("Invalid parameters: `n` must be either 24 or 23 to obtain a valid code."))
+        end
+        new(n)
+    end
+end
+
+function _circshift_row_golay(row::Vector{Int}, shift::Int, n::Int)
+    l = length(row)
+    return [row[mod((i - shift - 1), l) + 1] for i in 1:l]
+end
+
+# bordered reverse circulant matrix (see section 1.9.1, pg. 30-33) of [huffman2010fundamentals](@cite).
+function _create_A₂₄_golay(n::Int)
+    A = zeros(Int, n ÷ 2, n ÷ 2)
+    # Define the squared values modulo 11.
+    squares_mod₁₁ = [0, 1, 4, 9, 5, 3, 3, 5, 9, 4, 1]
+    A[1, 2:end] .= 1
+    A[2, 1] = 1
+    for i in squares_mod₁₁
+        A[2, i + 2] = 1
+    end
+    # Fill in the rest of the rows using the reverse circulant property.
+    for i in 3:n ÷ 2
+        A[i, 2:end] = _circshift_row_golay(A[i - 1, 2:end], -1, n ÷ 2)
+        A[i, 1] = 1 
+    end
+    return A
+end
+
+function generator(g::Golay)
+    if g.n == 24 
+        A₂₄ = _create_A₂₄_golay(24)
+        I₁₂ = Diagonal(ones(Int, g.n ÷ 2))
+        G₂₄ = hcat(I₁₂, (A₂₄)')
+        return G₂₄
+    else 
+        A₂₄ = _create_A₂₄_golay(24)
+        A₂₃ = A₂₄[:, 1:end - 1]
+        I₁₂ = Diagonal(ones(Int, g.n ÷ 2))
+        G₂₃ = hcat(I₁₂, (A₂₃)')
+        return G₂₃
+    end
+end
+
+function parity_checks(g::Golay)
+    if g.n == 24 
+        A₂₄ = _create_A₂₄_golay(24)
+        I₁₂ = Diagonal(ones(Int, g.n ÷ 2))
+        H₂₄ = hcat((A₂₄)', I₁₂)
+        return H₂₄
+    else 
+        A₂₄ = _create_A₂₄_golay(24)
+        A₂₃ = A₂₄[:, 1:end - 1]
+        I₁₂ = Diagonal(ones(Int, g.n ÷ 2))
+        H₂₃ = hcat((A₂₃)', I₁₂)
+        return H₂₃
+    end
+end
+
+code_n(g::Golay) = g.n
+code_k(g::Golay) = 12
+distance(g::Golay) = code_n(g::Golay) - code_k(g::Golay) - 4

test/runtests.jl

@@ -65,6 +65,7 @@ end
 @doset "ecc_gottesman"
 @doset "ecc_reedmuller"
 @doset "ecc_bch"
+@doset "ecc_golay"
 @doset "ecc_syndromes"
 @doset "ecc_throws"
 @doset "precompile"

test/test_ecc_golay.jl

@@ -0,0 +1,99 @@
+using Test
+using LinearAlgebra
+using QuantumClifford
+using QuantumClifford.ECC
+using QuantumClifford.ECC: AbstractECC, Golay, generator
+using Nemo: matrix, GF, echelon_form
+
+"""
+- Theorem: Let `C` be a binary linear code. If `C` is self-orthogonal and has a generator matrix `G` where each row has weight divisible by four, then every codeword of `C` has weight 
+divisible by four. 
+- `H₂₄` is self-dual because its generator matrix has all rows with weight divisible by four. By above theorem, all codewords of `H₂₄` must have weights divisible by four. Refer to pg. 30 to 33
+of Ch1 of Fundamentals of Error Correcting Codes by Huffman, Cary and Pless, Vera.
+""" 
+function code_weight_property(matrix)
+    for row in eachrow(matrix)
+        count = sum(row)
+        if count % 4 == 0
+            return true
+        end
+    end
+    return false
+end
+
+"""
+Test the equivalence of punctured and extended code by verifying that puncturing the binary parity check matrix H₂₄ in any coordinate and then extending it by
+adding an overall parity check in the same position yields the original matrix H₂₄.
+
+Steps:
+1. Puncture the Code: Remove the i-th column from H₂₄ to create a punctured matrix H₂₃. Note: H₂₃ = H[:, [1:i-1; i+1:end]]
+2. Extend the Code: Add a column in the same position to ensure each row has even parity. Note: H'₂₄ = [H₂₃[:, 1:i-1] c H₂₃[:, i:end]]. Here, c is a column vector added to ensure 
+each row in H'₂₄ has even parity.
+3. Equivalence Check: Verify that H'₂₄ = H₂₄.
+"""
+function puncture_code(mat, i)
+    return mat[:, [1:i-1; i+1:end]]
+end
+
+function extend_code(mat, i)
+    k, _ = size(mat)
+    extended_mat = hcat(mat[:, 1:i-1], zeros(Bool, k, 1), mat[:, i:end])
+    # Calculate the parity for each row
+    for row in 1:k
+        row_parity = sum(mat[row, :]) % 2 == 1
+        extended_mat[row, i] = row_parity
+    end
+    return extended_mat
+end
+
+@testset "Testing binary Golay codes properties" begin
+    test_cases = [(24, 12), (23, 12)]
+    for (n, k) in test_cases
+        H = parity_checks(Golay(n))
+        mat = matrix(GF(2), parity_checks(Golay(n)))
+        computed_rank = rank(mat)
+        @test computed_rank == n - k
+    end
+    # [[24, 12, 8]] binary Golay code is a self-dual code [huffman2010fundamentals](@cite).
+    H = parity_checks(Golay(24))
+    @test code_weight_property(H) == true
+    @test H[:, (12 + 1):end]  == H[:, (12 + 1):end]'
+    # Example taken from [huffman2010fundamentals](@cite).
+    @test parity_checks(Golay(24)) == [0  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0;
+                                       1  1  1  0  1  1  1  0  0  0  1  0  0  1  0  0  0  0  0  0  0  0  0  0;
+                                       1  1  0  1  1  1  0  0  0  1  0  1  0  0  1  0  0  0  0  0  0  0  0  0;
+                                       1  0  1  1  1  0  0  0  1  0  1  1  0  0  0  1  0  0  0  0  0  0  0  0;
+                                       1  1  1  1  0  0  0  1  0  1  1  0  0  0  0  0  1  0  0  0  0  0  0  0;
+                                       1  1  1  0  0  0  1  0  1  1  0  1  0  0  0  0  0  1  0  0  0  0  0  0;
+                                       1  1  0  0  0  1  0  1  1  0  1  1  0  0  0  0  0  0  1  0  0  0  0  0;
+                                       1  0  0  0  1  0  1  1  0  1  1  1  0  0  0  0  0  0  0  1  0  0  0  0;
+                                       1  0  0  1  0  1  1  0  1  1  1  0  0  0  0  0  0  0  0  0  1  0  0  0;
+                                       1  0  1  0  1  1  0  1  1  1  0  0  0  0  0  0  0  0  0  0  0  1  0  0;
+                                       1  1  0  1  1  0  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0;
+                                       1  0  1  1  0  1  1  1  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1]
+
+    # cross-verifying the canonical equivalence of bordered reverse circulant matrix (A) from [huffman2010fundamentals](@cite) with matrix A taken from [bhatia2018mceliece](@cite).
+    A = [1 1 0 1 1 1 0 0 0 1 0 1;
+         1 0 1 1 1 0 0 0 1 0 1 1;
+         0 1 1 1 0 0 0 1 0 1 1 1;
+         1 1 1 0 0 0 1 0 1 1 0 1;
+         1 1 0 0 0 1 0 1 1 0 1 1;
+         1 0 0 0 1 0 1 1 0 1 1 1;
+         0 0 0 1 0 1 1 0 1 1 1 1;
+         0 0 1 0 1 1 0 1 1 1 0 1;
+         0 1 0 1 1 0 1 1 1 0 0 1;
+         1 0 1 1 0 1 1 1 0 0 0 1;
+         0 1 1 0 1 1 1 0 0 0 1 1;
+         1 1 1 1 1 1 1 1 1 1 1 0]
+     @test echelon_form(matrix(GF(2), A)) == echelon_form(matrix(GF(2), H[1:12, 1:12]))
+     # test self-duality for extended Golay code, G == H
+     @test echelon_form(matrix(GF(2), generator(Golay(24)))) == echelon_form(matrix(GF(2), parity_checks(Golay(24))))
+    # All punctured and extended matrices are equivalent to the H₂₄.
+    # Test each column for puncturing and extending
+    for i in 1:24
+        punctured_mat = puncture_code(H, i)
+        extended_mat = extend_code(punctured_mat, i)
+        # Test equivalence
+        @test extended_mat == H
+    end
+end

test/test_ecc_throws.jl

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