[2ʳ-1, 2ʳ-1-r, 3] Hamming code

docs/src/references.bib

@@ -584,3 +584,14 @@ @article{bhatia2018mceliece
   journal={arXiv preprint arXiv:1811.06246},
   year={2018}
 }
+
+@article{hamming1950error,
+  title={Error detecting and error correcting codes},
+  author={Hamming, Richard W},
+  journal={The Bell system technical journal},
+  volume={29},
+  number={2},
+  pages={147--160},
+  year={1950},
+  publisher={Nokia Bell Labs}
+}

docs/src/references.md

@@ -59,6 +59,8 @@ For classical code construction routines:
 - [golay1949notes](@cite)
 - [huffman2010fundamentals](@cite)
 - [bhatia2018mceliece](@cite)
+- [hamming1950error](@cite)
+- [huffman2010fundamentals](@cite)
 
 # References
 

src/ecc/ECC.jl

@@ -389,6 +389,7 @@ 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
 include("codes/classical/lifted.jl")

src/ecc/codes/classical/hamming.jl

@@ -0,0 +1,61 @@
+using SparseArrays
+
+abstract type ClassicalCode end
+
+"""
+The family of `[2ʳ - 1, 2ʳ - 1 - r, 3]` Hamming binary codes were discovered
+by Richard W. Hamming in his 1950 paper [hamming1950error](@cite) as a way of
+automatically correcting errors introduced by punched card readers. In his
+original paper, Hamming elaborated his general idea, but specifically focused
+on the `Hamming(7, 4)` code which adds three parity bits to four bits of data.
+
+The Hamming matrix `H` is an `r × (2ʳ - 1)` binary matrix, where each column
+corresponds to the binary representation of the integers from 1 to 2ʳ - 1, with
+`r ≥ 2`. This matrix serves as the parity-check matrix for a binary Hamming code
+with parameters `[2ʳ − 1, 2ʳ − 1 − r, 3]` The minimum Hamming distance of this
+code is 3, as detailed in [huffman2010fundamentals](@cite).
+
+The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/hamming).
+"""
+struct Hamming <: ClassicalCode
+    r::Int
+    function Hamming(r)
+        if r < 2
+            throw(ArgumentError("Invalid parameters: `r` must be ≥ 2 to obtain a valid code."))
+        end
+        new(r)
+    end
+end
+
+function parity_checks(h::Hamming)
+    n = 2^h.r - 1 # Number of columns in H
+    max_elem = n * h.r # Max non-zero entries in H
+    # Pre-allocate arrays for sparse matrix indices and values
+    rows = Vector{Int}(undef, max_elem)
+    cols = Vector{Int}(undef, max_elem)
+    vals = Vector{Int}(undef, max_elem)
+    idx = 1 # Tracks position in arrays
+    @inbounds for j in 1:n
+        mask = 1 << (h.r - 1) # Initialize mask for MSB
+        @simd for i in 1:h.r
+            if j & mask != 0 # Check if the i-th bit is 1
+                rows[idx] = i
+                cols[idx] = j
+                vals[idx] = 1
+                idx += 1
+            end
+            mask >>= 1 # Shift mask to next bit
+        end
+    end
+    # Resize arrays to actual number of non-zero elements
+    rows = rows[1:idx-1]
+    cols = cols[1:idx-1]
+    vals = vals[1:idx-1]
+    return sparse(rows, cols, vals, h.r, n)
+end
+
+code_n(h::Hamming) = 2 ^ h.r - 1
+
+code_k(h::Hamming) = 2 ^ h.r - 1 - h.r
+
+distance(h::Hamming) = 3

test/test_ecc_hamming.jl

@@ -0,0 +1,45 @@
+@testitem "ECC Hamming" begin
+
+    using LinearAlgebra
+    using QuantumClifford
+    using QuantumClifford.ECC
+    using QuantumClifford.ECC: AbstractECC, Hamming
+    using Nemo: matrix, GF, echelon_form
+
+    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 Hamming codes properties" begin
+        for r in 3:20
+            n = 2 ^ r - 1
+            k = 2 ^ r - 1 - r
+            H = parity_checks(Hamming(r))
+            H = Matrix{Bool}(H)
+            mat = matrix(GF(2), parity_checks(Hamming(r)))
+            computed_rank = rank(mat)
+            @test computed_rank == n - k
+        end
+        # minimum distance test is expensive (NP-hard) so we test for r = [3,4]
+        r_vals = [3,4]
+        for r in r_vals
+            H = parity_checks(Hamming(r))
+            H = Matrix{Bool}(H)
+            @test minimum_distance(parity_checks(Hamming(r))) == 3
+        end
+        # Example taken from [huffman2010fundamentals](@cite).
+        @test Matrix{Bool}(parity_checks(Hamming(3))) == [0  0  0  1  1  1  1;
+                                                          0  1  1  0  0  1  1;
+                                                          1  0  1  0  1  0  1]
+    end
+end

test/test_ecc_throws.jl

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