diff --git a/.github/workflows/benchmark-comment.yml b/.github/workflows/benchmark-comment.yml index fc39cc293..3ef7e6d10 100644 --- a/.github/workflows/benchmark-comment.yml +++ b/.github/workflows/benchmark-comment.yml @@ -26,7 +26,7 @@ jobs: - uses: actions/checkout@v4 # restore records from the artifacts - - uses: dawidd6/action-download-artifact@v6 + - uses: dawidd6/action-download-artifact@v7 with: workflow: benchmark.yml name: performance-tracking diff --git a/.github/workflows/ci-julia-nightly.yml b/.github/workflows/ci-julia-nightly.yml index d018ed30c..17506a342 100644 --- a/.github/workflows/ci-julia-nightly.yml +++ b/.github/workflows/ci-julia-nightly.yml @@ -44,7 +44,7 @@ jobs: JULIA_NUM_THREADS: ${{ matrix.threads }} JET_TEST: ${{ matrix.jet }} - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v5 with: file: lcov.info token: ${{ secrets.CODECOV_TOKEN }} \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index cf751e1c2..705c36598 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -56,7 +56,7 @@ jobs: env: JULIA_NUM_THREADS: ${{ matrix.threads }} - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v5 with: file: lcov.info token: ${{ secrets.CODECOV_TOKEN }} diff --git a/CHANGELOG.md b/CHANGELOG.md index 7bf3a570d..8b17c8dec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,15 @@ # News +## v0.9.16 - 2024-12-29 + +- 100× faster unbiased `random_pauli`. +- Enhancements to `GF(2)` Linear Algebra: unexported, experimental `gf2_row_echelon_with_pivots!`, `gf2_nullspace`, `gf2_rowspace_basis`. + +## v0.9.15 - 2024-12-22 + +- `pftrajectories` now supports fast multiqubit measurements with `PauliMeasurement` in addition to the already supported single qubit measurements `sMX/Z/Y` and workarounds like `naive_syndrome_circuit`. + ## v0.9.14 - 2024-11-03 - **(fix)** `affectedqubits()` on `sMX`, `sMY`, and `sMR*` @@ -75,7 +84,7 @@ - Gate errors are now conveniently supported by the various ECC benchmark setups in the `ECC` module. - Significant improvements to the low-level circuit compiler (the sumtype compactifier), leading to faster Pauli frame simulation of noisy circuits. - Bump `QuantumOpticsBase.jl` package extension compat bound. -- **(fix)** Remove printing of spurious debug info from the PyBP decoder. +- **(fix)** Remove printing of spurious debug info from the PyBP decoder. - **(fix)** Failed compactification of gates now only raises a warning instead of throwing an error. Defaults to slower non-compactified gates. ## v0.9.3 - 2024-04-10 @@ -93,7 +102,7 @@ - Implemented `iscss` function to identify whether a given code is known to be a CSS (Calderbank-Shor-Steane) code. - Added the classical Reed-Muller code in the ECC module. - Added the surface code to the ECC module. - + ## v0.9.0 - 2024-03-19 - **(breaking)** The defaults in `random_pauli` are now `realphase=true` and `nophase=true`. diff --git a/Project.toml b/Project.toml index e49f22503..342f02c3c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumClifford" uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1" authors = ["Stefan Krastanov and QuantumSavory community members"] -version = "0.9.14" +version = "0.9.16" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -48,15 +48,15 @@ Combinatorics = "1.0" DataStructures = "0.18" DocStringExtensions = "0.9" Graphs = "1.9" -Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34" +Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35" HostCPUFeatures = "0.1.6" ILog2 = "0.2.3, 1, 2" InteractiveUtils = "1.9" LDPCDecoders = "0.3.1" LinearAlgebra = "1.9" MacroTools = "0.5.9" -Makie = "0.20, 0.21" -Nemo = "0.42.1, 0.43, 0.44, 0.45, 0.46, 0.47" +Makie = "0.20, 0.21, 0.22" +Nemo = "0.42.1, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48" Plots = "1.38.0" PrecompileTools = "1.2" PyQDecoders = "0.2.1" diff --git a/README.md b/README.md index f64fd7ce5..0b2484b50 100644 --- a/README.md +++ b/README.md @@ -41,21 +41,21 @@ To install it use: ``` Works efficiently with -[pure](https://quantumsavory.github.io/QuantumClifford.jl/dev/manual/#Stabilizers-1) and -[mixed stabilizer](https://quantumsavory.github.io/QuantumClifford.jl/dev/mixed/#Mixed-Stabilizer-States-1) +[pure](https://qc.quantumsavory.org/dev/stab-algebra-manual/#Stabilizers) and +[mixed stabilizer](https://qc.quantumsavory.org/dev/mixed/#Mixed-Stabilizer-States) states of thousands of qubits as well as -[sparse or dense Clifford operations](https://quantumsavory.github.io/QuantumClifford.jl/dev/manual/#Clifford-Operators-1) +[sparse or dense Clifford operations](https://qc.quantumsavory.org/dev/stab-algebra-manual/#Clifford-Operators) acting upon them. -Implements [Pauli frames](https://quantumsavory.github.io/QuantumClifford.jl/dev/ecc_example_sim/) for fast sampling. +Implements [Pauli frames](https://qc.quantumsavory.org/dev/ecc_example_sim/) for fast sampling. Provides -[canonicalization](https://quantumsavory.github.io/QuantumClifford.jl/dev/manual/#Canonicalization-of-Stabilizers-1), -[projection](https://quantumsavory.github.io/QuantumClifford.jl/dev/manual/#Projective-Measurements-1), and -[generation](https://quantumsavory.github.io/QuantumClifford.jl/dev/manual/#Generating-a-Pauli-Operator-with-Stabilizer-Generators-1) operations, +[canonicalization](https://qc.quantumsavory.org/dev/stab-algebra-manual/#Canonicalization-of-Stabilizers), +[projection](http://qc.quantumsavory.org/dev/stab-algebra-manual/#Projective-Measurements), and +[generation](https://qc.quantumsavory.org/dev/stab-algebra-manual/#Generating-a-Pauli-Operator-with-Stabilizer-Generators) operations, as well as -[partial traces](https://quantumsavory.github.io/QuantumClifford.jl/dev/manual/#Partial-Traces-1). +[partial traces](https://qc.quantumsavory.org/dev/stab-algebra-manual/#Partial-Traces). ```jldoctest julia> P"X" * P"Z" diff --git a/docs/src/references.bib b/docs/src/references.bib index d10cf011a..348895335 100644 --- a/docs/src/references.bib +++ b/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}, @@ -562,6 +562,29 @@ @article{haah2011local 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} +} + @article{hamming1950error, title={Error detecting and error correcting codes}, author={Hamming, Richard W}, @@ -572,10 +595,3 @@ @article{hamming1950error year={1950}, publisher={Nokia Bell Labs} } - -@book{huffman2010fundamentals, - title={Fundamentals of error-correcting codes}, - author={Huffman, W Cary and Pless, Vera}, - year={2010}, - publisher={Cambridge university press} -} diff --git a/docs/src/references.md b/docs/src/references.md index 9c30ff88a..c389e57dc 100644 --- a/docs/src/references.md +++ b/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) - [hamming1950error](@cite) - [huffman2010fundamentals](@cite) diff --git a/ext/QuantumCliffordHeckeExt/lifted_product.jl b/ext/QuantumCliffordHeckeExt/lifted_product.jl index 66bd35ed6..e5c957229 100644 --- a/ext/QuantumCliffordHeckeExt/lifted_product.jl +++ b/ext/QuantumCliffordHeckeExt/lifted_product.jl @@ -28,7 +28,7 @@ During the construction, we do arithmetic operations to get the group algebra el Here `x` is the generator of the group algebra, i.e., offset-1 cyclic permutation, and `GA(1)` is the unit element. ```jldoctest -julia> import Hecke: group_algebra, GF, abelian_group, gens; import LinearAlgebra: diagind; +julia> import Hecke: group_algebra, GF, abelian_group, gens; import LinearAlgebra: diagind; using QuantumClifford.ECC; julia> l = 63; GA = group_algebra(GF(2), abelian_group(l)); x = gens(GA)[]; @@ -52,10 +52,12 @@ julia> code_n(c1), code_k(c1) (882, 24) ``` -A [[175, 19, d ≤ 0]] code from Eq. (18) in Appendix A of [raveendran2022finite](@cite), +A [[175, 19, d ≤ 10]] code from Eq. (18) in Appendix A of [raveendran2022finite](@cite), following the 4th constructor. ```jldoctest +julia> import Hecke; using QuantumClifford.ECC; + julia> base_matrix = [0 0 0 0; 0 1 2 5; 0 6 3 1]; l = 7; julia> c2 = LPCode(base_matrix, l .- base_matrix', l); @@ -162,15 +164,15 @@ Here is an example of a [[56, 28, 2]] 2BGA code from Table 2 of [lin2024quantum] with direct product of `C₄ x C₂`. ```jldoctest -julia> import Hecke: group_algebra, GF, abelian_group, gens +julia> import Hecke: group_algebra, GF, abelian_group, gens; using QuantumClifford.ECC; julia> GA = group_algebra(GF(2), abelian_group([14,2])); julia> x, s = gens(GA); -julia> A = 1 + x^7 +julia> A = 1 + x^7; -julia> B = 1 + x^7 + s + x^8 + s*x^7 + x +julia> B = 1 + x^7 + s + x^8 + s*x^7 + x; julia> c = two_block_group_algebra_codes(A,B); @@ -189,7 +191,7 @@ The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/qcga A [[756, 16, ≤ 34]] code from Table 3 of [bravyi2024high](@cite): ```jldoctest -julia> import Hecke: group_algebra, GF, abelian_group, gens +julia> import Hecke: group_algebra, GF, abelian_group, gens; using QuantumClifford.ECC; julia> l=21; m=18; @@ -215,10 +217,14 @@ where `𝐺ᵣ = ℤ/l₁ × ℤ/l₂ × ... × ℤ/lᵣ`. A [[48, 4, 6]] Weight-6 TB-QLDPC code from Appendix A Table 2 of [voss2024multivariatebicyclecodes](@cite). ```jldoctest -julia> import Hecke: group_algebra, GF, abelian_group, gens; # hide +julia> import Hecke: group_algebra, GF, abelian_group, gens; using QuantumClifford.ECC; julia> l=4; m=6; +julia> GA = group_algebra(GF(2), abelian_group([l, m])); + +julia> x, y = gens(GA); + julia> z = x*y; julia> A = x^3 + y^5; @@ -238,10 +244,10 @@ where `𝑙` and `𝑚` are coprime, and can be expressed as univariate polynomi with generator `𝜋 = 𝑥𝑦`. They can be viewed as a special case of Lifted Product construction based on abelian group `ℤₗ x ℤₘ` where `ℤⱼ` cyclic group of order `j`. -[108, 12, 6]] coprime-bivariate bicycle (BB) code from Table 2 of [wang2024coprime](@cite). +[[108, 12, 6]] coprime-bivariate bicycle (BB) code from Table 2 of [wang2024coprime](@cite). ```jldoctest -julia> import Hecke: group_algebra, GF, abelian_group, gens; +julia> import Hecke: group_algebra, GF, abelian_group, gens; using QuantumClifford.ECC; julia> l=2; m=27; @@ -252,6 +258,10 @@ julia> 𝜋 = gens(GA)[1]; julia> A = 𝜋^2 + 𝜋^5 + 𝜋^44; julia> B = 𝜋^8 + 𝜋^14 + 𝜋^47; + +julia> c = two_block_group_algebra_codes(A, B); + +julia> code_n(c), code_k(c) (108, 12) ``` @@ -271,6 +281,8 @@ See also: [`two_block_group_algebra_codes`](@ref), [`bicycle_codes`](@ref). A [[254, 28, 14 ≤ d ≤ 20]] code from (A1) in Appendix B of [panteleev2021degenerate](@cite). ```jldoctest +julia> import Hecke; using QuantumClifford.ECC + julia> c = generalized_bicycle_codes([0, 15, 20, 28, 66], [0, 58, 59, 100, 121], 127); julia> code_n(c), code_k(c) @@ -281,6 +293,8 @@ An [[70, 8, 10]] *abelian* 2BGA code from Table 1 of [lin2024quantum](@cite), wi order `l = 35`, illustrates that *abelian* 2BGA codes can be viewed as GB codes. ```jldoctest +julia> import Hecke; using QuantumClifford.ECC + julia> l = 35; julia> c1 = generalized_bicycle_codes([0, 15, 16, 18], [0, 1, 24, 27], l); @@ -317,6 +331,8 @@ code with the group `G = ℤ₃ˣ³` corresponds to a cubic code. The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/haah_cubic). ```jldoctest +julia> import Hecke; using QuantumClifford.ECC; + julia> c = haah_cubic_codes([0, 15, 20, 28, 66], [0, 58, 59, 100, 121], 6); julia> code_n(c), code_k(c) diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl index 6b90af7a8..25f84a218 100644 --- a/src/QuantumClifford.jl +++ b/src/QuantumClifford.jl @@ -1168,6 +1168,69 @@ function gf2_H_to_G(H) G[:,invperm(sindx)] end +"""Performs in-place Gaussian elimination on a binary matrix and returns +its *row echelon form*,*rank*, the *transformation matrix*, and the *pivot +columns*. The transformation matrix that converts the original matrix into +the row echelon form. The `full` parameter controls the extent of elimination: +if `true`, only rows below the pivot are affected; if `false`, both above and +below the pivot are eliminated.""" +function gf2_row_echelon_with_pivots!(M::AbstractMatrix{Int}; full=false) + r, c = size(M) + N = Matrix{Int}(LinearAlgebra.I, r, r) + p = 1 + pivots = Int[] + for col in 1:c + @inbounds for row in p:r + if M[row, col] == 1 + if row != p + M[[row, p], :] .= M[[p, row], :] + N[[row, p], :] .= N[[p, row], :] + end + break + end + end + if M[p, col] == 1 + if !full + elim_range = p+1:r + else + elim_range = 1:r + end + @simd for j in elim_range + @inbounds if j != p && M[j, col] == 1 + M[j, :] .= (M[j, :] .+ M[p, :]) .% 2 + N[j, :] .= (N[j, :] .+ N[p, :]) .% 2 + end + end + p += 1 + push!(pivots, col) + end + if p > r + break + end + end + rank = p - 1 + return M, rank, N, pivots +end + +"""The nullspace of a binary matrix.""" +function gf2_nullspace(H::AbstractMatrix{Int}) + m = size(H',1) + _, matrix_rank, transformation_matrix, _ = gf2_row_echelon_with_pivots!(copy(H)') + if m == matrix_rank + # By the rank-nullity theorem, if rank(M) = m, then nullity(M) = 0 + return zeros(Bool, 1, m) + end + # Extract the nullspace from the transformation matrix + return transformation_matrix[matrix_rank+1:end, :] +end + +"""The basis for the row space of the binary matrix.""" +function gf2_rowspace_basis(H::AbstractMatrix{Int}) + pivots = gf2_row_echelon_with_pivots!(copy(H)')[4] + # Extract the rows corresponding to the pivot columns + H[pivots,:] +end + ############################## # Error classes ############################## diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl index bcc7e43b9..d249c80f0 100644 --- a/src/ecc/ECC.jl +++ b/src/ecc/ECC.jl @@ -388,6 +388,7 @@ include("codes/quantumreedmuller.jl") include("codes/classical/reedmuller.jl") include("codes/classical/recursivereedmuller.jl") include("codes/classical/bch.jl") +include("codes/classical/golay.jl") include("codes/classical/hamming.jl") # qLDPC diff --git a/src/ecc/codes/classical/golay.jl b/src/ecc/codes/classical/golay.jl new file mode 100644 index 000000000..4770a5d45 --- /dev/null +++ b/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 diff --git a/src/pauli_frames.jl b/src/pauli_frames.jl index af84bd63b..53a60d940 100644 --- a/src/pauli_frames.jl +++ b/src/pauli_frames.jl @@ -11,12 +11,14 @@ struct PauliFrame{T,S} <: AbstractQCState measurements::S # TODO check if when looping over this we are actually looping over the fast axis end -nqubits(f::PauliFrame) = nqubits(f.frame) +nqubits(f::PauliFrame) = nqubits(f.frame) - 1 # dont count the ancilla qubit Base.length(f::PauliFrame) = size(f.measurements, 1) Base.eachindex(f::PauliFrame) = 1:length(f) Base.copy(f::PauliFrame) = PauliFrame(copy(f.frame), copy(f.measurements)) Base.view(frame::PauliFrame, r) = PauliFrame(view(frame.frame, r), view(frame.measurements, r, :)) +tab(f::PauliFrame) = Tableau(f.frame.tab.phases, nqubits(f), f.frame.tab.xzs) + fastrow(s::PauliFrame) = PauliFrame(fastrow(s.frame), s.measurements) fastcolumn(s::PauliFrame) = PauliFrame(fastcolumn(s.frame), s.measurements) @@ -26,7 +28,8 @@ $(TYPEDSIGNATURES) Prepare an empty set of Pauli frames with the given number of `frames` and `qubits`. Preallocates spaces for `measurement` number of measurements. """ function PauliFrame(frames, qubits, measurements) - stab = fastcolumn(zero(Stabilizer, frames, qubits)) # TODO this should really be a Tableau + # one extra qubit for ancilla measurements + stab = fastcolumn(zero(Stabilizer, frames, qubits + 1)) # TODO this should really be a Tableau bits = zeros(Bool, frames, measurements) frame = PauliFrame(stab, bits) initZ!(frame) @@ -112,6 +115,24 @@ function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRY, and faster sMRX return frame end +function apply!(frame::PauliFrame, op::PauliMeasurement) + # this is inspired by ECC.naive_syndrome_circuit + n = nqubits(op.pauli) + for qubit in 1:n + if op.pauli[qubit] == (1, 0) + apply!(frame, sXCZ(qubit, n + 1)) + elseif op.pauli[qubit] == (0, 1) + apply!(frame, sCNOT(qubit, n + 1)) + elseif op.pauli[qubit] == (1, 1) + apply!(frame, sYCX(qubit, n + 1)) + end + end + op.pauli.phase[] == 0 || apply!(frame, sX(n + 1)) + apply!(frame, sMRZ(n + 1, op.bit)) + + return frame +end + function applynoise!(frame::PauliFrame,noise::UnbiasedUncorrelatedNoise,i::Int) p = noise.p xzs = tab(frame.frame).xzs diff --git a/src/randoms.jl b/src/randoms.jl index c87932fd4..efac990e4 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -1,4 +1,4 @@ -using Random: randperm, AbstractRNG, GLOBAL_RNG +using Random: randperm, AbstractRNG, GLOBAL_RNG, rand! using ILog2 import Nemo @@ -21,11 +21,24 @@ function random_pauli end """An in-place version of [`random_pauli`](@ref)""" function random_pauli! end + +# Modified from `Base` for `BitArray` +@inline _msk_end(::Type{T}, l::Int) where {T<:Unsigned} = ~T(0) >>> _mod(T, -l) + +# A mask for the non-coding bits in the last chunk. +@inline _msk_end(P::PauliOperator) = _msk_end(eltype(P.xz), length(P)) + +# Unset the leftover bits in the data array that don't code the Pauli operator. +@inline function _unset_noncoding_bits!(P::PauliOperator) + msk = _msk_end(P) + P.xz[end] &= msk + P.xz[end÷2] &= msk + nothing +end + function random_pauli!(rng::AbstractRNG, P::PauliOperator; nophase=true, realphase=true) - n = nqubits(P) - for i in 1:n - P[i] = rand(rng, (true, false)), rand(rng, (true,false)) - end + rand!(rng, P.xz) + _unset_noncoding_bits!(P) P.phase[] = nophase ? 0x0 : (realphase ? rand(rng,(0x0,0x2)) : rand(rng,0x0:0x3)) P end diff --git a/src/sumtypes.jl b/src/sumtypes.jl index 9a2cc4462..f496532b0 100644 --- a/src/sumtypes.jl +++ b/src/sumtypes.jl @@ -223,6 +223,12 @@ function concrete_typeparams(t::Type{NoiseOp}) ] end +function concrete_typeparams(::Type{PauliMeasurement}) + return [ + (Array{UInt8,0}, Vector{UInt64}), + ] +end + # XXX This has to happen after defining all the `concrete_typeparams` methods diff --git a/test/Project.toml b/test/Project.toml index eb254c893..4ba671eb1 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/test_bitpack.jl b/test/test_bitpack.jl index e5aae358c..be7ad212f 100644 --- a/test/test_bitpack.jl +++ b/test/test_bitpack.jl @@ -1,6 +1,6 @@ @testitem "Alternative bit packing" tags=[:bitpack] begin using Random - using QuantumClifford: Tableau + using QuantumClifford: Tableau, mul_left! @testset "alternative bit packing" begin for n in [1,3] # can not go higher than 4 (limitation from SIMD acting on transposed/strided arrays) @@ -20,6 +20,8 @@ after_cliff = stab_to_gf2(_after_clif); after_cliff_phases = phases(_after_clif); + after_mul = stab_to_gf2(mul_left!(copy(s64), p64)) + for int in [UInt8, UInt16, UInt32, UInt64] p = PauliOperator(p64.phase, N, collect(reinterpret(int,p64.xz))); xzs = collect(reinterpret(int, collect(xzs64))); @@ -45,8 +47,38 @@ @test after_cliff == stab_to_gf2(after_clifford) @test after_cliff_phases == phases(after_clifford) end + + @test after_mul == stab_to_gf2(mul_left!(copy(s), p)) end end end end + + @testset "fast column and fast row mul_left correctness" begin + reinterpret_stab(s) = Stabilizer(Tableau(copy(phases(s)), nqubits(s), collect(reinterpret(UInt8, collect(s.tab.xzs))))) + reinterpret_p(p) = PauliOperator(p.phase, nqubits(p), collect(reinterpret(UInt8, p.xz))) + for N in [7,8,9,33,61,62,63,64,65,66] + + s0 = random_stabilizer(2,N) + p = random_pauli(N) + s = copy(s0) + sr = QuantumClifford.fastrow(copy(s)) + sc = QuantumClifford.fastcolumn(copy(s)) + + mul_left!(s, p) + mul_left!(sr, p) + mul_left!(sc, p) + + s8 = reinterpret_stab(copy(s0)) + s8r = QuantumClifford.fastrow(copy(s8)) + s8c = QuantumClifford.fastcolumn(copy(s8)) + p8 = reinterpret_p(copy(p)) + mul_left!(s8, p8) + mul_left!(s8r, p8) + mul_left!(s8c, p8) + + @test stab_to_gf2(s) == stab_to_gf2(sr) == stab_to_gf2(sc) == stab_to_gf2(s8) == stab_to_gf2(s8r) == stab_to_gf2(s8c) + end + end + end end diff --git a/test/test_doctests.jl b/test/test_doctests.jl index 0ef1a0f23..b3b1ceb32 100644 --- a/test/test_doctests.jl +++ b/test/test_doctests.jl @@ -1,12 +1,19 @@ @testitem "Doctests" tags=[:doctests] begin using Documenter using QuantumClifford + using QuantumClifford.Experimental.NoisyCircuits + using QuantumInterface + + ENV["HECKE_PRINT_BANNER"] = "false" + import Hecke + const QuantumCliffordHeckeExt = Base.get_extension(QuantumClifford, :QuantumCliffordHeckeExt) ENV["LINES"] = 80 # for forcing `displaysize(io)` to be big enough ENV["COLUMNS"] = 80 DocMeta.setdocmeta!(QuantumClifford, :DocTestSetup, :(using QuantumClifford; using QuantumClifford.ECC); recursive=true) + modules = [QuantumClifford, QuantumClifford.Experimental.NoisyCircuits, QuantumClifford.ECC, QuantumInterface, QuantumCliffordHeckeExt] doctestfilters = [r"(QuantumClifford\.|)"] - doctest(QuantumClifford; + doctest(nothing, modules; doctestfilters #fix=true ) diff --git a/test/test_ecc_base.jl b/test/test_ecc_base.jl index 7c1d12ef5..39a1ab2e9 100644 --- a/test/test_ecc_base.jl +++ b/test/test_ecc_base.jl @@ -160,11 +160,14 @@ const code_instance_args = Dict( function all_testablable_code_instances(;maxn=nothing) codeinstances = [] + i = 1 for t in subtypes(QuantumClifford.ECC.AbstractECC) for c in get(code_instance_args, t.name.name, []) codeinstance = t(c...) !isnothing(maxn) && nqubits(codeinstance) > maxn && continue push!(codeinstances, codeinstance) + #@show i, t, code_n(codeinstance), code_k(codeinstance), code_s(codeinstance), code_n(codeinstance)-code_k(codeinstance) + i += 1 end end return codeinstances diff --git a/test/test_ecc_golay.jl b/test/test_ecc_golay.jl new file mode 100644 index 000000000..1f81b7b0e --- /dev/null +++ b/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 diff --git a/test/test_ecc_syndromes.jl b/test/test_ecc_syndromes.jl index 0f6f67ea0..01fbc5469 100644 --- a/test/test_ecc_syndromes.jl +++ b/test/test_ecc_syndromes.jl @@ -5,6 +5,11 @@ include("test_ecc_base.jl") + using QuantumClifford: Tableau + reinterpret_frame(frame) = PauliFrame(reinterpret_stab(frame.frame), copy(frame.measurements)) + reinterpret_stab(s) = Stabilizer(Tableau(copy(phases(s)), nqubits(s), collect(reinterpret(UInt8, collect(s.tab.xzs)))[[1:1+(nqubits(s)-1)÷8;end÷2+1:end÷2+1+(nqubits(s)-1)÷8],:])) + reinterpret_p(p) = PauliOperator(p.phase, nqubits(p), collect(reinterpret(UInt8, p.xz))[[1:1+(nqubits(p)-1)÷8;end÷2+1:end÷2+1+(nqubits(p)-1)÷8]]) + function pframe_naive_vs_shor_syndrome(code) ecirc = naive_encoding_circuit(code) naive_scirc, naive_ancillaries = naive_syndrome_circuit(code) @@ -30,19 +35,35 @@ pftrajectories(shor_frames, vcat(ecirc, shor_cat_scirc)) # manually injecting the same type of noise in the frames -- not really a user accessible API p = random_pauli(dataqubits, realphase=true) - pₙ = embed(naive_qubits, 1:dataqubits, p) - pₛ = embed(shor_qubits, 1:dataqubits, p) + pₙ = embed(naive_qubits+1, 1:dataqubits, p) # +1 to account for the buffer qubit hidden in pauli frames + pₛ = embed(shor_qubits+1, 1:dataqubits, p) # +1 to account for the buffer qubit hidden in pauli frames mul_left!(naive_frames.frame, pₙ) mul_left!(shor_frames.frame, pₛ) # run the syndrome circuits using the public API pftrajectories(naive_frames, naive_scirc) pftrajectories(shor_frames, shor_scirc) @test pfmeasurements(naive_frames) == pfmeasurements(shor_frames)[:,shor_bits] + + # just for completeness, let's also try bitpacking in UInt8 instead of the default UInt + _naive_frames = PauliFrame(nframes, naive_qubits, syndromebits) + _shor_frames = PauliFrame(nframes, shor_qubits, last(shor_bits)) + naive_uint8 = reinterpret_frame(_naive_frames) + shor_uint8 = reinterpret_frame(_shor_frames) + pftrajectories(naive_uint8, ecirc) + pftrajectories(shor_uint8, vcat(ecirc, shor_cat_scirc)) + p_uint8 = reinterpret_p(p) + pₙ_uint8 = embed(naive_qubits+1, 1:dataqubits, p_uint8) + pₛ_uint8 = embed(shor_qubits+1, 1:dataqubits, p_uint8) + mul_left!(naive_uint8.frame, pₙ_uint8) + mul_left!(shor_uint8.frame, pₛ_uint8) + pftrajectories(naive_uint8, naive_scirc) + pftrajectories(shor_uint8, shor_scirc) + @test pfmeasurements(shor_uint8)[:,shor_bits] == pfmeasurements(shor_frames)[:,shor_bits] == pfmeasurements(naive_frames) == pfmeasurements(naive_uint8) end end @testset "naive and shor measurement circuits" begin - for c in all_testablable_code_instances() + for (i,c) in enumerate(all_testablable_code_instances()) pframe_naive_vs_shor_syndrome(c) end end diff --git a/test/test_ecc_throws.jl b/test/test_ecc_throws.jl index b43a6f16e..f831b7ea2 100644 --- a/test/test_ecc_throws.jl +++ b/test/test_ecc_throws.jl @@ -1,6 +1,6 @@ @testitem "ECC throws" begin - using QuantumClifford.ECC: ReedMuller, BCH, RecursiveReedMuller, Hamming + using QuantumClifford.ECC: ReedMuller, BCH, RecursiveReedMuller, Golay, Hamming @test_throws ArgumentError ReedMuller(-1, 3) @test_throws ArgumentError ReedMuller(1, 0) @@ -13,5 +13,6 @@ @test_throws ArgumentError RecursiveReedMuller(1, 0) @test_throws ArgumentError RecursiveReedMuller(4, 2) + @test_throws ArgumentError Golay(21) @test_throws ArgumentError Hamming(1) end diff --git a/test/test_jet.jl b/test/test_jet.jl index 7f0248ab7..ef6374272 100644 --- a/test/test_jet.jl +++ b/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 diff --git a/test/test_pauliframe.jl b/test/test_pauliframe.jl index cc335934c..ece38ea30 100644 --- a/test/test_pauliframe.jl +++ b/test/test_pauliframe.jl @@ -88,4 +88,29 @@ @test all(0.25.*[1 0 0 0 1] .<= (sum(m, dims=1)[:,1:5])./n .<= 0.75.*[1 0 0 0 1]) end end + + @testset "PauliMeasurements" begin + n = 2000 + state = Register(one(MixedDestabilizer, 3), 5) + frame = PauliFrame(n, 3, 5) + + glassy_ghz_circuit = [ + sHadamard(1), sHadamard(2), sHadamard(3), + PauliMeasurement(P"ZZ_", 1), PauliMeasurement(P"_ZZ", 2), + sMZ(1, 3), sMZ(2, 4), sMZ(3, 5) + ] + for m in [pfmeasurements(pftrajectories(copy(frame), glassy_ghz_circuit)), + pfmeasurements(pftrajectories(glassy_ghz_circuit; trajectories=n, threads=false)), + pfmeasurements(pftrajectories(glassy_ghz_circuit; trajectories=n, threads=true))] + + # decode based on measurement outcomes + for r in eachrow(m) + r[4] ⊻= r[1] + r[5] ⊻= r[1] ⊻ r[2] + end + + # check that the correct correlations are present + @test all(m[:, 3] .== m[:, 4] .== m[:, 5]) + end + end end diff --git a/test/test_row_echelon_with_pivots.jl b/test/test_row_echelon_with_pivots.jl new file mode 100644 index 000000000..98d72429e --- /dev/null +++ b/test/test_row_echelon_with_pivots.jl @@ -0,0 +1,165 @@ +@testitem "Row echelon with pivots" begin + using Random + using Nemo + using Nemo: echelon_form, matrix, GF + using QuantumClifford + using QuantumClifford: gf2_row_echelon_with_pivots!, gf2_nullspace, gf2_rowspace_basis + test_sizes = [1,2,10,63,64,65,127,128,129] + + @testset "GF(2) row echelon form with transformation matrix, pivots etc." begin + for n in test_sizes + for rep in 1:10 + gf2_matrices = [rand(Bool, size, size) for size in test_sizes] + for (i, mat) in enumerate(gf2_matrices) + naive_echelon_form, _, transformation, _ = gf2_row_echelon_with_pivots!(Matrix{Int}(mat), full=true) # in-place + # Check the correctness of the transformation matrix + @test (transformation*mat) .%2 == naive_echelon_form + # Check the correctness of Gaussian elimination + @test naive_echelon_form == gf2_gausselim!(mat) + # Consistency check with Nemo.jl's echelon_form + nemo_mat = matrix(GF(2), Matrix{Int}(mat)) + @test echelon_form(nemo_mat) == matrix(GF(2), naive_echelon_form) + end + end + end + end + + function is_in_nullspace(A, x) + # Ensure x is the correct orientation + if size(x, 1) != size(A, 2) + x = transpose(x) + end + # Perform modulo 2 arithmetic: A * x must be zero mod 2 + if size(x, 2) == 1 # x is a single column vector + result = A * x + return all(result .% 2 .== 0) # Check if A * x = 0 mod 2 + else # x is a matrix, check each column vector + for i in 1:size(x, 2) + result = A * x[:, i] # Multiply A with the i-th column of x + if !all(result .% 2 .== 0) # Check if A * column = 0 mod 2 + return false + end + end + return true # All columns are in the null space mod 2 + end + end + + @testset "GF(2) nullspace of the binary matrix" begin + for n in test_sizes + for rep in 1:10 + gf2_matrices = [rand(Bool, size, size) for size in test_sizes] + for (i, matrix) in enumerate(gf2_matrices) + imat = Matrix{Int}(matrix) + ns = gf2_nullspace(imat) + @test is_in_nullspace(imat, ns) + end + end + end + end + + @testset "Consistency check with ldpc" begin + # sanity checks for comparison to https://github.com/quantumgizmos/ldpc + # results compared with 'from ldpc.mod2 import nullspace, row_basis, row_echelon' + # Consistency check 1 + H = [1 1 1; 1 1 1; 0 1 0] + echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place + @test echelon_form == [1 1 1; 0 1 0; 0 0 0] + @test rank == 2 + @test transformation == [1 0 0; 0 0 1; 1 1 0] + @test pivots == [1, 2] # in python, it's [0, 1] due to zero-indexing + @test mod.((transformation*copy(H)), 2) == echelon_form + @test gf2_nullspace(copy(H)) == [1 0 1] + @test gf2_rowspace_basis(copy(H)) == [1 1 1; 0 1 0] + # Consistency check 2 + H = [0 0 0 1 1 1 1; + 0 1 1 0 0 1 1; + 1 0 1 0 1 0 1] + echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place + @test echelon_form == [1 0 1 0 1 0 1; + 0 1 1 0 0 1 1; + 0 0 0 1 1 1 1] + @test rank == 3 + @test transformation == [0 0 1; + 0 1 0; + 1 0 0] + @test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing + @test mod.((transformation*copy(H)), 2) == echelon_form + @test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0; + 0 1 1 1 1 0 0; + 0 1 0 1 0 1 0; + 0 0 1 1 0 0 1] + @test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1; + 0 1 1 0 0 1 1; + 1 0 1 0 1 0 1] + # Consistency check 3 + H = [1 1 0; 0 1 1; 1 0 1] + echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place + @test echelon_form == [1 1 0; + 0 1 1; + 0 0 0] + @test rank == 2 + @test transformation == [1 0 0; + 0 1 0; + 1 1 1] + @test pivots == [1,2 ] # in python, it's [0, 1] due to zero-indexing + @test mod.((transformation*copy(H)), 2) == echelon_form + @test gf2_nullspace(copy(H)) == [1 1 1] + @test gf2_rowspace_basis(copy(H)) == [1 1 0; + 0 1 1] + # Consistency check 4 + H = [1 1 0; 0 1 0; 0 0 1] + echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place + @test echelon_form == [1 1 0; + 0 1 0; + 0 0 1] + @test rank == 3 + @test transformation == [1 0 0; + 0 1 0; + 0 0 1] + @test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing + @test mod.((transformation*copy(H)), 2) == echelon_form + @test gf2_nullspace(copy(H)) == [0 0 0] + @test gf2_rowspace_basis(copy(H)) == [1 1 0; + 0 1 0; + 0 0 1] + # Consistency check 5 + H = [1 1 0; 0 1 0; 0 0 1; 0 1 1] + echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place + @test echelon_form == [1 1 0; + 0 1 0; + 0 0 1; + 0 0 0] + @test rank == 3 + @test transformation == [1 0 0 0; + 0 1 0 0; + 0 0 1 0; + 0 1 1 1] + @test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing + @test mod.((transformation*copy(H)), 2) == echelon_form + @test gf2_nullspace(copy(H)) == [0 0 0] + @test gf2_rowspace_basis(copy(H)) == [1 1 0; + 0 1 0; + 0 0 1] + # Consistency check 6 + H = [0 0 0 1 1 1 1; + 0 1 1 0 0 1 1; + 1 0 1 0 1 0 1] + echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place + @test echelon_form == [1 0 1 0 1 0 1; + 0 1 1 0 0 1 1; + 0 0 0 1 1 1 1] + @test rank == 3 + @test transformation == [0 0 1; + 0 1 0; + 1 0 0] + @test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing + @test mod.((transformation*copy(H)), 2) == echelon_form + @test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0; + 0 1 1 1 1 0 0; + 0 1 0 1 0 1 0; + 0 0 1 1 0 0 1] + @test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1; + 0 1 1 0 0 1 1; + 1 0 1 0 1 0 1] + end +end