Bug fix to the parity_checks(ReedMuller(r, m)) along with `Recursi…

…veReedMuller` code for cross-reference (#277)


---------

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

CHANGELOG.md

@@ -5,6 +5,11 @@
 
 # News
 
+## dev
+
+- **(fix)** Bug fix to the `parity_checks(ReedMuller(r, m))` of classical Reed-Muller code (it was returning generator matrix).
+- `RecursiveReedMuller` code implementation as an alternative implementation of `ReedMuller`.
+
 ## v0.9.9 - 2024-08-05
 
 - `inv` is implemented for all Clifford operator types (symbolic, dense, sparse).

src/ecc/ECC.jl

@@ -1,14 +1,15 @@
 module ECC
 
 using LinearAlgebra
+using LinearAlgebra: I
 using QuantumClifford
 using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer
 import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
 using DocStringExtensions
 using Combinatorics: combinations
 using SparseArrays: sparse
 using Statistics: std
-using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible
+using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible, echelon_form
 
 abstract type AbstractECC end
 
@@ -70,6 +71,9 @@ In a [polynomial code](https://en.wikipedia.org/wiki/Polynomial_code), the gener
 """
 function generator_polynomial end
 
+"""The generator matrix of a code."""
+function generator end
+
 parity_checks(s::Stabilizer) = s
 Stabilizer(c::AbstractECC) = parity_checks(c)
 MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); kwarg...)
@@ -363,4 +367,5 @@ include("codes/concat.jl")
 include("codes/random_circuit.jl")
 include("codes/classical/reedmuller.jl")
 include("codes/classical/bch.jl")
+include("codes/classical/recursivereedmuller.jl")
 end #module

src/ecc/codes/classical/recursivereedmuller.jl

@@ -0,0 +1,64 @@
+"""
+A construction of the Reed-Muller class of codes using the recursive definition.
+
+The Plotkin `(u, u + v)` construction defines a recursive relation between generator matrices of Reed-Muller `(RM)` codes [abbe2020reed](@cite). To derive the generator matrix `G(m, r)` for `RM(r, m)`, the generator matrices of lower-order codes are utilized:
+- `G(r - 1, m - 1)`: Generator matrix of `RM(r - 1, m - 1)`
+- `G(r, m - 1)`: Generator matrix of `RM(r, m - 1)`
+
+The generator matrix `G(m, r)` of `RM(m, r)` is formulated as follows in matrix notation:
+
+```math
+G(m, r) = \begin{bmatrix}
+G(r, m - 1) & G(r, m - 1) \\
+0 & G(r - 1, m - 1)
+\end{bmatrix}
+```
+
+Here, the matrix 0 denotes an all-zero matrix with dimensions matching `G(r - 1, m - 1)`. This recursive approach facilitates the construction of higher-order Reed-Muller codes based on the generator matrices of lower-order codes.
+
+In addition, the dimension of `RM(m - r - 1, m)` equals the dimension of the dual of `RM(r, m)`. Thus, `RM(m - r - 1, m) = RM(r, m)^⊥` shows that the [dual code](https://en.wikipedia.org/wiki/Dual_code) of `RM(r, m)` is `RM(m − r − 1, m)`, indicating the parity check matrix of `RM(r, m)` is the generator matrix for `RM(m - r - 1, m)`.
+
+See also: `ReedMuller`
+"""
+struct RecursiveReedMuller <: ClassicalCode
+    r::Int
+    m::Int
+
+    function RecursiveReedMuller(r, m)
+        if r < 0 || r > m
+            throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to valid code."))
+        end
+        new(r, m)
+    end
+end
+
+function _recursiveReedMuller(r::Int, m::Int)
+    if r == 1 && m == 1
+        return Matrix{Int}([1 1; 0 1])
+    elseif r == m
+        return Matrix{Int}(I, 2^m, 2^m)
+    elseif r == 0
+        return Matrix{Int}(ones(1, 2^m))
+    else
+        Gᵣₘ₋₁ = _recursiveReedMuller(r, m - 1)
+        Gᵣ₋₁ₘ₋₁ = _recursiveReedMuller(r - 1, m - 1)
+        return vcat(hcat(Gᵣₘ₋₁, Gᵣₘ₋₁), hcat(zeros(Int, size(Gᵣ₋₁ₘ₋₁)...), Gᵣ₋₁ₘ₋₁))
+    end
+end
+
+function generator(c::RecursiveReedMuller)
+    return _recursiveReedMuller(c.r, c.m)
+end
+
+function parity_checks(c::RecursiveReedMuller)
+    H = generator(RecursiveReedMuller(c.m - c.r - 1, c.m))
+    return H
+end
+
+code_n(c::RecursiveReedMuller) = 2 ^ c.m
+
+code_k(c::RecursiveReedMuller) = sum(binomial.(c.m, 0:c.r))
+
+distance(c::RecursiveReedMuller) = 2 ^ (c.m - c.r)
+
+rate(c::RecursiveReedMuller) = code_k(c::RecursiveReedMuller) / code_n(c::RecursiveReedMuller)

src/ecc/codes/classical/reedmuller.jl

@@ -1,38 +1,60 @@
 """The family of Reed-Muller codes, as discovered by Muller in his 1954 paper [muller1954application](@cite) and Reed who proposed the first efficient decoding algorithm [reed1954class](@cite).
 
+Let `m` be a positive integer and `r` a nonnegative integer with `r ≤ m`. These linear codes, denoted as `RM(r, m)`, have order `r` (where `0 ≤ r ≤ m`) and codeword length `n` of `2ᵐ`.
+
+Two special cases of generator(RM(r, m)) exist:
+    1. `generator(RM(0, m))`: This is the `0ᵗʰ`-order `RM` code, similar to the binary repetition code with length `2ᵐ`. It's characterized by a single basis vector containing all ones.
+    2. `generator(RM(m, m))`: This is the `mᵗʰ`-order `RM` code. It encompasses the entire field `F(2ᵐ)`, representing all possible binary strings of length `2ᵐ`.
+
 You might be interested in consulting [raaphorst2003reed](@cite), [abbe2020reed](@cite), and [djordjevic2021quantum](@cite) as well.
 
-The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller)
-"""
+The dimension of `RM(m - r - 1, m)` equals the dimension of the dual of `RM(r, m)`. Thus, `RM(m - r - 1, m) = RM(r, m)^⊥` shows that the [dual code](https://en.wikipedia.org/wiki/Dual_code) of `RM(r, m)` is `RM(m − r − 1, m)`, indicating the parity check matrix of `RM(r, m)` is the generator matrix for `RM(m - r - 1, m)`.
 
-abstract type ClassicalCode end
+The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller).
 
+See also: `RecursiveReedMuller`
+"""
 struct ReedMuller <: ClassicalCode
     r::Int
     m::Int
 
     function ReedMuller(r, m)
-        if r < 0 || m < 1 || m >= 11
-            throw(ArgumentError("Invalid parameters: r must be non-negative and m must be positive and < 11 in order to obtain a valid code and to remain tractable"))
+        if r < 0 || r > m
+            throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to valid code."))
         end
         new(r, m)
     end
 end
 
-function variables_xi(m, i)
-    return repeat([fill(1, 2^(m - i - 1)); fill(0, 2^(m - i - 1))], outer = 2^i)
+function _variablesₓᵢ_rm(m, i)
+    return repeat([fill(1, 2 ^ (m - i - 1)); fill(0, 2 ^ (m - i - 1))], outer = 2 ^ i)
 end
 
-function vmult(vecs...)
+function _vmult_rm(vecs...)
     return [reduce(*, a, init=1) for a in zip(vecs...)]
 end
 
-function parity_checks(c::ReedMuller)
-    r=c.r
-    m=c.m
-    xi = [variables_xi(m, i) for i in 0:m - 1]
-    row_matrices = [reduce(vmult, [xi[i + 1] for i in S], init = ones(Int, 2^m)) for s in 0:r for S in combinations(0:m - 1, s)]
+function generator(c::ReedMuller)
+    r = c.r
+    m = c.m
+    xᵢ = [_variablesₓᵢ_rm(m, i) for i in 0:m - 1]
+    row_matrices = [reduce(_vmult_rm, [xᵢ[i + 1] for i in S], init = ones(Int, 2 ^ m)) for s in 0:r for S in combinations(0:m - 1, s)]
     rows = length(row_matrices)
     cols = length(row_matrices[1])
-    H = reshape(vcat(row_matrices...), cols, rows)'
-end
+    G = reshape(vcat(row_matrices...), cols, rows)'
+    G = Matrix{Bool}(G)
+    return G 
+end
+
+function parity_checks(c::ReedMuller)
+    H = generator(ReedMuller(c.m - c.r - 1, c.m))
+    return H
+end
+
+code_n(c::ReedMuller) = 2 ^ c.m
+
+code_k(c::ReedMuller) = sum(binomial.(c.m, 0:c.r))
+
+distance(c::ReedMuller) = 2 ^ (c.m - c.r)
+
+rate(c::ReedMuller) = code_k(c::ReedMuller) / code_n(c::ReedMuller)

test/test_ecc_reedmuller.jl

@@ -1,68 +1,130 @@
 @testitem "Reed-Muller" begin
-    using Nemo
-    using Combinatorics
+    using Test
+    using Nemo: echelon_form, matrix, GF
     using LinearAlgebra
+    using QuantumClifford
     using QuantumClifford.ECC
-    using QuantumClifford.ECC: AbstractECC, ReedMuller
+    using QuantumClifford.ECC: AbstractECC, ReedMuller, generator, RecursiveReedMuller
 
-    function binomial_coeff_sum(r, m)
-        total = 0
-        for i in 0:r
-            total += length(combinations(1:m, i))
+    function designed_distance(matrix, m, r)
+        distance = 2 ^ (m - r)
+        for row in eachrow(matrix)
+            count = sum(row)
+            if count < distance 
+                return false
+            end
         end
-        return total
+        return true
     end
 
-    @testset "Test RM(r, m) Matrix Rank" begin
-        for m in 2:5
-            for r in 0:m - 1
-                H = parity_checks(ReedMuller(r, m))
-                mat = Nemo.matrix(Nemo.GF(2), H)
+    # Generate binary matrix representing all possible binary strings of length `2ᵐ` encompassing the entire field `F(2ᵐ)`.
+    function check_RM_m_m(m::Int)
+        n = 2^m
+        matrix = Matrix{Bool}(undef, n, n)
+        for i in 0:(n-1)
+            for j in 0:(n-1)
+                matrix[i+1, j+1] = ((i & j) == j)
+            end
+        end
+        return matrix
+    end
+
+    @testset "Test RM(r, m) properties" begin
+        for m in 3:10
+            for r in 3:m-1
+                H = generator(ReedMuller(r, m))
+                mat = matrix(GF(2), H)
                 computed_rank = LinearAlgebra.rank(mat)
-                expected_rank = binomial_coeff_sum(r, m)
+                expected_rank = sum(binomial.(m, 0:r))
                 @test computed_rank == expected_rank
+                @test designed_distance(H, m, r) == true
+                @test rate(ReedMuller(r, m)) == sum(binomial.(m, 0:r)) / 2 ^ m == rate(RecursiveReedMuller(r, m))
+                @test code_n(ReedMuller(r, m)) == 2 ^ m == code_n(RecursiveReedMuller(r, m))
+                @test code_k(ReedMuller(r, m)) == sum(binomial.(m, 0:r)) == code_k(RecursiveReedMuller(r, m))
+                @test distance(ReedMuller(r, m)) == 2 ^ (m - r) == distance(RecursiveReedMuller(r, m))
+                H₁ = generator(RecursiveReedMuller(r, m))
+                # generator(RecursiveReedMuller(r, m)) is canonically equivalent to the generator(ReedMuller(r, m)) under reduced row echelon form.
+                @test echelon_form(matrix(GF(2), Matrix{Int64}(H))) == echelon_form(matrix(GF(2), Matrix{Int64}(H₁)))
+                H = parity_checks(ReedMuller(r, m))
+                H₁ = parity_checks(RecursiveReedMuller(r, m))
+                # parity_checks(RecursiveReedMuller(r, m)) is canonically equivalent to the parity_checks(ReedMuller(r, m)) under reduced row echelon form.
+                @test echelon_form(matrix(GF(2), Matrix{Int64}(H))) == echelon_form(matrix(GF(2), Matrix{Int64}(H₁)))
+                # dim(ReedMuller(m - r - 1, m)) = dim(ReedMuller(r, m)^⊥). 
+                # ReedMuller(m - r - 1, m) = ReedMuller(r, m)^⊥ ∴ parity check matrix (H) of ReedMuller(r, m) is the generator matrix (G) for ReedMuller(m - r - 1, m).
+                H₁ = parity_checks(ReedMuller(m - r - 1, m))
+                G₂ = generator(ReedMuller(r, m))
+                @test size(H₁) == size(G₂)
+                @test echelon_form(matrix(GF(2), Matrix{Int64}(H₁))) == echelon_form(matrix(GF(2), Matrix{Int64}(G₂)))
+                G₃ = generator(ReedMuller(m - r - 1, m))
+                H₄ = parity_checks(ReedMuller(r, m))
+                @test size(G₃) == size(H₄)
+                # dim(RecursiveReedMuller(m - r - 1, m)) = dim(RecursiveReedMuller(r, m)^⊥). 
+                # RecursiveReedMuller(m - r - 1, m) = RecursiveReedMuller(r, m)^⊥ ∴ parity check matrix (H) of RecursiveReedMuller(r, m) is the generator matrix (G) for RecursiveReedMuller(m - r - 1, m).
+                H₁ = parity_checks(RecursiveReedMuller(m - r - 1, m))
+                G₂ = generator(RecursiveReedMuller(r, m))
+                @test echelon_form(matrix(GF(2), Matrix{Int64}(H₄))) == echelon_form(matrix(GF(2), Matrix{Int64}(G₃)))
+                @test echelon_form(matrix(GF(2), Matrix{Int64}(H₁))) == echelon_form(matrix(GF(2), Matrix{Int64}(G₂)))
+                G₃ = generator(RecursiveReedMuller(m - r - 1, m))
+                H₄ = parity_checks(RecursiveReedMuller(r, m))
+                @test size(G₃) == size(H₄)
+                @test echelon_form(matrix(GF(2), Matrix{Int64}(H₄))) == echelon_form(matrix(GF(2), Matrix{Int64}(G₃)))
             end
         end
     end
 
-    @testset "Testing common examples of RM(r,m) codes [raaphorst2003reed](@cite), [djordjevic2021quantum](@cite), [abbe2020reed](@cite)" begin
+    @testset "Test special case 1: generator(RM(0, m))" begin
+        for m in 3:10
+            H = generator(ReedMuller(0, m))
+            expected = ones(Int64, 1, 2^m)
+            @test H == expected
+        end
+    end
 
-        #RM(0,3)  
-        @test parity_checks(ReedMuller(0,3)) == [1 1 1 1 1 1 1 1]
+    @testset "Test special case 2: generator(RM(m, m))" begin
+        for m in 3:10
+            H = generator(ReedMuller(m, m))
+            expected = check_RM_m_m(m)
+            @test echelon_form(matrix(GF(2), Matrix{Int64}(H))) == echelon_form(matrix(GF(2), expected))
+        end
+    end
 
-        #RM(1,3) 
-        @test parity_checks(ReedMuller(1,3)) == [1 1 1 1 1 1 1 1;
-                                                 1 1 1 1 0 0 0 0;
-                                                 1 1 0 0 1 1 0 0;
-                                                 1 0 1 0 1 0 1 0]
-        #RM(2,3)
-        @test parity_checks(ReedMuller(2,3)) == [1 1 1 1 1 1 1 1;
-                                                 1 1 1 1 0 0 0 0;
-                                                 1 1 0 0 1 1 0 0;
-                                                 1 0 1 0 1 0 1 0;
-                                                 1 1 0 0 0 0 0 0;
-                                                 1 0 1 0 0 0 0 0;
-                                                 1 0 0 0 1 0 0 0]
-        #RM(3,3)
-        @test parity_checks(ReedMuller(3,3)) == [1 1 1 1 1 1 1 1;
-                                                 1 1 1 1 0 0 0 0;
-                                                 1 1 0 0 1 1 0 0;
-                                                 1 0 1 0 1 0 1 0;
-                                                 1 1 0 0 0 0 0 0;
-                                                 1 0 1 0 0 0 0 0;
-                                                 1 0 0 0 1 0 0 0;
-                                                 1 0 0 0 0 0 0 0]
-        #RM(2,4)
-        @test parity_checks(ReedMuller(2,4)) == [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
-                                                 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0;
-                                                 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0;
-                                                 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0;
-                                                 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0;
-                                                 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
-                                                 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
-                                                 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0;
-                                                 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0;
-                                                 1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0;
-                                                 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0]
+    @testset "Test common examples of RM(r,m) codes" begin
+        # Examples from [raaphorst2003reed](@cite), [djordjevic2021quantum](@cite), [abbe2020reed](@cite).
+        # RM(0,3)
+        @test generator(ReedMuller(0,3)) == [1 1 1 1 1 1 1 1]
+        # RM(1,3)
+        @test generator(ReedMuller(1,3)) == [1 1 1 1 1 1 1 1;
+                                             1 1 1 1 0 0 0 0;
+                                             1 1 0 0 1 1 0 0;
+                                             1 0 1 0 1 0 1 0]
+        # RM(2,3)
+        @test generator(ReedMuller(2,3)) == [1 1 1 1 1 1 1 1;
+                                             1 1 1 1 0 0 0 0;
+                                             1 1 0 0 1 1 0 0;
+                                             1 0 1 0 1 0 1 0;
+                                             1 1 0 0 0 0 0 0;
+                                             1 0 1 0 0 0 0 0;
+                                             1 0 0 0 1 0 0 0]
+        # RM(3,3)
+        @test generator(ReedMuller(3,3)) == [1 1 1 1 1 1 1 1;
+                                             1 1 1 1 0 0 0 0;
+                                             1 1 0 0 1 1 0 0;
+                                             1 0 1 0 1 0 1 0;
+                                             1 1 0 0 0 0 0 0;
+                                             1 0 1 0 0 0 0 0;
+                                             1 0 0 0 1 0 0 0;
+                                             1 0 0 0 0 0 0 0]
+        # RM(2,4)
+        @test generator(ReedMuller(2,4)) == [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
+                                             1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0;
+                                             1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0;
+                                             1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0;
+                                             1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0;
+                                             1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
+                                             1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
+                                             1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0;
+                                             1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0;
+                                             1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0;
+                                             1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0]
     end
 end

test/test_ecc_throws.jl

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