galois
galois copied to clipboard
mhostetter
Pure-Python arithmetic occasionally has issues
When using overridden ufuncs (like np.bitwise_xor for addition in binary fields) sometimes there are type errors.
In [1]: import galois
In [2]: c = galois.irreducible_poly(2**100, 3); c
Out[2]: Poly(x^3 + 2, GF(2^100))
In [3]: lfsr = galois.FLFSR(c.reverse()); lfsr
Out[3]: <Fibonacci LFSR: f(x) = 2x^3 + 1 over GF(2^100)>
In [4]: print(lfsr)
Fibonacci LFSR:
field: GF(2^100)
feedback_poly: 2x^3 + 1
characteristic_poly: x^3 + 2
taps: [0, 0, 2]
order: 3
state: [1, 1, 1]
initial_state: [1, 1, 1]
In [5]: lfsr.step(10)
Out[5]: GF([1, 1, 1, 2, 2, 2, 4, 4, 4, 8], order=2^100)
In [6]: lfsr.step(10)
Out[6]: GF([ 8, 8, 16, 16, 16, 32, 32, 32, 64, 64], order=2^100)
In [7]: lfsr.step(10)
Out[7]: GF([ 64, 128, 128, 128, 256, 256, 256, 512, 512, 512], order=2^100)
In [8]: lfsr.step(10)
Out[8]: GF([1024, 1024, 1024, 2048, 2048, 2048, 4096, 4096, 4096, 8192], order=2^100)
In [9]: lfsr.step(10)
Out[9]:
GF([ 8192, 8192, 16384, 16384, 16384, 32768, 32768, 32768, 65536, 65536],
order=2^100)
In [10]: lfsr.step(10)
Out[10]:
GF([ 65536, 131072, 131072, 131072, 262144, 262144, 262144, 524288,
524288, 524288], order=2^100)
In [11]: lfsr.step(10)
Out[11]:
GF([1048576, 1048576, 1048576, 2097152, 2097152, 2097152, 4194304,
4194304, 4194304, 8388608], order=2^100)
In [12]: lfsr.step(10)
Out[12]:
GF([ 8388608, 8388608, 16777216, 16777216, 16777216, 33554432, 33554432,
33554432, 67108864, 67108864], order=2^100)
In [13]: lfsr.step(10)
Out[13]:
GF([ 67108864, 134217728, 134217728, 134217728, 268435456, 268435456,
268435456, 536870912, 536870912, 536870912], order=2^100)
In [14]: lfsr.step(10)
Out[14]:
GF([1073741824, 1073741824, 1073741824, 2147483648, 2147483648,
2147483648, 4294967296, 4294967296, 4294967296, 8589934592],
order=2^100)
In [15]: lfsr.step(10)
Out[15]:
GF([ 8589934592, 8589934592, 17179869184, 17179869184, 17179869184,
34359738368, 34359738368, 34359738368, 68719476736, 68719476736],
order=2^100)
In [16]: lfsr.step(10)
Out[16]:
GF([ 68719476736, 137438953472, 137438953472, 137438953472, 274877906944,
274877906944, 274877906944, 549755813888, 549755813888, 549755813888],
order=2^100)
In [17]: lfsr.step(10)
Out[17]:
GF([1099511627776, 1099511627776, 1099511627776, 2199023255552,
2199023255552, 2199023255552, 4398046511104, 4398046511104,
4398046511104, 8796093022208], order=2^100)
In [18]: lfsr.step(10)
Out[18]:
GF([ 8796093022208, 8796093022208, 17592186044416, 17592186044416,
17592186044416, 35184372088832, 35184372088832, 35184372088832,
70368744177664, 70368744177664], order=2^100)
In [19]: lfsr.step(10)
Out[19]:
GF([ 70368744177664, 140737488355328, 140737488355328, 140737488355328,
281474976710656, 281474976710656, 281474976710656, 562949953421312,
562949953421312, 562949953421312], order=2^100)
In [20]: lfsr.step(10)
Out[20]:
GF([1125899906842624, 1125899906842624, 1125899906842624,
2251799813685248, 2251799813685248, 2251799813685248,
4503599627370496, 4503599627370496, 4503599627370496,
9007199254740992], order=2^100)
In [21]: lfsr.step(10)
Out[21]:
GF([ 9007199254740992, 9007199254740992, 18014398509481984,
18014398509481984, 18014398509481984, 36028797018963968,
36028797018963968, 36028797018963968, 72057594037927936,
72057594037927936], order=2^100)
In [22]: lfsr.step(10)
Out[22]:
GF([ 72057594037927936, 144115188075855872, 144115188075855872,
144115188075855872, 288230376151711744, 288230376151711744,
288230376151711744, 576460752303423488, 576460752303423488,
576460752303423488], order=2^100)
In [23]: lfsr.step(10)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [23], in <cell line: 1>()
----> 1 lfsr.step(10)
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:519, in FLFSR.step(self, steps)
447 """
448 Produces the next `steps` output symbols.
449
(...)
516 lfsr.state
517 """
518 # pylint: disable=useless-super-delegation
--> 519 return super().step(steps)
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:106, in _LFSR.step(self, steps)
104 return self.field([])
105 elif steps > 0:
--> 106 return self._step_forward(steps)
107 else:
108 return self._step_backward(abs(steps))
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:114, in _LFSR._step_forward(self, steps)
111 assert steps > 0
113 if self._type == "fibonacci":
--> 114 y, state = fibonacci_lfsr_step_forward_jit(self.field)(self.taps, self.state, steps)
115 else:
116 y, state = galois_lfsr_step_forward_jit(self.field)(self.taps, self.state, steps)
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:734, in fibonacci_lfsr_step_forward_jit.__call__(self, taps, state, steps)
732 else:
733 state_ = state.view(np.ndarray) # NOTE: This will be modified
--> 734 y = self.python(taps.view(np.ndarray), state_, steps)
735 y = self.field._view(y)
737 return y, state_
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:755, in fibonacci_lfsr_step_forward_jit.implementation(taps, state, steps)
753 f = 0 # The feedback value
754 for j in range(n):
--> 755 f = ADD(f, MULTIPLY(state[j], taps[j]))
757 y[i] = state[-1] # Output is popped off the shift register
758 state[1:] = state[0:-1] # Shift state rightward
TypeError: ufunc 'bitwise_xor' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
It seems replacing
f = ADD(f, MULTIPLY(state[j], taps[j]))
with
f = ADD(f, MULTIPLY(state[j], taps[j], dtype=dtype), dtype=dtype)
resolves the issue. The issue is if the two inputs are small integers with type int, the output of np.bitwise_xor() will be np.int64. As the loop continues, eventually the np.int64 will overflow, since values can go to 2^100. Explicitly specifying the output dtype seems to resolve this.
This probably should be done everywhere ufuncs are used in JIT-ed functions.
Unfortunately the dtype solution doesn't work for custom ufuncs.
See below for details
In [6]: GF = galois.GF(2**10, compile="jit-calculate")
In [7]: c = galois.irreducible_poly(2**10, 3); c
Out[7]: Poly(x^3 + 2, GF(2^10))
In [8]: lfsr = galois.FLFSR(c.reverse()); lfsr
Out[8]: <Fibonacci LFSR: f(x) = 2x^3 + 1 over GF(2^10)>
In [9]: lfsr.step(10)
---------------------------------------------------------------------------
TypingError Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 lfsr.step(10)
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:519, in FLFSR.step(self, steps)
447 """
448 Produces the next `steps` output symbols.
449
(...)
516 lfsr.state
517 """
518 # pylint: disable=useless-super-delegation
--> 519 return super().step(steps)
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:106, in _LFSR.step(self, steps)
104 return self.field([])
105 elif steps > 0:
--> 106 return self._step_forward(steps)
107 else:
108 return self._step_backward(abs(steps))
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:114, in _LFSR._step_forward(self, steps)
111 assert steps > 0
113 if self._type == "fibonacci":
--> 114 y, state = fibonacci_lfsr_step_forward_jit(self.field)(self.taps, self.state, steps)
115 else:
116 y, state = galois_lfsr_step_forward_jit(self.field)(self.taps, self.state, steps)
File /mnt/c/Users/matth/repos/galois/galois/_lfsr.py:730, in fibonacci_lfsr_step_forward_jit.__call__(self, taps, state, steps)
728 if self.field.ufunc_mode != "python-calculate":
729 state_ = state.astype(np.int64) # NOTE: This will be modified
--> 730 y = self.jit(taps.astype(np.int64), state_, steps)
731 y = y.astype(state.dtype)
732 else:
File /mnt/c/Users/matth/repos/galois/galois/_domains/_function.py:78, in Function.jit(self)
76 if key_2 not in self._CACHE[key_1]:
77 self.set_globals() # Set the globals once before JIT compiling the function
---> 78 self._CACHE[key_1][key_2] = numba.jit(self._SIGNATURE.signature, nopython=True)(self.implementation)
80 return self._CACHE[key_1][key_2]
File ~/.local/lib/python3.8/site-packages/numba/core/decorators.py:219, in _jit.<locals>.wrapper(func)
217 with typeinfer.register_dispatcher(disp):
218 for sig in sigs:
--> 219 disp.compile(sig)
220 disp.disable_compile()
221 return disp
File ~/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:965, in Dispatcher.compile(self, sig)
963 with ev.trigger_event("numba:compile", data=ev_details):
964 try:
--> 965 cres = self._compiler.compile(args, return_type)
966 except errors.ForceLiteralArg as e:
967 def folded(args, kws):
File ~/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:129, in _FunctionCompiler.compile(self, args, return_type)
127 return retval
128 else:
--> 129 raise retval
File ~/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:139, in _FunctionCompiler._compile_cached(self, args, return_type)
136 pass
138 try:
--> 139 retval = self._compile_core(args, return_type)
140 except errors.TypingError as e:
141 self._failed_cache[key] = e
File ~/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:152, in _FunctionCompiler._compile_core(self, args, return_type)
149 flags = self._customize_flags(flags)
151 impl = self._get_implementation(args, {})
--> 152 cres = compiler.compile_extra(self.targetdescr.typing_context,
153 self.targetdescr.target_context,
154 impl,
155 args=args, return_type=return_type,
156 flags=flags, locals=self.locals,
157 pipeline_class=self.pipeline_class)
158 # Check typing error if object mode is used
159 if cres.typing_error is not None and not flags.enable_pyobject:
File ~/.local/lib/python3.8/site-packages/numba/core/compiler.py:693, in compile_extra(typingctx, targetctx, func, args, return_type, flags, locals, library, pipeline_class)
669 """Compiler entry point
670
671 Parameter
(...)
689 compiler pipeline
690 """
691 pipeline = pipeline_class(typingctx, targetctx, library,
692 args, return_type, flags, locals)
--> 693 return pipeline.compile_extra(func)
File ~/.local/lib/python3.8/site-packages/numba/core/compiler.py:429, in CompilerBase.compile_extra(self, func)
427 self.state.lifted = ()
428 self.state.lifted_from = None
--> 429 return self._compile_bytecode()
File ~/.local/lib/python3.8/site-packages/numba/core/compiler.py:497, in CompilerBase._compile_bytecode(self)
493 """
494 Populate and run pipeline for bytecode input
495 """
496 assert self.state.func_ir is None
--> 497 return self._compile_core()
File ~/.local/lib/python3.8/site-packages/numba/core/compiler.py:476, in CompilerBase._compile_core(self)
474 self.state.status.fail_reason = e
475 if is_final_pipeline:
--> 476 raise e
477 else:
478 raise CompilerError("All available pipelines exhausted")
File ~/.local/lib/python3.8/site-packages/numba/core/compiler.py:463, in CompilerBase._compile_core(self)
461 res = None
462 try:
--> 463 pm.run(self.state)
464 if self.state.cr is not None:
465 break
File ~/.local/lib/python3.8/site-packages/numba/core/compiler_machinery.py:353, in PassManager.run(self, state)
350 msg = "Failed in %s mode pipeline (step: %s)" % \
351 (self.pipeline_name, pass_desc)
352 patched_exception = self._patch_error(msg, e)
--> 353 raise patched_exception
File ~/.local/lib/python3.8/site-packages/numba/core/compiler_machinery.py:341, in PassManager.run(self, state)
339 pass_inst = _pass_registry.get(pss).pass_inst
340 if isinstance(pass_inst, CompilerPass):
--> 341 self._runPass(idx, pass_inst, state)
342 else:
343 raise BaseException("Legacy pass in use")
File ~/.local/lib/python3.8/site-packages/numba/core/compiler_lock.py:35, in _CompilerLock.__call__.<locals>._acquire_compile_lock(*args, **kwargs)
32 @functools.wraps(func)
33 def _acquire_compile_lock(*args, **kwargs):
34 with self:
---> 35 return func(*args, **kwargs)
File ~/.local/lib/python3.8/site-packages/numba/core/compiler_machinery.py:296, in PassManager._runPass(self, index, pss, internal_state)
294 mutated |= check(pss.run_initialization, internal_state)
295 with SimpleTimer() as pass_time:
--> 296 mutated |= check(pss.run_pass, internal_state)
297 with SimpleTimer() as finalize_time:
298 mutated |= check(pss.run_finalizer, internal_state)
File ~/.local/lib/python3.8/site-packages/numba/core/compiler_machinery.py:269, in PassManager._runPass.<locals>.check(func, compiler_state)
268 def check(func, compiler_state):
--> 269 mangled = func(compiler_state)
270 if mangled not in (True, False):
271 msg = ("CompilerPass implementations should return True/False. "
272 "CompilerPass with name '%s' did not.")
File ~/.local/lib/python3.8/site-packages/numba/core/typed_passes.py:105, in BaseTypeInference.run_pass(self, state)
99 """
100 Type inference and legalization
101 """
102 with fallback_context(state, 'Function "%s" failed type inference'
103 % (state.func_id.func_name,)):
104 # Type inference
--> 105 typemap, return_type, calltypes, errs = type_inference_stage(
106 state.typingctx,
107 state.targetctx,
108 state.func_ir,
109 state.args,
110 state.return_type,
111 state.locals,
112 raise_errors=self._raise_errors)
113 state.typemap = typemap
114 # save errors in case of partial typing
File ~/.local/lib/python3.8/site-packages/numba/core/typed_passes.py:83, in type_inference_stage(typingctx, targetctx, interp, args, return_type, locals, raise_errors)
81 infer.build_constraint()
82 # return errors in case of partial typing
---> 83 errs = infer.propagate(raise_errors=raise_errors)
84 typemap, restype, calltypes = infer.unify(raise_errors=raise_errors)
86 # Output all Numba warnings
File ~/.local/lib/python3.8/site-packages/numba/core/typeinfer.py:1086, in TypeInferer.propagate(self, raise_errors)
1083 force_lit_args = [e for e in errors
1084 if isinstance(e, ForceLiteralArg)]
1085 if not force_lit_args:
-> 1086 raise errors[0]
1087 else:
1088 raise reduce(operator.or_, force_lit_args)
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<numba._DUFunc 'calculate'>) found for signature:
>>> calculate(int64, int64, dtype=dtype(int64))
There are 2 candidate implementations:
- Of which 2 did not match due to:
Overload in function 'DUFunc._type_me': File: numba/np/ufunc/dufunc.py: Line 264.
With argument(s): '(int64, int64, dtype=dtype(int64))':
Rejected as the implementation raised a specific error:
AssertionError:
raised from /home/matt/.local/lib/python3.8/site-packages/numba/np/ufunc/dufunc.py:272
During: resolving callee type: Function(<numba._DUFunc 'calculate'>)
During: typing of call at /mnt/c/Users/matth/repos/galois/galois/_lfsr.py (756)
File "galois/_lfsr.py", line 756:
def implementation(taps, state, steps):
<source elided>
for j in range(n):
f = ADD(f, MULTIPLY(state[j], taps[j], dtype=dtype), dtype=dtype)
^
I also tried
class UFunc:
# ...
@property
def python_calculate(self) -> Callable:
"""
A pure-Python ufunc implemented using explicit calculation.
"""
if self.override:
return self.override # From
return lambda *args, **kwargs: self.override(*args, **kwargs, dtype=np.object_) # To
which works in newer NumPy versions and solves the issues. However, it apparently doesn't work in NumPy 1.20.3.
>>> np.bitwise_xor(2, 3, dtype=object)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: No loop matching the specified signature and casting was found for ufunc bitwise_xor
I believe it was added/modified in NumPy 1.21, see https://numpy.org/doc/stable/release/1.21.0-notes.html#ufunc-signature-and-dtype-generalization-and-casting
Treating the inputs to multiply as arrays (not scalars) works, but is ugly.
class fibonacci_lfsr_step_forward_jit(Function):
@staticmethod
def implementation(taps, state, steps):
n = taps.size
y = np.zeros(steps, dtype=state.dtype) # The output array
for i in range(steps):
f = np.zeros(1, dtype=state.dtype) # The feedback value
for j in range(n):
f[0:1] = ADD(f[0:1], MULTIPLY(state[j:j+1], taps[j:j+1]))
y[i] = state[-1] # Output is popped off the shift register
state[1:] = state[0:-1] # Shift state rightward
state[0] = f[0] # Insert feedback value at leftmost position
return y