Began work on hamming codes - Implemented rudimentary version in GF(2)

docs/_static/extra.css

@@ -16,6 +16,14 @@
   background-color: #344FA1;
 }
 
+/* @media screen and (max-width: 76.1875em) */
+.md-nav--primary .md-nav__title[for=__drawer] {
+    background-color: #344FA1;
+}
+.md-source {
+  background-color: #344FA1;
+}
+
 /* Fix line numbers in code blocks */
 .linenos {
   background-color: var(--md-default-bg-color--light);

docs/conf.py

@@ -281,13 +281,13 @@
         "name": "fast-performance",
         "title": "Faster performance",
         "icon": "material/speedometer",
-        "color": (40, 167, 69),  # Green: --sd-color-success
+        "color": (47, 177, 112),  # Green: --md-code-hl-string-color
     },
     {
         "name": "slow-performance",
         "title": "Slower performance",
         "icon": "material/speedometer-slow",
-        "color": (220, 53, 69),  # Red: --sd-color-danger
+        "color": (230, 105, 91),  # Red: --md-code-hl-number-color
     },
 ]
 

src/galois/_codes/__init__.py

@@ -2,4 +2,5 @@
 A subpackage containing forward error correction codes.
 """
 from ._bch import *
+from ._hamming import *
 from ._reed_solomon import *

src/galois/_codes/_hamming.py

@@ -0,0 +1,195 @@
+"""
+G is kxn
+H is (n-k)xn
+m is 1xk
+c is 1xn
+"""
+from __future__ import annotations
+
+from typing import overload
+
+import numpy as np
+from typing_extensions import Literal
+
+from .._fields import GF, FieldArray
+from .._helper import export, verify_isinstance, verify_issubclass
+from ..typing import ArrayLike
+from ._linear import _LinearCode, parity_check_to_generator_matrix
+
+
+@export
+class Hamming(_LinearCode):
+
+    def __init__(self, n: int | None = None, k: int | None = None, field: FieldArray = GF(2), r: int | None = None, extended: bool = False):
+        # Construct the Generator and Parity Check Matrix and then initialise the parent class
+
+        if n is not None and k is not None:
+            r = n - k
+            q = field.characteristic
+            expected_n = int((q ** r - 1)/(q - 1))
+            if n != expected_n:
+                r = n - k - 1
+                expected_n_ext = int((q ** r - 1)/(q - 1)) + 1
+                if n != expected_n_ext:
+                    raise Exception(
+                        f"Given n = {n}, k = {k} but expected n = {expected_n} or {expected_n_ext}")
+                else:
+                    extended = True
+        elif r is not None:
+            q = field.characteristic
+            n = int((q ** r - 1)/(q - 1))
+            k = n - r
+            if extended:
+                n += 1
+
+        self._r = r
+        self._q = field.characteristic
+        self._extended = extended
+
+        k = int((self._q ** self._r - 1)/(self._q - 1)) - r
+        n = k + r
+        if extended:
+            n += 1
+        d = 4 if extended else 3
+
+        H = self._generate_systematic_parity_check_matrix(r, field, extended)
+        G = parity_check_to_generator_matrix(H)
+
+        super().__init__(n, k, d, G, H, systematic=True)
+
+    @property
+    def r(self) -> int:
+        return self._r
+
+    @property
+    def q(self) -> int:
+        return self._q
+
+    @property
+    def extended(self) -> bool:
+        return self._extended
+
+    def encode(self, message: ArrayLike, output: Literal["codeword", "parity"] = "codeword") -> FieldArray:
+        # Call the parent class's encode method
+        return super().encode(message=message, output=output)
+
+    def detect(self, codeword: FieldArray) -> np.ndarray:
+        # Call the parent class's detect method
+        return super().detect(codeword=codeword)
+
+    @overload
+    def decode(
+        self,
+        codeword: ArrayLike,
+        output: Literal["message", "codeword"] = "message",
+        errors: Literal[False] = False,
+    ) -> FieldArray:
+        ...
+
+    @overload
+    def decode(
+        self,
+        codeword: ArrayLike,
+        output: Literal["message", "codeword"] = "message",
+        errors: Literal[True] = True,
+    ) -> tuple[FieldArray, int | np.ndarray]:
+        ...
+
+    def decode(self, codeword, output="message", errors=False):
+        # Call the parent class's decode method, but implement the _decode_codeword method
+        return super().decode(codeword=codeword, output=output, errors=errors)
+
+    def _decode_codeword(self, codeword: FieldArray) -> tuple[FieldArray, np.ndarray]:
+        # Using syndrome error correction we will find the position of the error and then fix it
+        n_errors = np.zeros(codeword.shape[0], dtype=int)
+        syndromes = codeword @ self.H.T
+        errors_detected = ~np.all(syndromes == 0, axis=1)
+
+        # We will return this variable finally after making corrections in it
+        decoded_codeword = codeword.copy()
+        # We iterate over every codeword's syndrome
+        for i, (error_detected, syndrome) in enumerate(zip(errors_detected, syndromes)):
+            if not error_detected:
+                continue
+            # If an error is detected, then we find the position of the error
+            # Since Hamming codes are single error correcting codes
+            # Hc = H*(v + e) = Hv + He = He
+            # Here c is corrupted codeword to which an error vector e is added
+            # Here e is a weight-1 vector
+            # So He(syndrome) is a column of the parity check matrix scaled by a constant
+            # The location of this column is the position where the error occurred
+            error_position = 0
+            constant_scale = 0
+            while error_position < self.n:
+                parity_column = self.H[:, error_position]
+                for a in self._field.elements:
+                    if np.all(syndrome == a * parity_column):
+                        constant_scale = a
+                        break
+                if constant_scale != 0:
+                    break
+                error_position += 1
+            if error_position < self.n:
+                decoded_codeword[i, error_position] -= constant_scale
+                n_errors[i] = 1
+            else:
+                n_errors[i] = -1
+
+        return (decoded_codeword, n_errors)
+
+    def _convert_codeword_to_message(self, codeword: FieldArray) -> FieldArray:
+        message = None
+        if self.is_systematic:
+            message = codeword[..., :self.k]
+        return message
+
+    def _convert_codeword_to_parity(self, codeword: FieldArray) -> FieldArray:
+        parity = None
+        if self.is_systematic:
+            parity = codeword[..., -(self.n - self.k):]
+        return parity
+
+    ##############################################
+    # Helper Functions
+    ##############################################
+
+    @staticmethod
+    def _generate_systematic_parity_check_matrix(r, field, extended):
+
+        q = field.characteristic
+        n = int((q ** r - 1)/(q - 1))
+        H = field(np.zeros((r,  n), dtype=int))
+
+        gf = GF(q**r)
+
+        # Add the parity columns first
+        col_idx = 0
+        for num in gf.elements:
+            if num == 0:
+                continue
+            vec = gf(num).vector()
+            vec_weight = np.count_nonzero(vec)
+            # If the vector is a weight-1 vector, then it will be added at the end
+            # in the identity matrix
+            vec_exists = True if vec_weight == 1 else False
+            for factor in range(1, q):
+                scaled_vec = vec * factor
+                if gf.Vector(scaled_vec) < num:
+                    vec_exists = True
+            if not vec_exists:
+                H[:, col_idx] = vec
+                col_idx += 1
+
+        # Add the identity matrix
+        for pow in range(r-1, -1, -1):
+            H[:, col_idx] = gf(q**pow).vector()
+            col_idx += 1
+
+        if extended:
+            # Concatenate a zeros column to the right of the matrix
+            H = np.concatenate((H, np.zeros((r, 1), dtype=int)), axis=1)
+            # Concatenate a ones row to the bottom of the matrix
+            H = np.concatenate((H, np.ones((1, n + 1), dtype=int)), axis=0)
+            H = H.row_reduce(eye="right")
+
+        return H

src/galois/_polys/_dense.py

@@ -344,7 +344,7 @@ def __call__(self, a: Array, b: int, c: Array | None = None) -> Array:
         # Convert the integer b into a vector of uint64 [MSWord, ..., LSWord] so arbitrarily-large exponents may be
         # passed into the JIT-compiled version
         b_vec = []  # Pop on LSWord -> MSWord
-        while b > 2**64:
+        while b >= 2**64:
             q, r = divmod(b, 2**64)
             b_vec.append(r)
             b = q

tests/polys/test_arithmetic.py

@@ -175,6 +175,8 @@ def test_modular_power_large_exponent_jit():
         g_coeffs = [193, 88, 107, 214, 72, 3]
         f = R([F.fetch_int(fi) for fi in f_coeffs[::-1]])
         g = R([F.fetch_int(gi) for gi in g_coeffs[::-1]])
+        print(pow(f, 2**64 + 0, g))
+        print(pow(f, 2**64 + 1234, g))
         print(pow(f, 2**70 + 0, g))
         print(pow(f, 2**70 + 1234, g))
         print(pow(f, 2**70 + 7654, g))
@@ -184,6 +186,9 @@ def test_modular_power_large_exponent_jit():
     f = galois.Poly([255, 228, 34, 121, 243, 189, 6, 131, 102, 168, 82], field=GF)
     g = galois.Poly([193, 88, 107, 214, 72, 3], field=GF)
 
+    assert pow(f, 2**64 + 0, g) == galois.Poly.Str("221*x^4 + 112*x^3 + 124*x^2 + 87*x + 4", field=GF)
+    assert pow(f, 2**64 + 1234, g) == galois.Poly.Str("80*x^4 + 101*x^3 + 17*x^2 + 243*x + 74", field=GF)
+
     assert pow(f, 2**70 + 0, g) == galois.Poly.Str("178*x^4 + 228*x^3 + 198*x^2 + 191*x + 211", field=GF)
     assert pow(f, 2**70 + 1234, g) == galois.Poly.Str("100*x^4 + 242*x^3 + 235*x^2 + 171*x + 43", field=GF)
     assert pow(f, 2**70 + 7654, g) == galois.Poly.Str("203*x^4 + 203*x^3 + 155*x^2 + 221*x + 151", field=GF)
@@ -199,6 +204,8 @@ def test_modular_power_large_exponent_python():
         g_coeffs = [193, 88, 107, 214, 72, 3]
         f = R([F.fetch_int(fi) for fi in f_coeffs[::-1]])
         g = R([F.fetch_int(gi) for gi in g_coeffs[::-1]])
+        print([fi.integer_representation() for fi in pow(f, 2**64 + 0, g).list()[::-1]])
+        print([fi.integer_representation() for fi in pow(f, 2**64 + 1234, g).list()[::-1]])
         print([fi.integer_representation() for fi in pow(f, 2**70 + 0, g).list()[::-1]])
         print([fi.integer_representation() for fi in pow(f, 2**70 + 1234, g).list()[::-1]])
         print([fi.integer_representation() for fi in pow(f, 2**70 + 7654, g).list()[::-1]])
@@ -208,6 +215,27 @@ def test_modular_power_large_exponent_python():
     f = galois.Poly([255, 228, 34, 121, 243, 189, 6, 131, 102, 168, 82], field=GF)
     g = galois.Poly([193, 88, 107, 214, 72, 3], field=GF)
 
+    assert pow(f, 2**64 + 0, g) == galois.Poly(
+        [
+            441863720507004970561004876820,
+            13672604657018030323534797793,
+            667292756625776662634439753318,
+            1088336758052065395414982737453,
+            594859692808903855438703377448,
+        ],
+        field=GF,
+    )
+    assert pow(f, 2**64 + 1234, g) == galois.Poly(
+        [
+            199499378257056846135676738371,
+            461241799178212820931813246616,
+            997379144596106156826875511246,
+            1102132140575827285204706059544,
+            949801160766724443620356316946,
+        ],
+        field=GF,
+    )
+
     assert pow(f, 2**70 + 0, g) == galois.Poly(
         [
             420013998870488935594333531316,