diff --git a/CHANGELOG.md b/CHANGELOG.md index 604e855c9..6d46dc8c4 100644 --- a/CHANGELOG.md +++ b/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). diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl index ed50dba5b..5498f2ad9 100644 --- a/src/ecc/ECC.jl +++ b/src/ecc/ECC.jl @@ -1,6 +1,7 @@ module ECC using LinearAlgebra +using LinearAlgebra: I using QuantumClifford using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits @@ -8,7 +9,7 @@ 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 diff --git a/src/ecc/codes/classical/recursivereedmuller.jl b/src/ecc/codes/classical/recursivereedmuller.jl new file mode 100644 index 000000000..ae30ae578 --- /dev/null +++ b/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) diff --git a/src/ecc/codes/classical/reedmuller.jl b/src/ecc/codes/classical/reedmuller.jl index 7f9c1f668..eea827661 100644 --- a/src/ecc/codes/classical/reedmuller.jl +++ b/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 \ No newline at end of file + 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) diff --git a/test/test_ecc_reedmuller.jl b/test/test_ecc_reedmuller.jl index 12031d71f..e93cc2648 100644 --- a/test/test_ecc_reedmuller.jl +++ b/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 diff --git a/test/test_ecc_throws.jl b/test/test_ecc_throws.jl index bd648327b..53055f8cf 100644 --- a/test/test_ecc_throws.jl +++ b/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