Add initial LFSR implementation

learn2grid · Jul 30, 2021 · 8254a2c · 8254a2c
1 parent 1cc5c35
commit 8254a2c
Show file tree

Hide file tree

Showing 6 changed files with 523 additions and 125 deletions.
diff --git a/docs/api/linear-sequences.rst b/docs/api/linear-sequences.rst
@@ -5,6 +5,16 @@ This section contains classes and functions for creating and analyzing linear se
 
 .. currentmodule:: galois
 
+Linear-feedback shift registers
+-------------------------------
+
+.. rubric::
+.. autosummary::
+   :template: class.rst
+   :toctree:
+
+   LFSR
+
 Sequence analysis functions
 ---------------------------
 

diff --git a/galois/_code/_bch.py b/galois/_code/_bch.py
@@ -4,9 +4,10 @@
 from numba import int64
 import numpy as np
 
+from .. import _lfsr
 from .._factor import factors
 from .._field import Field, Poly, GF2, matlab_primitive_poly
-from .._field._meta_function import UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, BERLEKAMP_MASSEY_CALCULATE_SIG
+from .._field._meta_function import UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG
 from .._overrides import set_module
 
 from ._cyclic import poly_to_generator_matrix, roots_to_parity_check_matrix
@@ -216,7 +217,7 @@ def __init__(self, *args, **kwargs):  # pylint: disable=unused-argument
         self._power_jit = self.field._calculate_jit("power")
 
         # Pre-compile the JIT functions
-        self._berlekamp_massey_jit = self.field._function("berlekamp_massey")
+        self._berlekamp_massey_jit = _lfsr.jit_calculate("berlekamp_massey")
         self._poly_roots_jit = self.field._function("poly_roots")
         self._poly_divmod_jit = GF2._function("poly_divmod")
 
@@ -679,7 +680,7 @@ def is_narrow_sense(self):
 # JIT-compiled implementation of the specified functions
 ###############################################################################
 
-DECODE_CALCULATE_SIG = numba.types.FunctionType(int64[:,:](int64[:,:], int64[:,:], int64, int64, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BERLEKAMP_MASSEY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, int64, int64, int64))
+DECODE_CALCULATE_SIG = numba.types.FunctionType(int64[:,:](int64[:,:], int64[:,:], int64, int64, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, _lfsr.BERLEKAMP_MASSEY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, int64, int64, int64))
 
 def decode_calculate(codeword, syndrome, t, primitive_element, ADD, SUBTRACT, MULTIPLY, RECIPROCAL, POWER, BERLEKAMP_MASSEY, POLY_ROOTS, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY):  # pragma: no cover
     """

diff --git a/galois/_code/_reed_solomon.py b/galois/_code/_reed_solomon.py
@@ -2,9 +2,10 @@
 from numba import int64
 import numpy as np
 
+from .. import _lfsr
 from .._factor import factors
 from .._field import Field, Poly, matlab_primitive_poly
-from .._field._meta_function import UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BERLEKAMP_MASSEY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, POLY_EVALUATE_CALCULATE_SIG, CONVOLVE_CALCULATE_SIG
+from .._field._meta_function import UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, POLY_EVALUATE_CALCULATE_SIG, CONVOLVE_CALCULATE_SIG
 from .._overrides import set_module
 
 from ._cyclic import poly_to_generator_matrix, roots_to_parity_check_matrix
@@ -120,7 +121,7 @@ def __init__(self, *args, **kwargs):  # pylint: disable=unused-argument
         self._power_jit = self.field._calculate_jit("power")
 
         # Pre-compile the JIT functions
-        self._berlekamp_massey_jit = self.field._function("berlekamp_massey")
+        self._berlekamp_massey_jit = _lfsr.jit_calculate("berlekamp_massey")
         self._poly_divmod_jit = self.field._function("poly_divmod")
         self._poly_roots_jit = self.field._function("poly_roots")
         self._poly_eval_jit = self.field._function("poly_evaluate")
@@ -587,7 +588,7 @@ def is_narrow_sense(self):
 # JIT-compiled implementation of the specified functions
 ###############################################################################
 
-DECODE_CALCULATE_SIG = numba.types.FunctionType(int64[:,:](int64[:,:], int64[:,:], int64, int64, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BERLEKAMP_MASSEY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, POLY_EVALUATE_CALCULATE_SIG, CONVOLVE_CALCULATE_SIG, int64, int64, int64))
+DECODE_CALCULATE_SIG = numba.types.FunctionType(int64[:,:](int64[:,:], int64[:,:], int64, int64, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, UNARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, _lfsr.BERLEKAMP_MASSEY_CALCULATE_SIG, POLY_ROOTS_CALCULATE_SIG, POLY_EVALUATE_CALCULATE_SIG, CONVOLVE_CALCULATE_SIG, int64, int64, int64))
 
 def decode_calculate(codeword, syndrome, t, primitive_element, ADD, SUBTRACT, MULTIPLY, RECIPROCAL, POWER, BERLEKAMP_MASSEY, POLY_ROOTS, POLY_EVAL, CONVOLVE, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY):  # pragma: no cover
     """

diff --git a/galois/_field/_meta_function.py b/galois/_field/_meta_function.py
@@ -295,30 +295,6 @@ def _poly_roots(cls, nonzero_degrees, nonzero_coeffs):
         idxs = np.argsort(roots)
         return roots[idxs]
 
-    def _berlekamp_massey(cls, sequence):
-        assert isinstance(sequence, cls)
-        field = cls
-        dtype = sequence.dtype
-
-        if cls._ufunc_mode != "python-calculate":
-            sequence = sequence.astype(np.int64)
-            add = cls._calculate_jit("add")
-            subtract = cls._calculate_jit("subtract")
-            multiply = cls._calculate_jit("multiply")
-            reciprocal = cls._calculate_jit("reciprocal")
-            coeffs = cls._function("berlekamp_massey")(sequence, add, subtract, multiply, reciprocal, cls.characteristic, cls.degree, cls._irreducible_poly_int)
-            coeffs = coeffs.astype(dtype)
-        else:
-            sequence = sequence.view(np.ndarray)
-            subtract = cls._python_func("subtract")
-            multiply = cls._python_func("multiply")
-            DIVIDE = cls._python_func("divide")
-            reciprocal = cls._python_func("reciprocal")
-            coeffs = cls._function("berlekamp_massey")(sequence, subtract, multiply, DIVIDE, reciprocal, cls.characteristic, cls.degree, cls._irreducible_poly_int)
-        coeffs = coeffs.view(field)
-
-        return coeffs
-
 
 ##############################################################################
 # Individual gufuncs
@@ -638,98 +614,9 @@ def poly_evaluate_calculate(coeffs, value, ADD, MULTIPLY, CHARACTERISTIC, DEGREE
     return result
 
 
-###############################################################################
-# Berlekamp-Massey algorithm
-###############################################################################
-
-BERLEKAMP_MASSEY_LOOKUP_SIG = numba.types.FunctionType(int64[:](int64[:], BINARY_LOOKUP_SIG, BINARY_LOOKUP_SIG, BINARY_LOOKUP_SIG, UNARY_LOOKUP_SIG, int64[:], int64[:], int64[:], int64))
-
-def berlekamp_massey_lookup(sequence, ADD, SUBTRACT, MULTIPLY, RECIPROCAL, EXP, LOG, ZECH_LOG, ZECH_E):  # pragma: no cover
-    args = EXP, LOG, ZECH_LOG, ZECH_E
-    dtype = sequence.dtype
-
-    N = sequence.size
-
-    s = sequence
-    c = np.zeros(N, dtype=dtype)
-    b = np.zeros(N, dtype=dtype)
-    c[0] = 1  # The polynomial c(x) = 1
-    b[0] = 1  # The polynomial b(x) = 1
-    L = 0
-    m = 1
-    bb = 1
-
-    for n in range(0, N):
-        d = 0
-        for i in range(0, L + 1):
-            d = ADD(d, MULTIPLY(s[n - i], c[i], *args), *args)
-
-        if d == 0:
-            m += 1
-        elif 2*L <= n:
-            t = np.copy(c)
-            d_bb = MULTIPLY(d, RECIPROCAL(bb, *args), *args)
-            for i in range(m, N):
-                c[i] = SUBTRACT(c[i], MULTIPLY(d_bb, b[i - m], *args), *args)
-            L = n + 1 - L
-            b = np.copy(t)
-            bb = d
-            m = 1
-        else:
-            d_bb = MULTIPLY(d, RECIPROCAL(bb, *args), *args)
-            for i in range(m, N):
-                c[i] = SUBTRACT(c[i], MULTIPLY(d_bb, b[i - m], *args), *args)
-            m += 1
-
-    return c[0:L + 1][::-1]
-
-
-BERLEKAMP_MASSEY_CALCULATE_SIG = numba.types.FunctionType(int64[:](int64[:], BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, BINARY_CALCULATE_SIG, UNARY_CALCULATE_SIG, int64, int64, int64))
-
-def berlekamp_massey_calculate(sequence, ADD, SUBTRACT, MULTIPLY, RECIPROCAL, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY):  # pragma: no cover
-    args = CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY
-    dtype = sequence.dtype
-
-    N = sequence.size
-
-    s = sequence
-    c = np.zeros(N, dtype=dtype)
-    b = np.zeros(N, dtype=dtype)
-    c[0] = 1  # The polynomial c(x) = 1
-    b[0] = 1  # The polynomial b(x) = 1
-    L = 0
-    m = 1
-    bb = 1
-
-    for n in range(0, N):
-        d = 0
-        for i in range(0, L + 1):
-            d = ADD(d, MULTIPLY(s[n - i], c[i], *args), *args)
-
-        if d == 0:
-            m += 1
-        elif 2*L <= n:
-            t = np.copy(c)
-            d_bb = MULTIPLY(d, RECIPROCAL(bb, *args), *args)
-            for i in range(m, N):
-                c[i] = SUBTRACT(c[i], MULTIPLY(d_bb, b[i - m], *args), *args)
-            L = n + 1 - L
-            b = np.copy(t)
-            bb = d
-            m = 1
-        else:
-            d_bb = MULTIPLY(d, RECIPROCAL(bb, *args), *args)
-            for i in range(m, N):
-                c[i] = SUBTRACT(c[i], MULTIPLY(d_bb, b[i - m], *args), *args)
-            m += 1
-
-    return c[0:L + 1][::-1]
-
-
 # Set globals so that when the application-specific python functions are called they will use the
 # pure-python arithmetic routines
 MATMUL = matmul_calculate
 CONVOLVE = convolve_calculate
 POLY_DIVMOD = poly_divmod_calculate
 POLY_ROOTS = poly_roots_calculate
-BERLEKAMP_MASSEY = berlekamp_massey_calculate