Binary Golay codes: [24, 12, 8] and [23, 12, 7] (#276)

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

docs/src/references.bib

@@ -419,13 +419,13 @@ @article{RevModPhys.87.307
 }
 
 @misc{goodenough2024bipartiteentanglementnoisystabilizer,
-  title={Bipartite entanglement of noisy stabilizer states through the lens of stabilizer codes}, 
+  title={Bipartite entanglement of noisy stabilizer states through the lens of stabilizer codes},
   author={Kenneth Goodenough and Aqil Sajjad and Eneet Kaur and Saikat Guha and Don Towsley},
   year={2024},
   eprint={2406.02427},
   archivePrefix={arXiv},
   primaryClass={quant-ph},
-  url={https://arxiv.org/abs/2406.02427}, 
+  url={https://arxiv.org/abs/2406.02427},
 }
 
 @article{panteleev2021degenerate,
@@ -532,7 +532,7 @@ @article{wang2024coprime
 }
 
 @misc{voss2024multivariatebicyclecodes,
-  title={Multivariate Bicycle Codes}, 
+  title={Multivariate Bicycle Codes},
   author={Lukas Voss and Sim Jian Xian and Tobias Haug and Kishor Bharti},
   year={2024},
   eprint={2406.19151},
@@ -561,3 +561,26 @@ @article{haah2011local
   pages={042330},
   year={2011},
 }
+
+@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

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

src/ecc/ECC.jl

@@ -384,10 +384,11 @@ include("codes/gottesman.jl")
 include("codes/surface.jl")
 include("codes/concat.jl")
 include("codes/random_circuit.jl")
+include("codes/quantumreedmuller.jl")
 include("codes/classical/reedmuller.jl")
 include("codes/classical/recursivereedmuller.jl")
 include("codes/classical/bch.jl")
-include("codes/quantumreedmuller.jl")
+include("codes/classical/golay.jl")
 
 # qLDPC
 include("codes/classical/lifted.jl")

src/ecc/codes/classical/golay.jl

@@ -0,0 +1,86 @@
+"""
+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:
+
+- Binary `[23, 12, 7]` Golay code: The perfect code with code length `n = 23`
+and dimension `k = 12`. By puncturing in any of the coordinates of parity check
+matrix `H` = `[24, 12, 8]`, we obtain a `[23, 12, 7]` Golay code.
+
+- 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₂₄'.
+
+The parity check matrix is defined as follows: `H₂₄ = [I₁₂ | A']` where `I₁₂` is the
+`12 × 12` identity matrix and `A` is a bordered reverse circulant matrix. 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
+
+# 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(A[i - 1, 2:end], -1)
+        A[i, 1] = 1
+    end
+    return A
+end
+
+function generator(g::Golay)
+    if g.n == 24
+        A₂₄ = _create_A₂₄_golay(24)
+        I₁₂ = LinearAlgebra.Diagonal(ones(Int, g.n ÷ 2))
+        G₂₄ = hcat(I₁₂, (A₂₄)')
+        return G₂₄
+    else
+        A₂₄ = _create_A₂₄_golay(24)
+        A₂₃ = A₂₄[:, 1:end - 1]
+        I₁₂ = LinearAlgebra.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₁₂ = LinearAlgebra.Diagonal(ones(Int, g.n ÷ 2))
+        H₂₄ = hcat((A₂₄)', I₁₂)
+        return H₂₄
+    else
+        A₂₄ = _create_A₂₄_golay(24)
+        A₂₃ = A₂₄[:, 1:end - 1]
+        I₁₂ = LinearAlgebra.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/Project.toml

@@ -27,5 +27,6 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
+StyledStrings = "f489334b-da3d-4c2e-b8f0-e476e12c162b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"

test/test_ecc_golay.jl

@@ -0,0 +1,126 @@
+@testitem "ECC Golay" begin
+
+    using LinearAlgebra
+    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.
+    # Thus, 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
+
+    function minimum_distance(H)
+        n = size(H, 2)
+        min_dist = n + 1
+        for x_bits in 1:(2^n - 1)
+            x = reverse(digits(x_bits, base=2, pad=n))
+            xᵀ = reshape(x, :, 1)
+            if all(mod.(H * xᵀ, 2) .== 0)
+                weight = sum(x)
+                min_dist = min(min_dist, weight)
+            end
+        end
+        return min_dist
+    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]
+
+        # minimum distance test
+        # [24, 12, 8]
+        H = parity_checks(Golay(24))
+        @test minimum_distance(H) == 8
+        # [23, 12, 7]
+        H = parity_checks(Golay(23))
+        @test minimum_distance(H) == 7
+
+        # 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]
+
+        H = parity_checks(Golay(24))
+        @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 extended_mat == H
+        end
+    end
+end

test/test_ecc_throws.jl

@@ -1,6 +1,6 @@
 @testitem "ECC throws" begin
 
-    using QuantumClifford.ECC: ReedMuller, BCH, RecursiveReedMuller
+    using QuantumClifford.ECC: ReedMuller, BCH, RecursiveReedMuller, Golay
 
     @test_throws ArgumentError ReedMuller(-1, 3)
     @test_throws ArgumentError ReedMuller(1, 0)
@@ -12,4 +12,6 @@
     @test_throws ArgumentError RecursiveReedMuller(-1, 3)
     @test_throws ArgumentError RecursiveReedMuller(1, 0)
     @test_throws ArgumentError RecursiveReedMuller(4, 2)
+
+    @test_throws ArgumentError Golay(21)
 end

test/test_jet.jl

@@ -10,6 +10,7 @@
     using AbstractAlgebra
     using Hecke
     using StaticArrays
+    using StyledStrings
 
     rep = report_package("QuantumClifford";
         ignored_modules=(
@@ -23,6 +24,7 @@
             AnyFrameModule(AbstractAlgebra),
             AnyFrameModule(Hecke),
             AnyFrameModule(StaticArrays),
+            AnyFrameModule(StyledStrings),
     ))
 
     @show rep