From fa0038873d28307723467b5be3c62e23b44e5f22 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 18 Aug 2021 22:23:03 -0400 Subject: [PATCH] Greatly clean-up Galois field arithmetic functions/classes --- galois/_codes/_bch.py | 10 +- galois/_codes/_reed_solomon.py | 10 +- galois/_fields/_array.py | 6 + galois/_fields/_calculate.py | 171 ++++++----- galois/_fields/_class.py | 28 +- galois/_fields/_functions.py | 61 ++-- galois/_fields/_gf2.py | 185 +++-------- galois/_fields/_gf2m.py | 402 ++++++++---------------- galois/_fields/_gfp.py | 381 +++++++---------------- galois/_fields/_gfpm.py | 544 +++++++++++---------------------- galois/_fields/_lookup.py | 539 ++++++++++++++------------------ galois/_fields/_properties.py | 9 +- galois/_fields/_ufuncs.py | 132 +------- galois/_lfsr.py | 31 +- 14 files changed, 871 insertions(+), 1638 deletions(-) diff --git a/galois/_codes/_bch.py b/galois/_codes/_bch.py index 4ef2e2a0ac..7416345a7b 100644 --- a/galois/_codes/_bch.py +++ b/galois/_codes/_bch.py @@ -231,11 +231,11 @@ def __new__(cls, n, k, primitive_poly=None, primitive_element=None, systematic=T def __init__(self, *args, **kwargs): # pylint: disable=unused-argument # Pre-compile the arithmetic methods - self._add_jit = self.field._calculate_jit("add") - self._subtract_jit = self.field._calculate_jit("subtract") - self._multiply_jit = self.field._calculate_jit("multiply") - self._reciprocal_jit = self.field._calculate_jit("reciprocal") - self._power_jit = self.field._calculate_jit("power") + self._add_jit = self.field._func_calculate("add") + self._subtract_jit = self.field._func_calculate("subtract") + self._multiply_jit = self.field._func_calculate("multiply") + self._reciprocal_jit = self.field._func_calculate("reciprocal") + self._power_jit = self.field._func_calculate("power") # Pre-compile the JIT functions self._berlekamp_massey_jit = _lfsr.jit_calculate("berlekamp_massey") diff --git a/galois/_codes/_reed_solomon.py b/galois/_codes/_reed_solomon.py index d5627fe86b..abd6fe7fdc 100644 --- a/galois/_codes/_reed_solomon.py +++ b/galois/_codes/_reed_solomon.py @@ -129,11 +129,11 @@ def __new__(cls, n, k, c=1, primitive_poly=None, primitive_element=None, systema def __init__(self, *args, **kwargs): # pylint: disable=unused-argument # Pre-compile the arithmetic methods - self._add_jit = self.field._calculate_jit("add") - self._subtract_jit = self.field._calculate_jit("subtract") - self._multiply_jit = self.field._calculate_jit("multiply") - self._reciprocal_jit = self.field._calculate_jit("reciprocal") - self._power_jit = self.field._calculate_jit("power") + self._add_jit = self.field._func_calculate("add") + self._subtract_jit = self.field._func_calculate("subtract") + self._multiply_jit = self.field._func_calculate("multiply") + self._reciprocal_jit = self.field._func_calculate("reciprocal") + self._power_jit = self.field._func_calculate("power") # Pre-compile the JIT functions self._berlekamp_massey_jit = _lfsr.jit_calculate("berlekamp_massey") diff --git a/galois/_fields/_array.py b/galois/_fields/_array.py index 484039662c..05227d328d 100644 --- a/galois/_fields/_array.py +++ b/galois/_fields/_array.py @@ -1076,6 +1076,12 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs): def __str__(self): return self.__repr__() + # formatter = type(self)._formatter(self) + + # with np.printoptions(formatter=formatter): + # string = super().__str__() + + # return string def __repr__(self): formatter = type(self)._formatter(self) diff --git a/galois/_fields/_calculate.py b/galois/_fields/_calculate.py index 54be5046de..5f9bf089ac 100644 --- a/galois/_fields/_calculate.py +++ b/galois/_fields/_calculate.py @@ -14,11 +14,13 @@ class CalculateMeta(PropertiesMeta): """ # pylint: disable=no-value-for-parameter - # Function signatures for JIT-compiled arithmetic functions + # Function signatures for JIT-compiled "calculate" arithmetic functions _UNARY_CALCULATE_SIG = numba.types.FunctionType(int64(int64, int64, int64, int64)) _BINARY_CALCULATE_SIG = numba.types.FunctionType(int64(int64, int64, int64, int64, int64)) - # Lookup table of ufuncs to unary/binary type needed for LookupMeta, CalculateMeta, etc + _FUNC_CACHE_CALCULATE = {} + _UFUNC_CACHE_CALCULATE = {} + _UFUNC_TYPE = { "add": "binary", "negative": "unary", @@ -30,102 +32,123 @@ class CalculateMeta(PropertiesMeta): "log": "binary", } - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### + def __init__(cls, name, bases, namespace, **kwargs): + super().__init__(name, bases, namespace, **kwargs) + cls._reset_globals() - def _calculate_jit(cls, name): + def _set_globals(cls, name): # pylint: disable=unused-argument,no-self-use """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. + Sets the global variables in the GF*Meta metaclass mixins that are needed for linking during JIT compilation. """ - raise NotImplementedError + return - def _python_func(cls, name): + def _reset_globals(cls): # pylint: disable=no-self-use """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. + Resets the global variables in the GF*Meta metaclass mixins so when the pure-python functions/ufuncs invoke these + globals, they reference the correct pure-python functions and not the JIT-compiled functions/ufuncs. """ - raise NotImplementedError + return - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### + def _func_calculate(cls, name, reset=True): + """ + Returns a JIT-compiled arithmetic function using explicit calculation. These functions are once-compiled and shared for all + Galois fields. The only difference between Galois fields are the characteristic, degree, and irreducible polynomial that are + passed in as inputs. + """ + key = (name,) - def _calculate_ufunc(cls, name): - raise NotImplementedError + if key not in cls._FUNC_CACHE_CALCULATE: + cls._set_globals(name) + function = getattr(cls, f"_{name}_calculate") - def _python_ufunc(cls, name): # pylint: disable=no-self-use - function = eval(f"cls._{name}_python") + if cls._UFUNC_TYPE[name] == "unary": + cls._FUNC_CACHE_CALCULATE[key] = numba.jit(["int64(int64, int64, int64, int64)"], nopython=True, cache=True)(function) + else: + cls._FUNC_CACHE_CALCULATE[key] = numba.jit(["int64(int64, int64, int64, int64, int64)"], nopython=True, cache=True)(function) + + if reset: + cls._reset_globals() + + return cls._FUNC_CACHE_CALCULATE[key] + + def _ufunc_calculate(cls, name): + """ + Returns a JIT-compiled arithmetic ufunc using explicit calculation. These ufuncs are compiled for each Galois field since + the characteristic, degree, and irreducible polynomial are compiled into the ufuncs as constants. + """ + key = (name, cls.characteristic, cls.degree, cls._irreducible_poly_int) + + if key not in cls._UFUNC_CACHE_CALCULATE: + cls._set_globals(name) + function = getattr(cls, f"_{name}_calculate") + + # These variables must be locals and not class properties for Numba to compile them as literals + CHARACTERISTIC = cls.characteristic + DEGREE = cls.degree + IRREDUCIBLE_POLY = cls._irreducible_poly_int + + if cls._UFUNC_TYPE[name] == "unary": + cls._UFUNC_CACHE_CALCULATE[key] = numba.vectorize(["int64(int64)"], nopython=True)(lambda a: function(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY)) + else: + cls._UFUNC_CACHE_CALCULATE[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(lambda a, b: function(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY)) + + cls._reset_globals() + + return cls._UFUNC_CACHE_CALCULATE[key] + + def _func_python(cls, name): + """ + Returns a pure-python arithmetic function using explicit calculation. This lambda function wraps the arithmetic functions in + GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta by passing in the field's characteristic, degree, and irreducible polynomial. + """ + return getattr(cls, f"_{name}_calculate") + # if cls._UFUNC_TYPE[name] == "unary": + # return lambda a: function(a, cls.characteristic, cls.degree, cls._irreducible_poly_int) + # else: + # return lambda a, b: function(a, b, cls.characteristic, cls.degree, cls._irreducible_poly_int) + + def _ufunc_python(cls, name): + """ + Returns a pure-python arithmetic ufunc using explicit calculation. + """ + function = getattr(cls, f"_{name}_calculate") if cls._UFUNC_TYPE[name] == "unary": - return np.frompyfunc(function, 1, 1) + return np.frompyfunc(lambda a: function(a, cls.characteristic, cls.degree, cls._irreducible_poly_int), 1, 1) else: - return np.frompyfunc(function, 2, 1) + return np.frompyfunc(lambda a, b: function(a, b, cls.characteristic, cls.degree, cls._irreducible_poly_int), 2, 1) ############################################################################### - # Pure-python arithmetic methods using explicit calculation + # Arithmetic functions using explicit calculation ############################################################################### - def _add_python(cls, a, b): + @staticmethod + def _add_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): raise NotImplementedError - def _negative_python(cls, a): + @staticmethod + def _negative_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): raise NotImplementedError - def _subtract_python(cls, a, b): + @staticmethod + def _subtract_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): raise NotImplementedError - def _multiply_python(cls, a, b): + @staticmethod + def _multiply_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): raise NotImplementedError - def _reciprocal_python(cls, a): + @staticmethod + def _reciprocal_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): raise NotImplementedError - def _divide_python(cls, a, b): - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if a == 0: - return 0 - else: - b_inv = cls._reciprocal_python(b) - return cls._multiply_python(a, b_inv) - - def _power_python(cls, a, b): - if a == 0 and b < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b == 0: - return 1 - elif b < 0: - a = cls._reciprocal_python(a) - b = abs(b) - - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while b > 1: - if b % 2 == 0: - result_s = cls._multiply_python(result_s, result_s) - b //= 2 - else: - result_m = cls._multiply_python(result_m, result_s) - b -= 1 - - result = cls._multiply_python(result_m, result_s) - - return result - - def _log_python(cls, a, b): - """ - TODO: Replace this with a more efficient algorithm - """ - if a == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") + @staticmethod + def _divide_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + raise NotImplementedError - # Naive algorithm - result = 1 - for i in range(0, cls.order - 1): - if result == a: - break - result = cls._multiply_python(result, b) + @staticmethod + def _power_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + raise NotImplementedError - return i + @staticmethod + def _log_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + raise NotImplementedError diff --git a/galois/_fields/_class.py b/galois/_fields/_class.py index 33e1abea75..69230612ef 100644 --- a/galois/_fields/_class.py +++ b/galois/_fields/_class.py @@ -23,7 +23,7 @@ class FieldClass(FunctionMeta, UfuncMeta, PropertiesMeta): This metaclass gives :obj:`galois.FieldArray` subclasses returned from the class factory :func:`galois.GF` methods and attributes related to its Galois field. """ - # pylint: disable=no-value-for-parameter + # pylint: disable=no-value-for-parameter,unsupported-membership-test,abstract-method def __new__(cls, name, bases, namespace, **kwargs): # pylint: disable=unused-argument return super().__new__(cls, name, bases, namespace) @@ -34,32 +34,6 @@ def __str__(cls): def __repr__(cls): return str(cls) - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _calculate_jit(cls, name): - """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. - """ - raise NotImplementedError - - def _python_func(cls, name): - """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. - """ - raise NotImplementedError - - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### - - def _calculate_ufunc(cls, name): - """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. - """ - raise NotImplementedError - ############################################################################### # Class methods ############################################################################### diff --git a/galois/_fields/_functions.py b/galois/_fields/_functions.py index a440b7e878..074d731970 100644 --- a/galois/_fields/_functions.py +++ b/galois/_fields/_functions.py @@ -9,6 +9,9 @@ from ._dtypes import DTYPES from ._ufuncs import UfuncMeta +ADD_JIT = lambda a, b, *args: a + b +MULTIPLY_JIT = lambda a, b, *args: a * b + class FunctionMeta(UfuncMeta): """ @@ -81,8 +84,8 @@ def _gufunc(cls, name): cls._gufuncs["poly_evaluate"] = np.vectorize(poly_evaluate_python_calculate, excluded=["coeffs"], otypes=[np.object_]) else: global ADD_JIT, MULTIPLY_JIT - ADD_JIT = cls._calculate_jit("add") - MULTIPLY_JIT = cls._calculate_jit("multiply") + ADD_JIT = cls._func_calculate("add") + MULTIPLY_JIT = cls._func_calculate("multiply") cls._gufuncs["poly_evaluate"] = numba.guvectorize("int64[:], int64[:], int64, int64, int64, int64[:]", "(m),(n),(),(),()->(n)", nopython=True)(poly_evaluate_gufunc_calculate) else: raise NotImplementedError @@ -113,8 +116,8 @@ def _poly_evaluate(cls, coeffs, x): if cls.ufunc_mode == "python-calculate": # For object dtypes, call the vectorized classmethod - add = cls._python_func("add") - multiply = cls._python_func("multiply") + add = cls._func_python("add") + multiply = cls._func_python("multiply") y = cls._gufunc("poly_evaluate")(coeffs=coeffs.view(np.ndarray), value=x.view(np.ndarray), ADD=add, MULTIPLY=multiply, CHARACTERISTIC=cls.characteristic, DEGREE=cls.degree, IRREDUCIBLE_POLY=cls._irreducible_poly_int, result=field.Zeros(x.shape)) else: # For integer dtypes, call the JIT-compiled gufunc @@ -129,8 +132,8 @@ def _poly_evaluate(cls, coeffs, x): def _poly_evaluate_python(cls, coeffs, x): assert coeffs.ndim == 1 - add = cls._python_func("add") - multiply = cls._python_func("multiply") + add = cls._func_python("add") + multiply = cls._func_python("multiply") coeffs = coeffs.view(np.ndarray).astype(cls.dtypes[-1]) x = int(x) y = poly_evaluate_python_calculate(coeffs, x, add, multiply, cls.characteristic, cls.degree, cls._irreducible_poly_int, 0) @@ -175,15 +178,15 @@ def _matmul(cls, A, B, out=None, **kwargs): # pylint: disable=unused-argument if cls.ufunc_mode != "python-calculate": A = A.astype(np.int64) B = B.astype(np.int64) - add = cls._calculate_jit("add") - multiply = cls._calculate_jit("multiply") + add = cls._func_calculate("add") + multiply = cls._func_calculate("multiply") C = cls._function("matmul")(A, B, add, multiply, cls.characteristic, cls.degree, cls._irreducible_poly_int) C = C.astype(dtype) else: A = A.view(np.ndarray) B = B.view(np.ndarray) - add = cls._python_func("add") - multiply = cls._python_func("multiply") + add = cls._func_python("add") + multiply = cls._func_python("multiply") C = cls._function("matmul")(A, B, add, multiply, cls.characteristic, cls.degree, cls._irreducible_poly_int) C = C.view(field) @@ -227,15 +230,15 @@ def _convolve(cls, a, b, mode="full"): if cls._ufunc_mode != "python-calculate": a = a.astype(np.int64) b = b.astype(np.int64) - add = cls._calculate_jit("add") - multiply = cls._calculate_jit("multiply") + add = cls._func_calculate("add") + multiply = cls._func_calculate("multiply") c = cls._function("convolve")(a, b, add, multiply, cls.characteristic, cls.degree, cls._irreducible_poly_int) c = c.astype(dtype) else: a = a.view(np.ndarray) b = b.view(np.ndarray) - add = cls._python_func("add") - multiply = cls._python_func("multiply") + add = cls._func_python("add") + multiply = cls._func_python("multiply") c = cls._function("convolve")(a, b, add, multiply, cls.characteristic, cls.degree, cls._irreducible_poly_int) c = c.view(field) @@ -255,17 +258,17 @@ def _poly_divmod(cls, a, b): if cls._ufunc_mode != "python-calculate": a = a.astype(np.int64) b = b.astype(np.int64) - subtract = cls._calculate_jit("subtract") - multiply = cls._calculate_jit("multiply") - divide = cls._calculate_jit("divide") + subtract = cls._func_calculate("subtract") + multiply = cls._func_calculate("multiply") + divide = cls._func_calculate("divide") qr = cls._function("poly_divmod")(a, b, subtract, multiply, divide, cls.characteristic, cls.degree, cls._irreducible_poly_int) qr = qr.astype(dtype) else: a = a.view(np.ndarray) b = b.view(np.ndarray) - subtract = cls._python_func("subtract") - multiply = cls._python_func("multiply") - divide = cls._python_func("divide") + subtract = cls._func_python("subtract") + multiply = cls._func_python("multiply") + divide = cls._func_python("divide") qr = cls._function("poly_divmod")(a, b, subtract, multiply, divide, cls.characteristic, cls.degree, cls._irreducible_poly_int) qr = qr.view(field) @@ -286,17 +289,17 @@ def _poly_roots(cls, nonzero_degrees, nonzero_coeffs): if cls._ufunc_mode != "python-calculate": nonzero_degrees = nonzero_degrees.astype(np.int64) nonzero_coeffs = nonzero_coeffs.astype(np.int64) - add = cls._calculate_jit("add") - multiply = cls._calculate_jit("multiply") - power = cls._calculate_jit("power") + add = cls._func_calculate("add") + multiply = cls._func_calculate("multiply") + power = cls._func_calculate("power") roots = cls._function("poly_roots")(nonzero_degrees, nonzero_coeffs, np.int64(cls.primitive_element), add, multiply, power, cls.characteristic, cls.degree, cls._irreducible_poly_int)[0,:] roots = roots.astype(dtype) else: nonzero_degrees = nonzero_degrees.view(np.ndarray) nonzero_coeffs = nonzero_coeffs.view(np.ndarray) - add = cls._python_func("add") - multiply = cls._python_func("multiply") - power = cls._python_func("power") + add = cls._func_python("add") + multiply = cls._func_python("multiply") + power = cls._func_python("power") roots = cls._function("poly_roots")(nonzero_degrees, nonzero_coeffs, int(cls.primitive_element), add, multiply, power, cls.characteristic, cls.degree, cls._irreducible_poly_int)[0,:] roots = roots.view(field) @@ -308,12 +311,6 @@ def _poly_roots(cls, nonzero_degrees, nonzero_coeffs): # Individual gufuncs ############################################################################## -ADD_JIT = lambda a, b, *args: a + b -MULTIPLY_JIT = lambda a, b, *args: a * b - -# pylint: disable=redefined-outer-name,unused-argument - - def poly_evaluate_gufunc_lookup(coeffs, values, ADD, MULTIPLY, EXP, LOG, ZECH_LOG, ZECH_E, results): # pragma: no cover args = EXP, LOG, ZECH_LOG, ZECH_E diff --git a/galois/_fields/_gf2.py b/galois/_fields/_gf2.py index af247ee403..92291eada0 100644 --- a/galois/_fields/_gf2.py +++ b/galois/_fields/_gf2.py @@ -19,6 +19,9 @@ class GF2Meta(FieldClass, DirMeta, CalculateMeta): """ # pylint: disable=no-value-for-parameter + # Need to have a unique cache of "calculate" function for GF(2) + _FUNC_CACHE_CALCULATE = {} + def __init__(cls, name, bases, namespace, **kwargs): super().__init__(name, bases, namespace, **kwargs) cls._prime_subfield = cls @@ -45,23 +48,6 @@ def _compile_ufuncs(cls): cls._ufuncs["reciprocal"] = np.positive cls._ufuncs["divide"] = np.bitwise_and - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _calculate_jit(cls, name): - return compile_jit(name) - - def _python_func(cls, name): - return eval(f"{name}") - - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### - - def _calculate_ufunc(cls, name): - return compile_ufunc(name, cls.characteristic, cls.degree, cls._irreducible_poly_int) - ############################################################################### # Override ufunc routines to use native numpy bitwise ufuncs for GF(2) # arithmetic, which is faster than custom ufuncs @@ -101,34 +87,54 @@ def _ufunc_routine_square(cls, ufunc, method, inputs, kwargs, meta): # pylint: return inputs[0] ############################################################################### - # Pure python arithmetic methods + # Arithmetic functions using explicit calculation ############################################################################### - def _add_python(cls, a, b): + @staticmethod + def _add_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.bitwise_xor()` is faster. + """ return a ^ b - def _negative_python(cls, a): + @staticmethod + def _negative_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.positive()` is faster. + """ return a - def _subtract_python(cls, a, b): + @staticmethod + def _subtract_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.bitwise_xor()` is faster. + """ return a ^ b - def _multiply_python(cls, a, b): + @staticmethod + def _multiply_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.bitwise_and()` is faster. + """ return a & b - def _reciprocal_python(cls, a): + @staticmethod + def _reciprocal_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): if a == 0: raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") return 1 - def _divide_python(cls, a, b): + @staticmethod + def _divide_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): if b == 0: raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") return a & b - def _power_python(cls, a, b): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _power_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): if a == 0 and b < 0: raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") @@ -137,7 +143,9 @@ def _power_python(cls, a, b): else: return a - def _log_python(cls, a, b): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _log_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): if a == 0: raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") if b != 1: @@ -146,131 +154,6 @@ def _log_python(cls, a, b): return 0 -############################################################################### -# Compile functions -############################################################################### - -CHARACTERISTIC = None # The prime characteristic `p` of the Galois field -DEGREE = None # The prime power `m` of the Galois field -IRREDUCIBLE_POLY = None # The field's primitive polynomial in integer form - -# pylint: disable=redefined-outer-name,unused-argument - - -def compile_jit(name): - """ - Compile a JIT arithmetic function. These can be cached. - """ - if name not in compile_jit.cache: - function = eval(f"{name}") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64)"], nopython=True, cache=True)(function) - else: - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64, int64)"], nopython=True, cache=True)(function) - return compile_jit.cache[name] - -compile_jit.cache = {} - - -def compile_ufunc(name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_): - """ - Compile an arithmetic ufunc. These cannot be cached as the field parameters are compiled into the binary. - """ - key = (name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_) - if key not in compile_ufunc.cache: - global CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY - CHARACTERISTIC = CHARACTERISTIC_ - DEGREE = DEGREE_ - IRREDUCIBLE_POLY = IRREDUCIBLE_POLY_ - - function = eval(f"{name}_ufunc") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_ufunc.cache[key] = numba.vectorize(["int64(int64)"], nopython=True)(function) - else: - compile_ufunc.cache[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(function) - - return compile_ufunc.cache[key] - -compile_ufunc.cache = {} - - -############################################################################### -# Arithmetic explicitly calculated -############################################################################### - -def add(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a ^ b - - -def negative(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a - - -def subtract(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a ^ b - - -def multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a & b - - -@numba.extending.register_jitable(inline="always") -def reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if a == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - return 1 - - -def reciprocal_ufunc(a): # pragma: no cover - return reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - return a & b - - -def divide_ufunc(a, b): # pragma: no cover - return divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if a == 0 and b < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b == 0: - return 1 - else: - return a - - -def power_ufunc(a, b): # pragma: no cover - return power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if a == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") - if b != 1: - raise ArithmeticError("In GF(2), 1 is the only multiplicative generator.") - - return 0 - - -def log_ufunc(a, b): # pragma: no cover - return log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -############################################################################### -# A pre-generated FieldArray subclass for GF(2) -############################################################################### - @set_module("galois") class GF2(FieldArray, metaclass=GF2Meta, characteristic=2, degree=1, order=2, primitive_element=1, compile="jit-calculate"): r""" diff --git a/galois/_fields/_gf2m.py b/galois/_fields/_gf2m.py index d2657776ca..f06270a31a 100644 --- a/galois/_fields/_gf2m.py +++ b/galois/_fields/_gf2m.py @@ -6,6 +6,9 @@ from ._class import FieldClass, DirMeta +MULTIPLY = lambda a, b, *args: a * b +RECIPROCAL = lambda a, *args: 1 / a + class GF2mMeta(FieldClass, DirMeta): """ @@ -13,6 +16,9 @@ class GF2mMeta(FieldClass, DirMeta): """ # pylint: disable=no-value-for-parameter + # Need to have a unique cache of "calculate" function for GF(2^m) + _FUNC_CACHE_CALCULATE = {} + def __init__(cls, name, bases, namespace, **kwargs): super().__init__(name, bases, namespace, **kwargs) cls._prime_subfield = kwargs["prime_subfield"] @@ -24,46 +30,65 @@ def __init__(cls, name, bases, namespace, **kwargs): cls._is_primitive_poly = cls._poly_evaluate_python(cls._irreducible_poly.coeffs, cls.primitive_element) == 0 def _compile_ufuncs(cls): - cls._ufuncs = {} # Reset the dictionary so each ufunc will get recompiled - if cls.ufunc_mode == "jit-lookup": - cls._build_lookup_tables() + super()._compile_ufuncs() # Some explicit calculation functions are faster than using lookup tables. See https://github.com/mhostetter/galois/pull/92#issuecomment-835552639. cls._ufuncs["add"] = np.bitwise_xor cls._ufuncs["negative"] = np.positive cls._ufuncs["subtract"] = np.bitwise_xor - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _calculate_jit(cls, name): - return compile_jit(name) - - def _python_func(cls, name): - return eval(f"{name}") - - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### + def _set_globals(cls, name): + global MULTIPLY, RECIPROCAL + if name in ["reciprocal", "divide", "power", "log"]: + MULTIPLY = cls._func_calculate("multiply") + if name in ["divide", "power"]: + RECIPROCAL = cls._func_calculate("reciprocal") - def _calculate_ufunc(cls, name): - return compile_ufunc(name, cls.characteristic, cls.degree, cls._irreducible_poly_int) + def _reset_globals(cls): + global MULTIPLY, RECIPROCAL + MULTIPLY = cls._multiply_calculate + RECIPROCAL = cls._reciprocal_calculate ############################################################################### - # Pure python arithmetic methods + # Arithmetic functions using explicit calculation ############################################################################### - def _add_python(cls, a, b): + @staticmethod + def _add_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.bitwise_xor()` is faster. + """ return a ^ b - def _negative_python(cls, a): + @staticmethod + def _negative_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.positive()` is faster. + """ return a - def _subtract_python(cls, a, b): + @staticmethod + def _subtract_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Not actually used. `np.bitwise_xor()` is faster. + """ return a ^ b - def _multiply_python(cls, a, b): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _multiply_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + a in GF(2^m), can be represented as a degree m-1 polynomial a(x) in GF(2)[x] + b in GF(2^m), can be represented as a degree m-1 polynomial b(x) in GF(2)[x] + p(x) in GF(2)[x] with degree m is the irreducible polynomial of GF(2^m) + + a * b = c + = (a(x) * b(x)) % p(x) in GF(2) + = c(x) + = c + """ + ORDER = CHARACTERISTIC**DEGREE + # Re-order operands such that a > b so the while loop has less loops if b > a: a, b = b, a @@ -75,280 +100,113 @@ def _multiply_python(cls, a, b): b >>= 1 # Divide b(x) by x a <<= 1 # Multiply a(x) by x - if a >= cls._order: - a ^= cls._irreducible_poly_int # Compute a(x) % p(x) + if a >= ORDER: + a ^= IRREDUCIBLE_POLY # Compute a(x) % p(x) return c - def _reciprocal_python(cls, a): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _reciprocal_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + From Fermat's Little Theorem: + a in GF(p^m) + a^(p^m - 1) = 1 + + a * a^-1 = 1 + a * a^-1 = a^(p^m - 1) + a^-1 = a^(p^m - 2) + """ if a == 0: raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - exponent = cls.order - 2 + ORDER = CHARACTERISTIC**DEGREE + exponent = ORDER - 2 result_s = a # The "squaring" part result_m = 1 # The "multiplicative" part while exponent > 1: if exponent % 2 == 0: - result_s = cls._multiply_python(result_s, result_s) + result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) exponent //= 2 else: - result_m = cls._multiply_python(result_m, result_s) + result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) exponent -= 1 - result = cls._multiply_python(result_m, result_s) + result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) return result + @staticmethod + @numba.extending.register_jitable(inline="always") + def _divide_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + if b == 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") -############################################################################### -# Compile functions -############################################################################### - -CHARACTERISTIC = None # The prime characteristic `p` of the Galois field -DEGREE = None # The prime power `m` of the Galois field -IRREDUCIBLE_POLY = None # The field's primitive polynomial in integer form - -MULTIPLY = lambda a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: a * b -RECIPROCAL = lambda a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: 1 / a - -# pylint: disable=redefined-outer-name,unused-argument - - -def compile_jit(name, reset=True): - """ - Compile a JIT arithmetic function. These can be cached. - """ - if name not in compile_jit.cache: - global MULTIPLY, RECIPROCAL - - if name in ["reciprocal", "divide", "power", "log"]: - MULTIPLY = compile_jit("multiply", reset=False) - if name in ["divide", "power"]: - RECIPROCAL = compile_jit("reciprocal", reset=False) - - function = eval(f"{name}") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64)"], nopython=True, cache=True)(function) - else: - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64, int64)"], nopython=True, cache=True)(function) - - if reset: - reset_globals() - - return compile_jit.cache[name] - -compile_jit.cache = {} - - -def compile_ufunc(name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_): - """ - Compile an arithmetic ufunc. These cannot be cached as the field parameters are compiled into the binary. - """ - key = (name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_) - if key not in compile_ufunc.cache: - global CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY, MULTIPLY, RECIPROCAL - CHARACTERISTIC = CHARACTERISTIC_ - DEGREE = DEGREE_ - IRREDUCIBLE_POLY = IRREDUCIBLE_POLY_ - - if name in ["reciprocal", "divide", "power", "log"]: - MULTIPLY = compile_jit("multiply", reset=False) - if name in ["divide", "power"]: - RECIPROCAL = compile_jit("reciprocal", reset=False) - - function = eval(f"{name}_ufunc") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_ufunc.cache[key] = numba.vectorize(["int64(int64)"], nopython=True)(function) - else: - compile_ufunc.cache[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(function) - - reset_globals() - - return compile_ufunc.cache[key] - -compile_ufunc.cache = {} - - -############################################################################### -# Arithmetic explicitly calculated -############################################################################### - -def add(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a ^ b - - -def negative(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a - - -def subtract(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return a ^ b - - -@numba.extending.register_jitable(inline="always") -def multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - a in GF(2^m), can be represented as a degree m-1 polynomial in GF(2)[x] - b in GF(2^m), can be represented as a degree m-1 polynomial in GF(2)[x] - p(x) in GF(2)[x] with degree m is the primitive polynomial of GF(2^m) - - a * b = c - = (a(x) * b(x)) % p(x), in GF(2) - = c(x) - = c - """ - ORDER = CHARACTERISTIC**DEGREE - - # Re-order operands such that a > b so the while loop has less loops - if b > a: - a, b = b, a - - c = 0 - while b > 0: - if b & 0b1: - c ^= a # Add a(x) to c(x) - - b >>= 1 # Divide b(x) by x - a <<= 1 # Multiply a(x) by x - if a >= ORDER: - a ^= IRREDUCIBLE_POLY # Compute a(x) % p(x) - - return c - - -def multiply_ufunc(a, b): # pragma: no cover - return multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - From Fermat's Little Theorem: - a^(p^m - 1) = 1 (mod p^m), for a in GF(p^m) - - a * a^-1 = 1 - a * a^-1 = a^(p^m - 1) - a^-1 = a^(p^m - 2) - """ - if a == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - ORDER = CHARACTERISTIC**DEGREE - exponent = ORDER - 2 - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while exponent > 1: - if exponent % 2 == 0: - result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - exponent //= 2 - else: - result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - exponent -= 1 - - result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - return result - - -def reciprocal_ufunc(a): # pragma: no cover - return reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if a == 0: - return 0 - else: - b_inv = RECIPROCAL(b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - return MULTIPLY(a, b_inv, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -def divide_ufunc(a, b): # pragma: no cover - return divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - Square and Multiply Algorithm - - a^13 = (1) * (a)^13 - = (a) * (a)^12 - = (a) * (a^2)^6 - = (a) * (a^4)^3 - = (a * a^4) * (a^4)^2 - = (a * a^4) * (a^8) - = result_m * result_s - """ - if a == 0 and b < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b == 0: - return 1 - elif b < 0: - a = RECIPROCAL(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b = abs(b) - - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while b > 1: - if b % 2 == 0: - result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b //= 2 + if a == 0: + return 0 else: - result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b -= 1 - - result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - return result - - -def power_ufunc(a, b): # pragma: no cover - return power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - TODO: Replace this with more efficient algorithm - - a in GF(p^m) - b in GF(p^m) and generates field + b_inv = RECIPROCAL(b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + return MULTIPLY(a, b_inv, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _power_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Square and Multiply Algorithm + + a^13 = (1) * (a)^13 + = (a) * (a)^12 + = (a) * (a^2)^6 + = (a) * (a^4)^3 + = (a * a^4) * (a^4)^2 + = (a * a^4) * (a^8) + = result_m * result_s + """ + if a == 0 and b < 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - c = log(a, b), such that b^a = c - """ - if a == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") + if b == 0: + return 1 + elif b < 0: + a = RECIPROCAL(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b = abs(b) - # Naive algorithm - ORDER = CHARACTERISTIC**DEGREE - result = 1 - for i in range(0, ORDER - 1): - if result == a: - break - result = MULTIPLY(result, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + result_s = a # The "squaring" part + result_m = 1 # The "multiplicative" part - return i + while b > 1: + if b % 2 == 0: + result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b //= 2 + else: + result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b -= 1 + result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) -def log_ufunc(a, b): # pragma: no cover - return log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + return result + @staticmethod + @numba.extending.register_jitable(inline="always") + def _log_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + TODO: Replace this with more efficient algorithm + a = α^m + b is a primitive element of the field + + c = log(a, b) + a = b^c + """ + if a == 0: + raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") -def reset_globals(): - """ - Reset the global variable so when the pure-python ufuncs call these routines, they reference - the correct pure-python functions (not JIT functions or JIT-compiled ufuncs). - """ - global MULTIPLY, RECIPROCAL - MULTIPLY = multiply - RECIPROCAL = reciprocal + # Naive algorithm + ORDER = CHARACTERISTIC**DEGREE + result = 1 + for i in range(0, ORDER - 1): + if result == a: + break + result = MULTIPLY(result, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) -reset_globals() + return i diff --git a/galois/_fields/_gfp.py b/galois/_fields/_gfp.py index 5b1f1a95db..a42fc87ae8 100644 --- a/galois/_fields/_gfp.py +++ b/galois/_fields/_gfp.py @@ -7,6 +7,8 @@ from ._class import FieldClass, DirMeta from ._dtypes import DTYPES +RECIPROCAL = lambda a, *args: 1 / a + class GFpMeta(FieldClass, DirMeta): """ @@ -14,6 +16,9 @@ class GFpMeta(FieldClass, DirMeta): """ # pylint: disable=no-value-for-parameter + # Need to have a unique cache of "calculate" function for GF(p) + _FUNC_CACHE_CALCULATE = {} + def __init__(cls, name, bases, namespace, **kwargs): super().__init__(name, bases, namespace, **kwargs) cls._prime_subfield = cls @@ -22,65 +27,76 @@ def __init__(cls, name, bases, namespace, **kwargs): @property def dtypes(cls): + """ + The only valid dtypes are ones that can hold x*x for x in [0, order). + """ max_dtype = DTYPES[-1] d = [dtype for dtype in DTYPES if np.iinfo(dtype).max >= cls.order - 1 and np.iinfo(max_dtype).max >= (cls.order - 1)**2] if len(d) == 0: d = [np.object_] return d - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _calculate_jit(cls, name): - return compile_jit(name) - - def _python_func(cls, name): - return eval(f"{name}") - - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### - def _ufunc(cls, name): # Some explicit calculation functions are faster than using lookup tables. See https://github.com/mhostetter/galois/pull/92#issuecomment-835548405. if name not in cls._ufuncs and cls.ufunc_mode == "jit-lookup" and name in ["add", "negative", "subtract"]: - cls._ufuncs[name] = cls._calculate_ufunc(name) + cls._ufuncs[name] = cls._ufunc_calculate(name) return super()._ufunc(name) - def _calculate_ufunc(cls, name): - return compile_ufunc(name, cls.characteristic, cls.degree, cls._irreducible_poly_int) + def _set_globals(cls, name): + global RECIPROCAL + if name in ["divide", "power"]: + RECIPROCAL = cls._func_calculate("reciprocal", reset=False) + + def _reset_globals(cls): + global RECIPROCAL + RECIPROCAL = cls._reciprocal_calculate ############################################################################### - # Pure python arithmetic methods + # Arithmetic functions using explicit calculation ############################################################################### - def _add_python(cls, a, b): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _add_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): c = a + b - if c >= cls.order: - c -= cls.order + if c >= CHARACTERISTIC: + c -= CHARACTERISTIC return c - def _negative_python(cls, a): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _negative_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): if a == 0: return 0 else: - return cls.order - a + return CHARACTERISTIC - a - def _subtract_python(cls, a, b): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _subtract_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): if a >= b: return a - b else: - return cls.order + a - b - - def _multiply_python(cls, a, b): - return (a * b) % cls.order - - def _reciprocal_python(cls, a): + return CHARACTERISTIC + a - b + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _multiply_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + return (a * b) % CHARACTERISTIC + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _reciprocal_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + s*x + t*y = gcd(x, y) = 1 + x = p + y = a in GF(p) + t = a**-1 in GF(p) + """ if a == 0: raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - r2, r1 = cls.order, a + r2, r1 = CHARACTERISTIC, a t2, t1 = 0, 1 while r1 != 0: @@ -89,245 +105,80 @@ def _reciprocal_python(cls, a): t2, t1 = t1, t2 - q*t1 if t2 < 0: - t2 += cls.order + t2 += CHARACTERISTIC return t2 + @staticmethod + @numba.extending.register_jitable(inline="always") + def _divide_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + if b == 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") -############################################################################### -# Compile functions -############################################################################### - -CHARACTERISTIC = None # The prime characteristic `p` of the Galois field -DEGREE = None # The prime power `m` of the Galois field -IRREDUCIBLE_POLY = None # The field's primitive polynomial in integer form - -RECIPROCAL = lambda a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: 1 / a - -# pylint: disable=redefined-outer-name,unused-argument - - -def compile_jit(name, reset=True): - """ - Compile a JIT arithmetic function. These can be cached. - """ - if name not in compile_jit.cache: - global RECIPROCAL - - if name in ["divide", "power"]: - RECIPROCAL = compile_jit("reciprocal", reset=False) - - function = eval(f"{name}") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64)"], nopython=True, cache=True)(function) - else: - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64, int64)"], nopython=True, cache=True)(function) - - if reset: - reset_globals() - - return compile_jit.cache[name] - -compile_jit.cache = {} - - -def compile_ufunc(name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_): - """ - Compile an arithmetic ufunc. These cannot be cached as the field parameters are compiled into the binary. - """ - key = (name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_) - if key not in compile_ufunc.cache: - global CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY, RECIPROCAL - CHARACTERISTIC = CHARACTERISTIC_ - DEGREE = DEGREE_ - IRREDUCIBLE_POLY = IRREDUCIBLE_POLY_ - - if name in ["divide", "power"]: - RECIPROCAL = compile_jit("reciprocal", reset=False) - - function = eval(f"{name}_ufunc") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_ufunc.cache[key] = numba.vectorize(["int64(int64)"], nopython=True)(function) - else: - compile_ufunc.cache[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(function) - - reset_globals() - - return compile_ufunc.cache[key] - -compile_ufunc.cache = {} - - -############################################################################### -# Arithmetic explicitly calculated -############################################################################### - -@numba.extending.register_jitable(inline="always") -def add(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - c = a + b - if c >= CHARACTERISTIC: - c -= CHARACTERISTIC - return c - - -def add_ufunc(a, b): # pragma: no cover - return add(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def negative(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if a == 0: - return 0 - else: - return CHARACTERISTIC - a - - -def negative_ufunc(a): # pragma: no cover - return negative(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def subtract(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if a >= b: - return a - b - else: - return CHARACTERISTIC + a - b - - -def subtract_ufunc(a, b): # pragma: no cover - return subtract(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - return (a * b) % CHARACTERISTIC - - -def multiply_ufunc(a, b): # pragma: no cover - return multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - s*x + t*y = gcd(x, y) = 1 - x = p - y = a in GF(p) - t = a**-1 in GF(p) - """ - if a == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - r2, r1 = CHARACTERISTIC, a - t2, t1 = 0, 1 - - while r1 != 0: - q = r2 // r1 - r2, r1 = r1, r2 - q*r1 - t2, t1 = t1, t2 - q*t1 - - if t2 < 0: - t2 += CHARACTERISTIC - - return t2 - - -def reciprocal_ufunc(a): # pragma: no cover - return reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if a == 0: - return 0 - else: - b_inv = RECIPROCAL(b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - return (a * b_inv) % CHARACTERISTIC - - -def divide_ufunc(a, b): # pragma: no cover - return divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - Square and Multiply Algorithm - - a^13 = (1) * (a)^13 - = (a) * (a)^12 - = (a) * (a^2)^6 - = (a) * (a^4)^3 - = (a * a^4) * (a^4)^2 - = (a * a^4) * (a^8) - = result_m * result_s - """ - if a == 0 and b < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b == 0: - return 1 - elif b < 0: - a = RECIPROCAL(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b = abs(b) - - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while b > 1: - if b % 2 == 0: - result_s = (result_s * result_s) % CHARACTERISTIC - b //= 2 + if a == 0: + return 0 else: - result_m = (result_m * result_s) % CHARACTERISTIC - b -= 1 - - result = (result_m * result_s) % CHARACTERISTIC - - return result - - -def power_ufunc(a, b): # pragma: no cover - return power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - TODO: Replace this with more efficient algorithm - - a in GF(p^m) - b in GF(p^m) and generates field - - c = log(a, b), such that b^a = c - """ - if a == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") - - # Naive algorithm - ORDER = CHARACTERISTIC**DEGREE - result = 1 - for i in range(0, ORDER - 1): - if result == a: - break - result = (result * b) % CHARACTERISTIC - - return i - - -def log_ufunc(a, b): # pragma: no cover - return log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b_inv = RECIPROCAL(b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + return (a * b_inv) % CHARACTERISTIC + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _power_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Square and Multiply Algorithm + + a^13 = (1) * (a)^13 + = (a) * (a)^12 + = (a) * (a^2)^6 + = (a) * (a^4)^3 + = (a * a^4) * (a^4)^2 + = (a * a^4) * (a^8) + = result_m * result_s + """ + if a == 0 and b < 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") + if b == 0: + return 1 + elif b < 0: + a = RECIPROCAL(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b = abs(b) + + result_s = a # The "squaring" part + result_m = 1 # The "multiplicative" part + + while b > 1: + if b % 2 == 0: + result_s = (result_s * result_s) % CHARACTERISTIC + b //= 2 + else: + result_m = (result_m * result_s) % CHARACTERISTIC + b -= 1 + + result = (result_m * result_s) % CHARACTERISTIC + + return result + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _log_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + TODO: Replace this with more efficient algorithm + a = α^m + b is a primitive element of the field + + c = log(a, b) + a = b^c + """ + if a == 0: + raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") -def reset_globals(): - """ - Reset the global variable so when the pure-python ufuncs call these routines, they reference - the correct pure-python functions (not JIT functions or JIT-compiled ufuncs). - """ - global RECIPROCAL - RECIPROCAL = reciprocal + # Naive algorithm + ORDER = CHARACTERISTIC**DEGREE + result = 1 + for i in range(0, ORDER - 1): + if result == a: + break + result = (result * b) % CHARACTERISTIC -reset_globals() + return i diff --git a/galois/_fields/_gfpm.py b/galois/_fields/_gfpm.py index 0118f8858f..e63fbcde7d 100644 --- a/galois/_fields/_gfpm.py +++ b/galois/_fields/_gfpm.py @@ -7,6 +7,12 @@ from ._class import FieldClass, DirMeta from ._dtypes import DTYPES +DTYPE = np.int64 +INT_TO_POLY = lambda a, CHARACTERISTIC, DEGREE: [0,]*DEGREE +POLY_TO_INT = lambda a_vec, CHARACTERISTIC, DEGREE: 0 +MULTIPLY = lambda a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: a * b +RECIPROCAL = lambda a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: 1 / a + class GFpmMeta(FieldClass, DirMeta): """ @@ -14,6 +20,9 @@ class GFpmMeta(FieldClass, DirMeta): """ # pylint: disable=no-value-for-parameter + # Need to have a unique cache of "calculate" function for GF(p^m) + _FUNC_CACHE_CALCULATE = {} + def __init__(cls, name, bases, namespace, **kwargs): super().__init__(name, bases, namespace, **kwargs) cls._irreducible_poly_coeffs = np.array(cls._irreducible_poly.coeffs, dtype=cls.dtypes[-1]) @@ -27,28 +36,47 @@ def __init__(cls, name, bases, namespace, **kwargs): @property def dtypes(cls): + """ + The only valid dtypes are ones that can hold x*x for x in [0, order). + """ + # TODO: Is this correct? max_dtype = DTYPES[-1] d = [dtype for dtype in DTYPES if np.iinfo(dtype).max >= cls.order - 1 and np.iinfo(max_dtype).max >= (cls.order - 1)**2] if len(d) == 0: d = [np.object_] return d - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _calculate_jit(cls, name): - return compile_jit(name) - - def _python_func(cls, name): - return eval(f"{name}") - - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### + def _set_globals(cls, name): + global DTYPE, INT_TO_POLY, POLY_TO_INT, MULTIPLY, RECIPROCAL + DTYPE = np.int64 + INT_TO_POLY = cls._func_calculate("int_to_poly", reset=False) + POLY_TO_INT = cls._func_calculate("poly_to_int", reset=False) + if name in ["reciprocal", "divide", "power", "log"]: + MULTIPLY = cls._func_calculate("multiply", reset=False) + if name in ["divide", "power"]: + RECIPROCAL = cls._func_calculate("reciprocal", reset=False) + + def _reset_globals(cls): + global DTYPE, INT_TO_POLY, POLY_TO_INT, MULTIPLY, RECIPROCAL + DTYPE = np.object_ + INT_TO_POLY = cls._int_to_poly + POLY_TO_INT = cls._poly_to_int + MULTIPLY = cls._multiply_calculate + RECIPROCAL = cls._reciprocal_calculate + + def _func_calculate(cls, name, reset=True): + key = (name,) + + if key not in cls._FUNC_CACHE_CALCULATE: + # Generate extra JIT functions specific to the GF(p^m) field + if name == "int_to_poly": + cls._FUNC_CACHE_CALCULATE[key] = numba.jit(["int64[:](int64, int64, int64)"], nopython=True, cache=True)(cls._int_to_poly) + elif name == "poly_to_int": + cls._FUNC_CACHE_CALCULATE[key] = numba.jit(["int64(int64[:], int64, int64)"], nopython=True, cache=True)(cls._poly_to_int) + else: + super()._func_calculate(name, reset=reset) - def _calculate_ufunc(cls, name): - return compile_ufunc(name, cls.characteristic, cls.degree, cls._irreducible_poly_int) + return cls._FUNC_CACHE_CALCULATE[key] ############################################################################### # Ufunc routines @@ -117,400 +145,186 @@ def _ufunc_routine_subtract(cls, ufunc, method, inputs, kwargs, meta): return output ############################################################################### - # Pure python arithmetic methods + # Arithmetic functions using explicit calculation ############################################################################### - def _int_to_poly(cls, a): - a_vec = np.zeros(cls.degree, dtype=cls.dtypes[-1]) - for i in range(0, cls.degree): - q = a // cls.characteristic**(cls.degree - 1 - i) - a -= q*cls.characteristic**(cls.degree - 1 - i) + @staticmethod + @numba.extending.register_jitable(inline="always") + def _int_to_poly(a, CHARACTERISTIC, DEGREE): + """ + Convert the integer representation to vector/polynomial representation + """ + a_vec = np.zeros(DEGREE, dtype=DTYPE) + for i in range(0, DEGREE): + q = a // CHARACTERISTIC**(DEGREE - 1 - i) + a -= q*CHARACTERISTIC**(DEGREE - 1 - i) a_vec[i] = q return a_vec - def _poly_to_int(cls, a_vec): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _poly_to_int(a_vec, CHARACTERISTIC, DEGREE): + """ + Convert the integer representation to vector/polynomial representation + """ a = 0 - for i in range(0, cls.degree): - a += a_vec[i]*cls.characteristic**(cls.degree - 1 - i) + for i in range(0, DEGREE): + a += a_vec[i]*CHARACTERISTIC**(DEGREE - 1 - i) return a - def _add_python(cls, a, b): - a_vec = cls._int_to_poly(a) - b_vec = cls._int_to_poly(b) - c_vec = (a_vec + b_vec) % cls.characteristic - return cls._poly_to_int(c_vec) - - def _negative_python(cls, a): - a_vec = cls._int_to_poly(a) - c_vec = (-a_vec) % cls.characteristic - return cls._poly_to_int(c_vec) - - def _subtract_python(cls, a, b): - a_vec = cls._int_to_poly(a) - b_vec = cls._int_to_poly(b) - c_vec = (a_vec - b_vec) % cls.characteristic - return cls._poly_to_int(c_vec) - - def _multiply_python(cls, a, b): - a_vec = cls._int_to_poly(a) - b_vec = cls._int_to_poly(b) - - c_vec = np.zeros(cls.degree, dtype=cls.dtypes[-1]) - for _ in range(cls.degree): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _add_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) + b_vec = INT_TO_POLY(b, CHARACTERISTIC, DEGREE) + c_vec = (a_vec + b_vec) % CHARACTERISTIC + return POLY_TO_INT(c_vec, CHARACTERISTIC, DEGREE) + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _negative_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) + a_vec = (-a_vec) % CHARACTERISTIC + return POLY_TO_INT(a_vec, CHARACTERISTIC, DEGREE) + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _subtract_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) + b_vec = INT_TO_POLY(b, CHARACTERISTIC, DEGREE) + c_vec = (a_vec - b_vec) % CHARACTERISTIC + return POLY_TO_INT(c_vec, CHARACTERISTIC, DEGREE) + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _multiply_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) + b_vec = INT_TO_POLY(b, CHARACTERISTIC, DEGREE) + + # The irreducible polynomial with the x^degree term removed + irreducible_poly_vec = INT_TO_POLY(IRREDUCIBLE_POLY - CHARACTERISTIC**DEGREE, CHARACTERISTIC, DEGREE) + + c_vec = np.zeros(DEGREE, dtype=DTYPE) + for _ in range(DEGREE): if b_vec[-1] > 0: - c_vec = (c_vec + b_vec[-1]*a_vec) % cls.characteristic + c_vec = (c_vec + b_vec[-1]*a_vec) % CHARACTERISTIC # Multiply a(x) by x - q = a_vec[0] # Don't need to divide by the leading coefficient of p(x) because it must be 1 + q = a_vec[0] a_vec[:-1] = a_vec[1:] a_vec[-1] = 0 - # Reduce a(x) modulo the irreducible polynomial p(x) + # Reduce a(x) modulo the irreducible polynomial if q > 0: - a_vec = (a_vec - q*cls._irreducible_poly_coeffs[1:]) % cls.characteristic + a_vec = (a_vec - q*irreducible_poly_vec) % CHARACTERISTIC # Divide b(x) by x b_vec[1:] = b_vec[:-1] b_vec[0] = 0 - return cls._poly_to_int(c_vec) + return POLY_TO_INT(c_vec, CHARACTERISTIC, DEGREE) - def _reciprocal_python(cls, a): + @staticmethod + @numba.extending.register_jitable(inline="always") + def _reciprocal_calculate(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + From Fermat's Little Theorem: + a^(p^m - 1) = 1 (mod p^m), for a in GF(p^m) + + a * a^-1 = 1 + a * a^-1 = a^(p^m - 1) + a^-1 = a^(p^m - 2) + """ if a == 0: raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - exponent = cls.order - 2 + ORDER = CHARACTERISTIC**DEGREE + exponent = ORDER - 2 result_s = a # The "squaring" part result_m = 1 # The "multiplicative" part while exponent > 1: if exponent % 2 == 0: - result_s = cls._multiply_python(result_s, result_s) + result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) exponent //= 2 else: - result_m = cls._multiply_python(result_m, result_s) + result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) exponent -= 1 - result = cls._multiply_python(result_m, result_s) + result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) return result + @staticmethod + @numba.extending.register_jitable(inline="always") + def _divide_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + if b == 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") -############################################################################### -# Compile functions -############################################################################### - -CHARACTERISTIC = None # The prime characteristic `p` of the Galois field -DEGREE = None # The prime power `m` of the Galois field -IRREDUCIBLE_POLY = None # The field's primitive polynomial in integer form - -DTYPE = np.int64 -INT_TO_POLY = lambda a, CHARACTERISTIC, DEGREE: [0,]*DEGREE -POLY_TO_INT = lambda a_vec, CHARACTERISTIC, DEGREE: 0 -MULTIPLY = lambda a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: a * b -RECIPROCAL = lambda a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY: 1 / a - -# pylint: disable=redefined-outer-name,unused-argument - - -def compile_jit(name, reset=True): - """ - Compile a JIT arithmetic function. These can be cached. - """ - if name not in compile_jit.cache: - function = eval(f"{name}") - - if name == "int_to_poly": - compile_jit.cache[name] = numba.jit(["int64[:](int64, int64, int64)"], nopython=True, cache=True)(function) - elif name == "poly_to_int": - compile_jit.cache[name] = numba.jit(["int64(int64[:], int64, int64)"], nopython=True, cache=True)(function) + if a == 0: + return 0 else: - global DTYPE, INT_TO_POLY, POLY_TO_INT, MULTIPLY, RECIPROCAL + b_inv = RECIPROCAL(b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + return MULTIPLY(a, b_inv, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _power_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + Square and Multiply Algorithm + + a^13 = (1) * (a)^13 + = (a) * (a)^12 + = (a) * (a^2)^6 + = (a) * (a^4)^3 + = (a * a^4) * (a^4)^2 + = (a * a^4) * (a^8) + = result_m * result_s + """ + if a == 0 and b < 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - DTYPE = np.int64 - INT_TO_POLY = compile_jit("int_to_poly", reset=False) - POLY_TO_INT = compile_jit("poly_to_int", reset=False) + if b == 0: + return 1 + elif b < 0: + a = RECIPROCAL(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b = abs(b) - if name in ["reciprocal", "divide", "power", "log"]: - MULTIPLY = compile_jit("multiply", reset=False) - if name in ["divide", "power"]: - RECIPROCAL = compile_jit("reciprocal", reset=False) + result_s = a # The "squaring" part + result_m = 1 # The "multiplicative" part - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64)"], nopython=True, cache=True)(function) + while b > 1: + if b % 2 == 0: + result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b //= 2 else: - compile_jit.cache[name] = numba.jit(["int64(int64, int64, int64, int64, int64)"], nopython=True, cache=True)(function) - - if reset: - reset_globals() - - return compile_jit.cache[name] - -compile_jit.cache = {} - - -def compile_ufunc(name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_): - """ - Compile an arithmetic ufunc. These cannot be cached as the field parameters are compiled into the binary. - """ - key = (name, CHARACTERISTIC_, DEGREE_, IRREDUCIBLE_POLY_) - if key not in compile_ufunc.cache: - global CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY, DTYPE, INT_TO_POLY, POLY_TO_INT, MULTIPLY, RECIPROCAL - CHARACTERISTIC = CHARACTERISTIC_ - DEGREE = DEGREE_ - IRREDUCIBLE_POLY = IRREDUCIBLE_POLY_ - - DTYPE = np.int64 - INT_TO_POLY = compile_jit("int_to_poly", reset=False) - POLY_TO_INT = compile_jit("poly_to_int", reset=False) - - if name in ["reciprocal", "divide", "power", "log"]: - MULTIPLY = compile_jit("multiply", reset=False) - if name in ["divide", "power"]: - RECIPROCAL = compile_jit("reciprocal", reset=False) - - function = eval(f"{name}_ufunc") - if FieldClass._UFUNC_TYPE[name] == "unary": - compile_ufunc.cache[key] = numba.vectorize(["int64(int64)"], nopython=True)(function) - else: - compile_ufunc.cache[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(function) - - reset_globals() - - return compile_ufunc.cache[key] - -compile_ufunc.cache = {} + result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + b -= 1 + result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) -############################################################################### -# Helper functions -############################################################################### - -def int_to_poly(a, CHARACTERISTIC, DEGREE): - """ - Convert the integer representation to vector/polynomial representation - """ - a_vec = np.zeros(DEGREE, dtype=DTYPE) - for i in range(0, DEGREE): - q = a // CHARACTERISTIC**(DEGREE - 1 - i) - a -= q*CHARACTERISTIC**(DEGREE - 1 - i) - a_vec[i] = q - return a_vec - - -def poly_to_int(a_vec, CHARACTERISTIC, DEGREE): - """ - Convert the integer representation to vector/polynomial representation - """ - a = 0 - for i in range(0, DEGREE): - a += a_vec[i]*CHARACTERISTIC**(DEGREE - 1 - i) - return a - - -############################################################################### -# Arithmetic explicitly calculated -############################################################################### - -@numba.extending.register_jitable(inline="always") -def add(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) - b_vec = INT_TO_POLY(b, CHARACTERISTIC, DEGREE) - c_vec = (a_vec + b_vec) % CHARACTERISTIC - return POLY_TO_INT(c_vec, CHARACTERISTIC, DEGREE) - - -def add_ufunc(a, b): # pragma: no cover - return add(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def negative(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) - a_vec = (-a_vec) % CHARACTERISTIC - return POLY_TO_INT(a_vec, CHARACTERISTIC, DEGREE) - - -def negative_ufunc(a): # pragma: no cover - return negative(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def subtract(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) - b_vec = INT_TO_POLY(b, CHARACTERISTIC, DEGREE) - c_vec = (a_vec - b_vec) % CHARACTERISTIC - return POLY_TO_INT(c_vec, CHARACTERISTIC, DEGREE) - - -def subtract_ufunc(a, b): # pragma: no cover - return subtract(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - a_vec = INT_TO_POLY(a, CHARACTERISTIC, DEGREE) - b_vec = INT_TO_POLY(b, CHARACTERISTIC, DEGREE) - - # The irreducible polynomial with the x^degree term removed - irreducible_poly_vec = INT_TO_POLY(IRREDUCIBLE_POLY - CHARACTERISTIC**DEGREE, CHARACTERISTIC, DEGREE) - - c_vec = np.zeros(DEGREE, dtype=DTYPE) - for _ in range(DEGREE): - if b_vec[-1] > 0: - c_vec = (c_vec + b_vec[-1]*a_vec) % CHARACTERISTIC - - # Multiply a(x) by x - q = a_vec[0] - a_vec[:-1] = a_vec[1:] - a_vec[-1] = 0 - - # Reduce a(x) modulo the irreducible polynomial - if q > 0: - a_vec = (a_vec - q*irreducible_poly_vec) % CHARACTERISTIC - - # Divide b(x) by x - b_vec[1:] = b_vec[:-1] - b_vec[0] = 0 - - return POLY_TO_INT(c_vec, CHARACTERISTIC, DEGREE) - - -def multiply_ufunc(a, b): # pragma: no cover - return multiply(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - From Fermat's Little Theorem: - a^(p^m - 1) = 1 (mod p^m), for a in GF(p^m) - - a * a^-1 = 1 - a * a^-1 = a^(p^m - 1) - a^-1 = a^(p^m - 2) - """ - if a == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - ORDER = CHARACTERISTIC**DEGREE - exponent = ORDER - 2 - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while exponent > 1: - if exponent % 2 == 0: - result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - exponent //= 2 - else: - result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - exponent -= 1 - - result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - return result - - -def reciprocal_ufunc(a): # pragma: no cover - return reciprocal(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if a == 0: - return 0 - else: - b_inv = RECIPROCAL(b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - return MULTIPLY(a, b_inv, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -def divide_ufunc(a, b): # pragma: no cover - return divide(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - Square and Multiply Algorithm - - a^13 = (1) * (a)^13 - = (a) * (a)^12 - = (a) * (a^2)^6 - = (a) * (a^4)^3 - = (a * a^4) * (a^4)^2 - = (a * a^4) * (a^8) - = result_m * result_s - """ - if a == 0 and b < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b == 0: - return 1 - elif b < 0: - a = RECIPROCAL(a, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b = abs(b) - - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while b > 1: - if b % 2 == 0: - result_s = MULTIPLY(result_s, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b //= 2 - else: - result_m = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - b -= 1 - - result = MULTIPLY(result_m, result_s, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - return result - - -def power_ufunc(a, b): # pragma: no cover - return power(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - -@numba.extending.register_jitable(inline="always") -def log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): # pragma: no cover - """ - TODO: Replace this with more efficient algorithm - - a in GF(p^m) - b in GF(p^m) and generates field - - c = log(a, b), such that b^a = c - """ - if a == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") - - # Naive algorithm - ORDER = CHARACTERISTIC**DEGREE - result = 1 - for i in range(0, ORDER - 1): - if result == a: - break - result = MULTIPLY(result, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - - return i - + return result -def log_ufunc(a, b): # pragma: no cover - return log(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) + @staticmethod + @numba.extending.register_jitable(inline="always") + def _log_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY): + """ + TODO: Replace this with more efficient algorithm + a = α^m + b is a primitive element of the field + + c = log(a, b) + a = b^c + """ + if a == 0: + raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") + # Naive algorithm + ORDER = CHARACTERISTIC**DEGREE + result = 1 + for i in range(0, ORDER - 1): + if result == a: + break + result = MULTIPLY(result, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) -def reset_globals(): - """ - Reset the global variable so when the pure-python ufuncs call these routines, they reference - the correct pure-python functions (not JIT functions or JIT-compiled ufuncs). - """ - global DTYPE, INT_TO_POLY, POLY_TO_INT, MULTIPLY, RECIPROCAL - DTYPE = np.object_ - INT_TO_POLY = int_to_poly - POLY_TO_INT = poly_to_int - MULTIPLY = multiply - RECIPROCAL = reciprocal - -reset_globals() + return i diff --git a/galois/_fields/_lookup.py b/galois/_fields/_lookup.py index e446601040..ecee7f4578 100644 --- a/galois/_fields/_lookup.py +++ b/galois/_fields/_lookup.py @@ -5,29 +5,21 @@ from numba import int64 import numpy as np -from ._properties import PropertiesMeta +from ._calculate import CalculateMeta -class LookupMeta(PropertiesMeta): +class LookupMeta(CalculateMeta): """ A mixin class that provides Galois field arithmetic using lookup tables. """ - # pylint: disable=no-value-for-parameter + # pylint: disable=no-value-for-parameter,abstract-method,unused-argument + # Function signatures for JIT-compiled "lookup" arithmetic functions _UNARY_LOOKUP_SIG = numba.types.FunctionType(int64(int64, int64[:], int64[:], int64[:], int64)) _BINARY_LOOKUP_SIG = numba.types.FunctionType(int64(int64, int64, int64[:], int64[:], int64[:], int64)) - # Lookup table of ufuncs to unary/binary type needed for LookupMeta, CalculateMeta, etc - _UFUNC_TYPE = { - "add": "binary", - "negative": "unary", - "subtract": "binary", - "multiply": "binary", - "reciprocal": "unary", - "divide": "binary", - "power": "binary", - "log": "binary", - } + _FUNC_CACHE_LOOKUP = {} + _UFUNC_CACHE_LOOKUP = {} def __init__(cls, name, bases, namespace, **kwargs): super().__init__(name, bases, namespace, **kwargs) @@ -36,17 +28,18 @@ def __init__(cls, name, bases, namespace, **kwargs): cls._ZECH_LOG = np.array([], dtype=np.int64) cls._ZECH_E = 0 - def _compile_ufuncs(cls): - cls._ufuncs = {} # Reset the dictionary so each ufunc will get recompiled - if cls.ufunc_mode == "jit-lookup": - cls._build_lookup_tables() - def _build_lookup_tables(cls): + """ + Construct EXP, LOG, and ZECH_LOG lookup tables to be used in the "lookup" arithmetic functions + """ + # Only construct the LUTs for this field once if cls._EXP.size > 0: return order = cls.order primitive_element = int(cls.primitive_element) + add = lambda a, b: cls._func_python("add")(a, b, cls.characteristic, cls.degree, cls._irreducible_poly_int) + multiply = lambda a, b: cls._func_python("multiply")(a, b, cls.characteristic, cls.degree, cls._irreducible_poly_int) cls._EXP = np.zeros(2*order, dtype=np.int64) cls._LOG = np.zeros(order, dtype=np.int64) @@ -61,7 +54,7 @@ def _build_lookup_tables(cls): cls._LOG[0] = 0 # Technically -Inf for i in range(1, order): # Increment by multiplying by the primitive element, which is a multiplicative generator of the field - element = cls._multiply_python(element, primitive_element) + element = multiply(element, primitive_element) cls._EXP[i] = element # Assign to the log lookup table but skip indices greater than or equal to `order - 1` @@ -71,7 +64,7 @@ def _build_lookup_tables(cls): # Compute Zech log lookup table for i in range(0, order): - one_plus_element = cls._add_python(1, cls._EXP[i]) + one_plus_element = add(1, cls._EXP[i]) cls._ZECH_LOG[i] = cls._LOG[one_plus_element] if not cls._EXP[order - 1] == 1: @@ -84,300 +77,238 @@ def _build_lookup_tables(cls): # Double the EXP table to prevent computing a `% (order - 1)` on every multiplication lookup cls._EXP[order:2*order] = cls._EXP[1:1 + order] - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _lookup_jit(cls, name): # pylint: disable=no-self-use - return compile_jit(name) + def _func_lookup(cls, name): # pylint: disable=no-self-use + """ + Returns an arithmetic function using lookup tables. These functions are once-compiled and shared for all Galois fields. The only difference + between Galois fields are the lookup tables that are passed in as inputs. + """ + key = (name,) + + if key not in cls._FUNC_CACHE_LOOKUP: + function = getattr(cls, f"_{name}_lookup") + if cls._UFUNC_TYPE[name] == "unary": + cls._FUNC_CACHE_LOOKUP[key] = numba.jit("int64(int64, int64[:], int64[:], int64[:], int64)", nopython=True, cache=True)(function) + else: + cls._FUNC_CACHE_LOOKUP[key] = numba.jit("int64(int64, int64, int64[:], int64[:], int64[:], int64)", nopython=True, cache=True)(function) + + return cls._FUNC_CACHE_LOOKUP[key] + + def _ufunc_lookup(cls, name): + """ + Returns an arithmetic ufunc using lookup tables. These ufuncs are compiled for each Galois field since the lookup tables are compiled + into the ufuncs as constants. + """ + key = (name, cls.characteristic, cls.degree, cls._irreducible_poly_int) + + if key not in cls._UFUNC_CACHE_LOOKUP: + EXP = cls._EXP + LOG = cls._LOG + ZECH_LOG = cls._ZECH_LOG + ZECH_E = cls._ZECH_E + + function = getattr(cls, f"_{name}_lookup") + if cls._UFUNC_TYPE[name] == "unary": + cls._UFUNC_CACHE_LOOKUP[key] = numba.vectorize(["int64(int64)"], nopython=True)(lambda a: function(a, EXP, LOG, ZECH_LOG, ZECH_E)) + else: + cls._UFUNC_CACHE_LOOKUP[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(lambda a, b: function(a, b, EXP, LOG, ZECH_LOG, ZECH_E)) + + return cls._UFUNC_CACHE_LOOKUP[key] ############################################################################### - # Individual ufuncs, compiled on-demand + # Arithmetic functions using lookup tables ############################################################################### - def _lookup_ufunc(cls, name): - return compile_ufunc(name, cls.characteristic, cls.degree, cls._irreducible_poly_int, cls._EXP, cls._LOG, cls._ZECH_LOG, cls._ZECH_E) - - -############################################################################### -# Compile arithmetic functions using lookup tables -############################################################################### - -EXP = [] # EXP[i] = α^i -LOG = [] # LOG[i] = x, such that α^x = i -ZECH_LOG = [] # ZECH_LOG[i] = log(1 + α^i) -ZECH_E = 0 # α^ZECH_E = -1, ZECH_LOG[ZECH_E] = -Inf - -# pylint: disable=redefined-outer-name,unused-argument - - -def compile_jit(name): - """ - Compile a JIT arithmetic function. These can be cached. - """ - if name not in compile_jit.cache: - function = eval(f"{name}") - if LookupMeta._UFUNC_TYPE[name] == "unary": - compile_jit.cache[name] = numba.jit("int64(int64, int64[:], int64[:], int64[:], int64)", nopython=True, cache=True)(function) - else: - compile_jit.cache[name] = numba.jit("int64(int64, int64, int64[:], int64[:], int64[:], int64)", nopython=True, cache=True)(function) - - return compile_jit.cache[name] - -compile_jit.cache = {} - - -def compile_ufunc(name, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY, EXP_, LOG_, ZECH_LOG_, ZECH_E_): - """ - Compile an arithmetic ufunc. These cannot be cached as the lookup tables are compiled into the binary. - """ - key = (name, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY) - if key not in compile_ufunc.cache: - global EXP, LOG, ZECH_LOG, ZECH_E - EXP = EXP_ - LOG = LOG_ - ZECH_LOG = ZECH_LOG_ - ZECH_E = ZECH_E_ - - function = eval(f"{name}_ufunc") - if LookupMeta._UFUNC_TYPE[name] == "unary": - compile_ufunc.cache[key] = numba.vectorize(["int64(int64)"], nopython=True)(function) - else: - compile_ufunc.cache[key] = numba.vectorize(["int64(int64, int64)"], nopython=True)(function) - - return compile_ufunc.cache[key] - -compile_ufunc.cache = {} - - -############################################################################### -# Arithmetic using lookup tables -############################################################################### - -@numba.extending.register_jitable(inline="always") -def add(a, b, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - b in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - a + b = α^m + α^n - = α^m * (1 + α^(n - m)) # If n is larger, factor out α^m - = α^m * α^ZECH_LOG(n - m) - = α^(m + ZECH_LOG(n - m)) - """ - # LOG[0] = -Inf, so catch these conditions - if a == 0: - return b - elif b == 0: - return a - - m = LOG[a] - n = LOG[b] - - if m > n: - # We want to factor out α^m, where m is smaller than n, such that `n - m` is always positive. If - # m is larger than n, switch a and b in the addition. - m, n = n, m - - if n - m == ZECH_E: - # ZECH_LOG[ZECH_E] = -Inf and α^(-Inf) = 0 - return 0 - else: - return EXP[m + ZECH_LOG[n - m]] - - -def add_ufunc(a, b): # pragma: no cover - return add(a, b, EXP, LOG, ZECH_LOG, ZECH_E) - - -@numba.extending.register_jitable(inline="always") -def negative(a, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - -a = -α^n - = -1 * α^n - = α^e * α^n - = α^(e + n) - """ - if a == 0: # LOG[0] = -Inf, so catch this condition - return 0 - else: - n = LOG[a] - return EXP[ZECH_E + n] - - -def negative_ufunc(a): # pragma: no cover - return negative(a, EXP, LOG, ZECH_LOG, ZECH_E) - - -@numba.extending.register_jitable(inline="always") -def subtract(a, b, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - b in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - a - b = α^m - α^n - = α^m + (-α^n) - = α^m + (-1 * α^n) - = α^m + (α^e * α^n) - = α^m + α^(e + n) - """ - ORDER = LOG.size - - # Same as addition if n = LOG[b] + e - m = LOG[a] - n = LOG[b] + ZECH_E + @staticmethod + @numba.extending.register_jitable(inline="always") + def _add_lookup(a, b, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + b = α^n + + a + b = α^m + α^n + = α^m * (1 + α^(n - m)) # If n is larger, factor out α^m + = α^m * α^ZECH_LOG(n - m) + = α^(m + ZECH_LOG(n - m)) + """ + if a == 0: + return b + elif b == 0: + return a - # LOG[0] = -Inf, so catch these conditions - if b == 0: - return a - elif a == 0: - return EXP[n] - - if m > n: - # We want to factor out α^m, where m is smaller than n, such that `n - m` is always positive. If - # m is larger than n, switch a and b in the addition. - m, n = n, m - - z = n - m - if z == ZECH_E: - # ZECH_LOG[ZECH_E] = -Inf and α^(-Inf) = 0 - return 0 - if z >= ORDER - 1: - # Reduce index of ZECH_LOG by the multiplicative order of the field, i.e. `order - 1` - z -= ORDER - 1 - - return EXP[m + ZECH_LOG[z]] - - -def subtract_ufunc(a, b): # pragma: no cover - return subtract(a, b, EXP, LOG, ZECH_LOG, ZECH_E) - - -@numba.extending.register_jitable(inline="always") -def multiply(a, b, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - b in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - a * b = α^m * α^n - = α^(m + n) - """ - if a == 0 or b == 0: # LOG[0] = -Inf, so catch these conditions - return 0 - else: m = LOG[a] n = LOG[b] - return EXP[m + n] - -def multiply_ufunc(a, b): # pragma: no cover - return multiply(a, b, EXP, LOG, ZECH_LOG, ZECH_E) - - -@numba.extending.register_jitable(inline="always") -def reciprocal(a, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - 1 / a = 1 / α^m - = α^(-m) - = 1 * α^(-m) - = α^(ORDER - 1) * α^(-m) - = α^(ORDER - 1 - m) - """ - if a == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - ORDER = LOG.size - m = LOG[a] - return EXP[(ORDER - 1) - m] - - -def reciprocal_ufunc(a): # pragma: no cover - return reciprocal(a, EXP, LOG, ZECH_LOG, ZECH_E) - - -@numba.extending.register_jitable(inline="always") -def divide(a, b, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - b in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - a / b = α^m / α^n - = α^(m - n) - = 1 * α^(m - n) - = α^(ORDER - 1) * α^(m - n) - = α^(ORDER - 1 + m - n) - """ - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") + if m > n: + # We want to factor out α^m, where m is smaller than n, such that `n - m` is always positive. If m is + # larger than n, switch a and b in the addition. + m, n = n, m - if a == 0: # LOG[0] = -Inf, so catch this condition - return 0 - else: + if n - m == ZECH_E: + # zech_log(zech_e) = -Inf and α^(-Inf) = 0 + return 0 + else: + return EXP[m + ZECH_LOG[n - m]] + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _negative_lookup(a, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + + -a = -α^m + = -1 * α^m + = α^e * α^m + = α^(e + m) + """ + if a == 0: + return 0 + else: + m = LOG[a] + return EXP[ZECH_E + m] + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _subtract_lookup(a, b, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + b = α^n + + a - b = α^m - α^n + = α^m + (-α^n) + = α^m + (-1 * α^n) + = α^m + (α^e * α^n) + = α^m + α^(e + n) + """ ORDER = LOG.size - m = LOG[a] - n = LOG[b] - return EXP[(ORDER - 1) + m - n] # We add `ORDER - 1` to guarantee the index is non-negative - - -def divide_ufunc(a, b): # pragma: no cover - return divide(a, b, EXP, LOG, ZECH_LOG, ZECH_E) + # Same as addition if n = log(b) + e + m = LOG[a] + n = LOG[b] + ZECH_E + + if b == 0: + return a + elif a == 0: + return EXP[n] + + if m > n: + # We want to factor out α^m, where m is smaller than n, such that `n - m` is always positive. If m is + # larger than n, switch a and b in the addition. + m, n = n, m + + z = n - m + if z == ZECH_E: + # zech_log(zech_e) = -Inf and α^(-Inf) = 0 + return 0 + if z >= ORDER - 1: + # Reduce index of ZECH_LOG by the multiplicative order of the field `ORDER - 1` + z -= ORDER - 1 + + return EXP[m + ZECH_LOG[z]] + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _multiply_lookup(a, b, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + b = α^n + + a * b = α^m * α^n + = α^(m + n) + """ + if a == 0 or b == 0: + return 0 + else: + m = LOG[a] + n = LOG[b] + return EXP[m + n] + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _reciprocal_lookup(a, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + + 1 / a = 1 / α^m + = α^(-m) + = 1 * α^(-m) + = α^(ORDER - 1) * α^(-m) + = α^(ORDER - 1 - m) + """ + if a == 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") -@numba.extending.register_jitable(inline="always") -def power(a, b_int, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - b_int in Z - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - a ** b_int = α^m ** b_int - = α^(m * b_int) - = α^(m * ((b_int // (ORDER - 1))*(ORDER - 1) + b_int % (ORDER - 1))) - = α^(m * ((b_int // (ORDER - 1))*(ORDER - 1)) * α^(m * (b_int % (ORDER - 1))) - = 1 * α^(m * (b_int % (ORDER - 1))) - = α^(m * (b_int % (ORDER - 1))) - """ - if a == 0 and b_int < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b_int == 0: - return 1 - elif a == 0: # LOG[0] = -Inf, so catch this condition - return 0 - else: ORDER = LOG.size m = LOG[a] - return EXP[(m * b_int) % (ORDER - 1)] - - -def power_ufunc(a, b_int): # pragma: no cover - return power(a, b_int, EXP, LOG, ZECH_LOG, ZECH_E) - - -@numba.extending.register_jitable(inline="always") -def log(beta, alpha, EXP, LOG, ZECH_LOG, ZECH_E): # pragma: no cover - """ - a in GF(p^m) - α is a primitive element of GF(p^m), such that GF(p^m) = {0, 1, α^1, ..., α^(p^m - 2)} - - log(beta, α) = log(α^m, α) - = m - """ - if beta == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") - - return LOG[beta] - - -def log_ufunc(beta, alpha): # pragma: no cover - return log(beta, alpha, EXP, LOG, ZECH_LOG, ZECH_E) + return EXP[(ORDER - 1) - m] + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _divide_lookup(a, b, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + b = α^n + + a / b = α^m / α^n + = α^(m - n) + = 1 * α^(m - n) + = α^(ORDER - 1) * α^(m - n) + = α^(ORDER - 1 + m - n) + """ + if b == 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") + + if a == 0: + return 0 + else: + ORDER = LOG.size + m = LOG[a] + n = LOG[b] + return EXP[(ORDER - 1) + m - n] # We add `ORDER - 1` to guarantee the index is non-negative + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _power_lookup(a, b, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + b in Z + + a ** b = α^m ** b + = α^(m * b) + = α^(m * ((b // (ORDER - 1))*(ORDER - 1) + b % (ORDER - 1))) + = α^(m * ((b // (ORDER - 1))*(ORDER - 1)) * α^(m * (b % (ORDER - 1))) + = 1 * α^(m * (b % (ORDER - 1))) + = α^(m * (b % (ORDER - 1))) + """ + if a == 0 and b < 0: + raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") + + if b == 0: + return 1 + elif a == 0: + return 0 + else: + ORDER = LOG.size + m = LOG[a] + return EXP[(m * b) % (ORDER - 1)] # TODO: Do b % (ORDER - 1) first? b could be very large and overflow int64 + + @staticmethod + @numba.extending.register_jitable(inline="always") + def _log_lookup(beta, alpha, EXP, LOG, ZECH_LOG, ZECH_E): + """ + α is a primitive element of GF(p^m) + a = α^m + + log(beta, α) = log(α^m, α) + = m + """ + if beta == 0: + raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") + + return LOG[beta] diff --git a/galois/_fields/_properties.py b/galois/_fields/_properties.py index e78475d8c7..4b437d770e 100644 --- a/galois/_fields/_properties.py +++ b/galois/_fields/_properties.py @@ -378,10 +378,11 @@ def ufunc_modes(cls): galois.GF(31).ufunc_modes galois.GF(2**100).ufunc_modes """ - if cls.dtypes == [np.object_]: - return ["python-calculate"] - else: - return ["jit-lookup", "jit-calculate"] + # if cls.dtypes == [np.object_]: + # return ["python-calculate"] + # else: + # return ["jit-lookup", "jit-calculate"] + return ["jit-lookup", "jit-calculate", "python-calculate"] @property def default_ufunc_mode(cls): diff --git a/galois/_fields/_ufuncs.py b/galois/_fields/_ufuncs.py index 6b6b42f4c9..65b4298cdd 100644 --- a/galois/_fields/_ufuncs.py +++ b/galois/_fields/_ufuncs.py @@ -11,7 +11,7 @@ class UfuncMeta(LookupMeta, CalculateMeta): """ A mixin class that provides the basics for compiling ufuncs. """ - # pylint: disable=no-value-for-parameter + # pylint: disable=no-value-for-parameter,abstract-method _UNSUPPORTED_UFUNCS_UNARY = [ np.invert, np.sqrt, @@ -62,46 +62,27 @@ def __init__(cls, name, bases, namespace, **kwargs): cls._ufuncs = {} def _compile_ufuncs(cls): + """ + Compile/re-compile the ufuncs based on the `ufunc_mode`. This may be supplemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. + """ cls._ufuncs = {} # Reset the dictionary so each ufunc will get recompiled + if cls.ufunc_mode == "jit-lookup": cls._build_lookup_tables() - ############################################################################### - # Individual JIT arithmetic functions, pre-compiled (cached) - ############################################################################### - - def _calculate_jit(cls, name): - """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. - """ - raise NotImplementedError - - def _python_func(cls, name): + def _ufunc(cls, name): """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. + Returns the ufunc for the specific type of arithmetic. The ufunc compilation is based on `ufunc_mode`. """ - raise NotImplementedError - - ############################################################################### - # Individual ufuncs, compiled on-demand - ############################################################################### - - def _ufunc(cls, name): if name not in cls._ufuncs: if cls.ufunc_mode == "jit-lookup": - cls._ufuncs[name] = cls._lookup_ufunc(name) + cls._ufuncs[name] = cls._ufunc_lookup(name) elif cls.ufunc_mode == "jit-calculate": - cls._ufuncs[name] = cls._calculate_ufunc(name) + cls._ufuncs[name] = cls._ufunc_calculate(name) else: - cls._ufuncs[name] = cls._python_ufunc(name) + cls._ufuncs[name] = cls._ufunc_python(name) return cls._ufuncs[name] - def _calculate_ufunc(cls, name): - """ - To be implemented in GF2Meta, GF2mMeta, GFpMeta, and GFpmMeta. - """ - raise NotImplementedError - ############################################################################### # Input/output conversion functions ############################################################################### @@ -280,96 +261,3 @@ def _ufunc_routine_log(cls, ufunc, method, inputs, kwargs, meta): # pylint: dis def _ufunc_routine_matmul(cls, ufunc, method, inputs, kwargs, meta): # pylint: disable=unused-argument cls._verify_method_only_call(ufunc, method) return cls._matmul(*inputs, **kwargs) - - ############################################################################### - # Pure-python arithmetic methods using explicit calculation - ############################################################################### - - def _add_python(cls, a, b): - raise NotImplementedError - - def _negative_python(cls, a): - raise NotImplementedError - - def _subtract_python(cls, a, b): - raise NotImplementedError - - def _multiply_python(cls, a, b): - raise NotImplementedError - - def _reciprocal_python(cls, a): - raise NotImplementedError - - def _divide_python(cls, a, b): - if b == 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if a == 0: - return 0 - else: - b_inv = cls._reciprocal_python(b) - return cls._multiply_python(a, b_inv) - - def _power_python(cls, a, b): - if a == 0 and b < 0: - raise ZeroDivisionError("Cannot compute the multiplicative inverse of 0 in a Galois field.") - - if b == 0: - return 1 - elif b < 0: - a = cls._reciprocal_python(a) - b = abs(b) - - result_s = a # The "squaring" part - result_m = 1 # The "multiplicative" part - - while b > 1: - if b % 2 == 0: - result_s = cls._multiply_python(result_s, result_s) - b //= 2 - else: - result_m = cls._multiply_python(result_m, result_s) - b -= 1 - - result = cls._multiply_python(result_m, result_s) - - return result - - def _log_python(cls, a, b): - """ - TODO: Replace this with a more efficient algorithm - """ - if a == 0: - raise ArithmeticError("Cannot compute the discrete logarithm of 0 in a Galois field.") - - # Naive algorithm - result = 1 - for i in range(0, cls.order - 1): - if result == a: - break - result = cls._multiply_python(result, b) - - return i - - # N = cls.order - # alpha = cls._primitive_element_int - # n = N - 1 # Multiplicative order of the group - # x, a, b = 1, 0, 0 - # X, A, B = x, a, b - - # def update(x, a, b): - # if x % 3 == 0: - # return (x*x) % N, (a*2) % n, (b*2) % n - # elif x % 3 == 1: - # return (x*alpha) % N, (a + 1) % n, b - # else: - # return (x*beta) % N, a, (b + 1) % n - - # for i in range(1, n): - # x, a, b = update(x, a, b) - # X, A, B = update(X, A, B) - # X, A, B = update(X, A, B) - # if x == X: - # break - - # return cls(a - A) / cls(B - b) diff --git a/galois/_lfsr.py b/galois/_lfsr.py index 852faced08..cc669c5435 100644 --- a/galois/_lfsr.py +++ b/galois/_lfsr.py @@ -63,6 +63,13 @@ class LFSR: In the Galois configuration, the next output is :math:`y = s_{n-1}` and the next state is computed by :math:`s_k = s_{n-1} g_n g_k + s_{k-1}`. In the case of :math:`s_0` there is no previous state added. + + References + ---------- + * https://core.ac.uk/download/pdf/288371609.pdf + * https://www.wseas.org/multimedia/journals/control/2018/a945903-022.pdf + * https://jhafranco.com/2014/02/15/n-ary-m-sequence-generator-in-python/ + * https://www.cs.uky.edu/~klapper/pdf/galois.pdf """ def __init__(self, poly, state=1, config="fibonacci"): @@ -88,12 +95,12 @@ def __init__(self, poly, state=1, config="fibonacci"): # Pre-compile the arithmetic functions and JIT routines if self.field._ufunc_mode != "python-calculate": - self._add = self.field._calculate_jit("add") - self._multiply = self.field._calculate_jit("multiply") + self._add = self.field._func_calculate("add") + self._multiply = self.field._func_calculate("multiply") self._step = jit_calculate(f"{self.config}_lfsr_step") else: - self._add = self.field._python_func("add") - self._multiply = self.field._python_func("multiply") + self._add = self.field._func_python("add") + self._multiply = self.field._func_python("multiply") self._step = python_func(f"{self.config}_lfsr_step") def __str__(self): @@ -331,18 +338,18 @@ def berlekamp_massey(sequence, config="fibonacci", state=False): if field.ufunc_mode != "python-calculate": sequence = sequence.astype(np.int64) - add = field._calculate_jit("add") - subtract = field._calculate_jit("subtract") - multiply = field._calculate_jit("multiply") - reciprocal = field._calculate_jit("reciprocal") + add = field._func_calculate("add") + subtract = field._func_calculate("subtract") + multiply = field._func_calculate("multiply") + reciprocal = field._func_calculate("reciprocal") coeffs = jit_calculate("berlekamp_massey")(sequence, add, subtract, multiply, reciprocal, field.characteristic, field.degree, field._irreducible_poly_int) coeffs = coeffs.astype(dtype) else: sequence = sequence.view(np.ndarray) - subtract = field._python_func("subtract") - multiply = field._python_func("multiply") - divide = field._python_func("divide") - reciprocal = field._python_func("reciprocal") + subtract = field._func_python("subtract") + multiply = field._func_python("multiply") + divide = field._func_python("divide") + reciprocal = field._func_python("reciprocal") coeffs = python_func("berlekamp_massey")(sequence, subtract, multiply, divide, reciprocal, field.characteristic, field.degree, field._irreducible_poly_int) if config == "fibonacci":