Merge pull request #66 from J08nY/fix/sympy-1.13

Fix sympy compat and enhance formula eval

8 files changed, 340 insertions, 125 deletions

diff --git a/pyecsca/ec/curve.py b/pyecsca/ec/curve.py
index 957aafc..ff89819 100644
--- a/pyecsca/ec/curve.py
+++ b/pyecsca/ec/curve.py
@@ -119,13 +119,14 @@ class EllipticCurve:
                 assumption_string = unparse(assumption).strip()
                 lhs, rhs = assumption_string.split(" = ")
                 expr = sympify(f"{rhs} - {lhs}")
-                for curve_param, value in self.parameters.items():
-                    expr = expr.subs(curve_param, k(value))
+                for symbol in expr.free_symbols:
+                    if (val := self.parameters.get(str(symbol), None)) is not None:
+                        expr = expr.xreplace({symbol: val})
                 if len(expr.free_symbols) > 0:
                     raise ValueError(
                         f"Missing necessary coordinate model parameter ({assumption_string})."
                     )
-                if k(expr) != 0:
+                if k.from_sympy(expr) != 0:
                     raise_unsatisified_assumption(
                         getconfig().ec.unsatisfied_coordinate_assumption_action,
                         f"Coordinate model {self.coordinate_model} has an unsatisifed assumption on the {param} parameter (0 = {expr})."
diff --git a/pyecsca/ec/formula/base.py b/pyecsca/ec/formula/base.py
index 704ba11..fee39aa 100644
--- a/pyecsca/ec/formula/base.py
+++ b/pyecsca/ec/formula/base.py
@@ -8,7 +8,7 @@ from astunparse import unparse
 from typing import List, Any, ClassVar, MutableMapping, Tuple, Union, Dict
 
 from public import public
-from sympy import FF, symbols, Poly, Rational
+from sympy import FF, symbols, Poly
 
 from pyecsca.ec.context import ResultAction
 from pyecsca.ec import context
@@ -16,7 +16,7 @@ from pyecsca.ec.error import UnsatisfiedAssumptionError, raise_unsatisified_assu
 from pyecsca.ec.mod import Mod, SymbolicMod
 from pyecsca.ec.op import CodeOp, OpType
 from pyecsca.misc.cfg import getconfig
-from pyecsca.misc.cache import sympify, simplify
+from pyecsca.misc.cache import sympify
 
 
 @public
@@ -108,6 +108,9 @@ class FormulaAction(ResultAction):
         return f"{self.__class__.__name__}({self.formula}, {self.input_points}) = {self.output_points}"
 
 
+_assumption_cache: Dict[Tuple[str, str, FF, Tuple[Mod, ...]], Mod] = {}
+
+
 @public
 class Formula(ABC):
     """Formula operating on points."""
@@ -157,76 +160,92 @@ class Formula(ABC):
                     )
                 params[coord + str(i + 1)] = value
 
+    def __validate_assumption_point(self, assumption, params):
+        # Handle an assumption check on value of input points.
+        alocals: Dict[str, Union[Mod, int]] = {**params}
+        compiled = compile(assumption, "", mode="eval")
+        holds = eval(compiled, None, alocals)  # eval is OK here, skipcq: PYL-W0123
+        return holds
+
+    def __validate_assumption_simple(self, lhs, rhs, field, params):
+        # Handle a simple parameter assignment (lhs is an unassigned parameter of the formula).
+        expr = sympify(rhs, evaluate=False)
+        used_symbols = sorted(expr.free_symbols)
+        used_params = []
+        for symbol in used_symbols:
+            if (value := params.get(str(symbol), None)) is not None:
+                used_params.append(value)
+                if isinstance(value, SymbolicMod):
+                    expr = expr.xreplace({symbol: value.x})
+                else:
+                    expr = expr.xreplace({symbol: int(value)})
+            else:
+                return False
+        cache_key = (lhs, rhs, field, tuple(used_params))
+        if cache_key in _assumption_cache:
+            params[lhs] = _assumption_cache[cache_key]
+        else:
+            if any(isinstance(x, SymbolicMod) for x in params.values()):
+                params[lhs] = SymbolicMod(expr, field)
+            else:
+                domain = FF(field)
+                numerator, denominator = expr.as_numer_denom()
+                val = int(domain.from_sympy(numerator) / domain.from_sympy(denominator))
+                params[lhs] = Mod(val, field)
+            _assumption_cache[cache_key] = params[lhs]
+        return True
+
+    def __validate_assumption_generic(self, lhs, rhs, field, params, assumption_string):
+        # Handle a generic parameter assignment (parameter may be anyway in the assumption).
+        expr = sympify(f"{rhs} - {lhs}", evaluate=False)
+        remaining = []
+        for symbol in expr.free_symbols:
+            if (value := params.get(str(symbol), None)) is not None:
+                if isinstance(value, SymbolicMod):
+                    expr = expr.xreplace({symbol: value.x})
+                else:
+                    expr = expr.xreplace({symbol: int(value)})
+            else:
+                remaining.append(symbol)
+        if len(remaining) > 1 or (param := str(remaining[0])) not in self.parameters:
+            raise ValueError(
+                f"This formula couldn't be executed due to an unsupported assumption ({assumption_string})."
+            )
+        numerator, _ = expr.as_numer_denom()
+        domain = FF(field)
+        poly = Poly(numerator, symbols(param), domain=domain)
+        roots = poly.ground_roots()
+        for root in roots:
+            params[param] = Mod(int(domain.from_sympy(root)), field)
+            return
+        raise UnsatisfiedAssumptionError(
+            f"Unsatisfied assumption in the formula ({assumption_string}).\n"
+            f"'{expr}' has no roots in the base field GF({field})."
+        )
+
     def __validate_assumptions(self, field, params):
         # Validate assumptions and compute formula parameters.
         # TODO: Should this also validate coordinate assumptions and compute their parameters?
-        is_symbolic = any(isinstance(x, SymbolicMod) for x in params.values())
         for assumption, assumption_string in zip(
             self.assumptions, self.assumptions_str
         ):
             lhs, rhs = assumption_string.split(" == ")
             if lhs in params:
-                # Handle an assumption check on value of input points.
-                alocals: Dict[str, Union[Mod, int]] = {**params}
-                compiled = compile(assumption, "", mode="eval")
-                holds = eval(compiled, None, alocals)  # eval is OK here, skipcq: PYL-W0123
-                if not holds:
-                    # The assumption doesn't hold, see what is the current configured action and do it.
+                if not self.__validate_assumption_point(assumption, params):
                     raise_unsatisified_assumption(
                         getconfig().ec.unsatisfied_formula_assumption_action,
                         f"Unsatisfied assumption in the formula ({assumption_string}).",
                     )
-            elif lhs in self.parameters and is_symbolic:
-                # Handle a symbolic assignment to a new parameter.
-                k = FF(field)
-                expr = sympify(rhs, evaluate=False)
-                for curve_param, value in params.items():
-                    if isinstance(value, SymbolicMod):
-                        expr = expr.subs(curve_param, value.x)
-                    else:
-                        expr = expr.subs(curve_param, k(value))
-                params[lhs] = SymbolicMod(expr, field)
-            else:
-                k = FF(field)
-                expr = sympify(f"{rhs} - {lhs}", evaluate=False)
-                for curve_param, value in params.items():
-                    if isinstance(value, SymbolicMod):
-                        expr = expr.subs(curve_param, value.x)
-                    else:
-                        expr = expr.subs(curve_param, k(value))
-                if (
-                    len(expr.free_symbols) > 1
-                    or (param := str(expr.free_symbols.pop())) not in self.parameters
-                ):
-                    raise ValueError(
-                        f"This formula couldn't be executed due to an unsupported assumption ({assumption_string})."
-                    )
-
-                def resolve(expression, k):
-                    if not expression.args:
-                        return expression
-                    args = []
-                    for arg in expression.args:
-                        if isinstance(arg, Rational):
-                            a = arg.p
-                            b = arg.q
-                            res = k(a) / k(b)
-                        else:
-                            res = resolve(arg, k)
-                        args.append(res)
-                    return expression.func(*args)
-
-                expr = resolve(simplify(expr), k)
-                poly = Poly(expr, symbols(param), domain=k)
-                roots = poly.ground_roots()
-                for root in roots:
-                    params[param] = Mod(int(root), field)
-                    break
-                else:
-                    raise UnsatisfiedAssumptionError(
-                        f"Unsatisfied assumption in the formula ({assumption_string}).\n"
-                        f"'{expr}' has no roots in the base field {k}."
+            elif lhs in self.parameters:
+                if not self.__validate_assumption_simple(lhs, rhs, field, params):
+                    raise_unsatisified_assumption(
+                        getconfig().ec.unsatisfied_formula_assumption_action,
+                        f"Unsatisfied assumption in the formula ({assumption_string}).",
                     )
+            else:
+                self.__validate_assumption_generic(
+                    lhs, rhs, field, params, assumption_string
+                )
 
     def __call__(self, field: int, *points: Any, **params: Mod) -> Tuple[Any, ...]:
         """
diff --git a/pyecsca/ec/formula/unroll.py b/pyecsca/ec/formula/unroll.py
index 83a7517..3c63e00 100644
--- a/pyecsca/ec/formula/unroll.py
+++ b/pyecsca/ec/formula/unroll.py
@@ -39,7 +39,7 @@ def unroll_formula_expr(formula: Formula) -> List[Tuple[str, Expr]]:
             # Handle a symbolic assignment to a new parameter.
             expr = sympify(rhs, evaluate=False)
             for curve_param, value in params.items():
-                expr = expr.subs(curve_param, value)
+                expr = expr.xreplace({curve_param: value})
             params[lhs] = expr
 
     locls = {**params, **inputs}
diff --git a/pyecsca/ec/key_generation.py b/pyecsca/ec/key_generation.py
index 13bf6d8..c583160 100644
--- a/pyecsca/ec/key_generation.py
+++ b/pyecsca/ec/key_generation.py
@@ -54,7 +54,7 @@ class KeyGeneration:
         """
         with KeygenAction(self.params) as action:
             privkey = Mod.random(self.params.order)
-            pubkey = self.mult.multiply(privkey.x)
+            pubkey = self.mult.multiply(int(privkey.x))
             if self.affine:
                 pubkey = pubkey.to_affine()
             return action.exit((privkey, pubkey))
diff --git a/pyecsca/ec/mod.py b/pyecsca/ec/mod.py
index 07b2734..4221a93 100644
--- a/pyecsca/ec/mod.py
+++ b/pyecsca/ec/mod.py
@@ -9,13 +9,19 @@ dispatches to the implementation chosen by the runtime configuration of the libr
 """
 import random
 import secrets
+import warnings
 from functools import wraps, lru_cache
 from typing import Type, Dict, Any, Tuple, Union
 
 from public import public
-from sympy import Expr, FF
+from sympy import Expr
 
-from pyecsca.ec.error import raise_non_invertible, raise_non_residue
+from pyecsca.ec.error import (
+    raise_non_invertible,
+    raise_non_residue,
+    NonResidueError,
+    NonResidueWarning,
+)
 from pyecsca.ec.context import ResultAction
 from pyecsca.misc.cfg import getconfig
 
@@ -28,8 +34,21 @@ except ImportError:
     gmpy2 = None
 
 
+has_flint = False
+try:
+    import flint
+
+    _major, _minor, *_ = flint.__version__.split(".")
+    if (int(_major), int(_minor)) >= (0, 5):
+        has_flint = True
+    else:
+        flint = None
+except ImportError:
+    flint = None
+
+
 @public
-def gcd(a, b):
+def gcd(a: int, b: int) -> int:
     """Euclid's greatest common denominator algorithm."""
     if abs(a) < abs(b):
         return gcd(b, a)
@@ -42,7 +61,7 @@ def gcd(a, b):
 
 
 @public
-def extgcd(a, b):
+def extgcd(a: int, b: int) -> Tuple[int, int, int]:
     """Compute the extended Euclid's greatest common denominator algorithm."""
     if abs(b) > abs(a):
         x, y, d = extgcd(b, a)
@@ -138,6 +157,7 @@ class RandomModAction(ResultAction):
 
 
 _mod_classes: Dict[str, Type] = {}
+_mod_order = ["gmp", "flint", "python"]
 
 
 @public
@@ -156,7 +176,10 @@ class Mod:
         selected_class = getconfig().ec.mod_implementation
         if selected_class not in _mod_classes:
             # Fallback to something
-            selected_class = next(iter(_mod_classes.keys()))
+            for fallback in _mod_order:
+                if fallback in _mod_classes:
+                    selected_class = fallback
+                    break
         return _mod_classes[selected_class].__new__(
             _mod_classes[selected_class], *args, **kwargs
         )
@@ -231,11 +254,6 @@ class Mod:
     def __rfloordiv__(self, other) -> "Mod":
         return ~self * other
 
-    @_check
-    def __divmod__(self, divisor) -> Tuple["Mod", "Mod"]:
-        q, r = divmod(self.x, divisor.x)
-        return self.__class__(q, self.n), self.__class__(r, self.n)
-
     def __bytes__(self) -> bytes:
         raise NotImplementedError
 
@@ -432,9 +450,6 @@ class Undefined(Mod):
     def __rfloordiv__(self, other):
         return NotImplemented
 
-    def __divmod__(self, divisor):
-        return NotImplemented
-
     def __bytes__(self):
         raise NotImplementedError
 
@@ -457,27 +472,6 @@ class Undefined(Mod):
         return NotImplemented
 
 
-@lru_cache
-def __ff_cache(n):
-    return FF(n)
-
-
-def _symbolic_check(func):
-    @wraps(func)
-    def method(self, other):
-        if type(self) is not type(other):
-            if type(other) is int:
-                other = self.__class__(__ff_cache(self.n)(other), self.n)
-            else:
-                other = self.__class__(other, self.n)
-        else:
-            if self.n != other.n:
-                raise ValueError
-        return func(self, other)
-
-    return method
-
-
 @public
 class SymbolicMod(Mod):
     """A symbolic element x of ℤₙ (implemented using sympy)."""
@@ -493,19 +487,19 @@ class SymbolicMod(Mod):
         self.x = x
         self.n = n
 
-    @_symbolic_check
+    @_check
     def __add__(self, other) -> "SymbolicMod":
         return self.__class__((self.x + other.x), self.n)
 
-    @_symbolic_check
+    @_check
     def __radd__(self, other) -> "SymbolicMod":
         return self + other
 
-    @_symbolic_check
+    @_check
     def __sub__(self, other) -> "SymbolicMod":
         return self.__class__((self.x - other.x), self.n)
 
-    @_symbolic_check
+    @_check
     def __rsub__(self, other) -> "SymbolicMod":
         return -self + other
 
@@ -527,33 +521,30 @@ class SymbolicMod(Mod):
     def __invert__(self) -> "SymbolicMod":
         return self.inverse()
 
-    @_symbolic_check
+    @_check
     def __mul__(self, other) -> "SymbolicMod":
         return self.__class__(self.x * other.x, self.n)
 
-    @_symbolic_check
+    @_check
     def __rmul__(self, other) -> "SymbolicMod":
         return self * other
 
-    @_symbolic_check
+    @_check
     def __truediv__(self, other) -> "SymbolicMod":
         return self * ~other
 
-    @_symbolic_check
+    @_check
     def __rtruediv__(self, other) -> "SymbolicMod":
         return ~self * other
 
-    @_symbolic_check
+    @_check
     def __floordiv__(self, other) -> "SymbolicMod":
         return self * ~other
 
-    @_symbolic_check
+    @_check
     def __rfloordiv__(self, other) -> "SymbolicMod":
         return ~self * other
 
-    def __divmod__(self, divisor) -> "SymbolicMod":
-        return NotImplemented
-
     def __bytes__(self):
         return int(self.x).to_bytes((self.n.bit_length() + 7) // 8, byteorder="big")
 
@@ -585,7 +576,7 @@ _mod_classes["symbolic"] = SymbolicMod
 if has_gmp:
 
     @lru_cache
-    def _is_prime(x) -> bool:
+    def _gmpy_is_prime(x) -> bool:
         return gmpy2.is_prime(x)
 
     @public
@@ -628,7 +619,7 @@ if has_gmp:
             return GMPMod(res, self.n, ensure=False)
 
         def is_residue(self) -> bool:
-            if not _is_prime(self.n):
+            if not _gmpy_is_prime(self.n):
                 raise NotImplementedError
             if self.x == 0:
                 return True
@@ -637,7 +628,7 @@ if has_gmp:
             return gmpy2.legendre(self.x, self.n) == 1
 
         def sqrt(self) -> "GMPMod":
-            if not _is_prime(self.n):
+            if not _gmpy_is_prime(self.n):
                 raise NotImplementedError
             if self.x == 0:
                 return GMPMod(gmpy2.mpz(0), self.n, ensure=False)
@@ -688,11 +679,6 @@ if has_gmp:
         def __mul__(self, other) -> "GMPMod":
             return GMPMod((self.x * other.x) % self.n, self.n, ensure=False)
 
-        @_check
-        def __divmod__(self, divisor) -> Tuple["GMPMod", "GMPMod"]:
-            q, r = gmpy2.f_divmod(self.x, divisor.x)
-            return GMPMod(q, self.n, ensure=False), GMPMod(r, self.n, ensure=False)
-
         def __bytes__(self):
             return int(self.x).to_bytes((self.n.bit_length() + 7) // 8, byteorder="big")
 
@@ -729,3 +715,210 @@ if has_gmp:
             )
 
     _mod_classes["gmp"] = GMPMod
+
+
+if has_flint:
+
+    @lru_cache
+    def _fmpz_ctx(n: Union[int, flint.fmpz_mod_ctx]) -> flint.fmpz_mod_ctx:
+        if type(n) is flint.fmpz_mod_ctx:
+            return n
+        return flint.fmpz_mod_ctx(n)
+
+    @lru_cache
+    def _fmpz_is_prime(x: flint.fmpz) -> bool:
+        return x.is_probable_prime()
+
+    def _flint_check(func):
+        @wraps(func)
+        def method(self, other):
+            if self.__class__ is not type(other):
+                other = self.__class__(other, self.n)
+            elif self._ctx != other._ctx:
+                raise ValueError
+            return func(self, other)
+
+        return method
+
+    @public
+    class FlintMod(Mod):
+        """An element x of ℤₙ. Implemented by GMP."""
+
+        x: flint.fmpz_mod
+        _ctx: flint.fmpz_mod_ctx
+        __slots__ = ("x", "_ctx")
+
+        def __new__(cls, *args, **kwargs):
+            return object.__new__(cls)
+
+        def __init__(
+            self,
+            x: Union[int, flint.fmpz_mod],
+            n: Union[int, flint.fmpz_mod_ctx],
+            ensure: bool = True,
+        ):
+            if ensure:
+                self._ctx = _fmpz_ctx(n)
+                self.x = self._ctx(x)
+            else:
+                self._ctx = n
+                self.x = x
+
+        @property
+        def n(self) -> flint.fmpz:
+            return self._ctx.modulus()
+
+        def bit_length(self):
+            return int(self.x).bit_length()
+
+        def inverse(self) -> "FlintMod":
+            if self.x == 0:
+                raise_non_invertible()
+            if self.x == 1:
+                return FlintMod(self._ctx(1), self._ctx, ensure=False)
+            try:
+                res = self.x.inverse()
+            except ZeroDivisionError:
+                raise_non_invertible()
+                res = self._ctx(0)
+            return FlintMod(res, self._ctx, ensure=False)
+
+        def is_residue(self) -> bool:
+            try:
+                with warnings.catch_warnings(record=True) as warns:
+                    self.sqrt()
+                if warns and isinstance(warns[0], NonResidueWarning):
+                    return False
+            except NonResidueError:
+                return False
+            return True
+
+        def sqrt(self) -> "FlintMod":
+            mod = self.n
+            if not _fmpz_is_prime(mod):
+                raise NotImplementedError
+            try:
+                res = flint.fmpz(int(self.x)).sqrtmod(mod)
+                return FlintMod(self._ctx(res), self._ctx, ensure=False)
+            except ValueError:
+                raise_non_residue()
+
+            if mod % 4 == 3:
+                return self ** int((mod + 1) // 4)
+            q = mod - 1
+            s = 0
+            while q % 2 == 0:
+                q //= 2
+                s += 1
+
+            z = self._ctx(2)
+            while FlintMod(z, self._ctx, ensure=False).is_residue():
+                z += 1
+
+            m = s
+            c = FlintMod(z, self._ctx, ensure=False) ** int(q)
+            t = self ** int(q)
+            r_exp = (q + 1) // 2
+            r = self ** int(r_exp)
+
+            while t != 1:
+                i = 1
+                while not (t ** (2**i)) == 1:
+                    i += 1
+                two_exp = m - (i + 1)
+                b = c ** int(FlintMod(self._ctx(2), self._ctx, ensure=False) ** two_exp)
+                m = int(FlintMod(self._ctx(i), self._ctx, ensure=False))
+                c = b**2
+                t *= c
+                r *= b
+            return r
+
+        @_flint_check
+        def __add__(self, other) -> "FlintMod":
+            return FlintMod(self.x + other.x, self._ctx, ensure=False)
+
+        @_flint_check
+        def __radd__(self, other) -> "Mod":
+            return self + other
+
+        @_flint_check
+        def __sub__(self, other) -> "FlintMod":
+            return FlintMod(self.x - other.x, self._ctx, ensure=False)
+
+        @_flint_check
+        def __rsub__(self, other) -> "Mod":
+            return -self + other
+
+        def __neg__(self) -> "FlintMod":
+            return FlintMod(-self.x, self._ctx, ensure=False)
+
+        @_flint_check
+        def __mul__(self, other) -> "FlintMod":
+            return FlintMod(self.x * other.x, self._ctx, ensure=False)
+
+        @_flint_check
+        def __rmul__(self, other) -> "Mod":
+            return self * other
+
+        @_flint_check
+        def __truediv__(self, other) -> "Mod":
+            return self * ~other
+
+        @_flint_check
+        def __rtruediv__(self, other) -> "Mod":
+            return ~self * other
+
+        @_flint_check
+        def __floordiv__(self, other) -> "Mod":
+            return self * ~other
+
+        @_flint_check
+        def __rfloordiv__(self, other) -> "Mod":
+            return ~self * other
+
+        def __bytes__(self):
+            return int(self.x).to_bytes(
+                (int(self.n).bit_length() + 7) // 8, byteorder="big"
+            )
+
+        def __int__(self):
+            return int(self.x)
+
+        def __eq__(self, other):
+            if type(other) is int:
+                return self.x == other
+            if type(other) is not FlintMod:
+                return False
+            try:
+                return self.x == other.x
+            except ValueError:
+                return False
+
+        def __ne__(self, other):
+            return not self == other
+
+        def __repr__(self):
+            return str(int(self.x))
+
+        def __hash__(self):
+            return hash(("FlintMod", self.x, self.n))
+
+        def __pow__(self, n) -> "FlintMod":
+            if type(n) not in (int, flint.fmpz):
+                raise TypeError
+            if n == 0:
+                return FlintMod(self._ctx(1), self._ctx, ensure=False)
+            if n < 0:
+                return self.inverse() ** (-n)
+            if n == 1:
+                return FlintMod(self.x, self._ctx, ensure=False)
+            return FlintMod(self.x**n, self._ctx, ensure=False)
+
+        def __getstate__(self):
+            return {"x": int(self.x), "n": int(self.n)}
+
+        def __setstate__(self, state):
+            self._ctx = _fmpz_ctx(state["n"])
+            self.x = self._ctx(state["x"])
+
+    _mod_classes["flint"] = FlintMod
diff --git a/pyecsca/ec/params.py b/pyecsca/ec/params.py
index d9e05e5..6693581 100644
--- a/pyecsca/ec/params.py
+++ b/pyecsca/ec/params.py
@@ -220,7 +220,7 @@ def _create_params(curve, coords, infty):
                 lhs, rhs = assumption_string.split(" = ")
                 expr = sympify(f"{rhs} - {lhs}")
                 for curve_param, value in params.items():
-                    expr = expr.subs(curve_param, k(value))
+                    expr = expr.subs(curve_param, value)
                 if (
                     len(expr.free_symbols) > 1
                     or (param := str(expr.free_symbols.pop()))
@@ -229,10 +229,11 @@ def _create_params(curve, coords, infty):
                     raise ValueError(
                         f"This coordinate model couldn't be loaded due to an unsupported assumption ({assumption_string})."
                     )
-                poly = Poly(expr, symbols(param), domain=k)
+                numerator, _ = expr.as_numer_denom()
+                poly = Poly(numerator, symbols(param), domain=k)
                 roots = poly.ground_roots()
                 for root in roots:
-                    params[param] = Mod(int(root), field)
+                    params[param] = Mod(int(k.from_sympy(root)), field)
                     break
                 else:
                     raise_unsatisified_assumption(
diff --git a/pyecsca/misc/cfg.py b/pyecsca/misc/cfg.py
index 05f89b8..037074d 100644
--- a/pyecsca/misc/cfg.py
+++ b/pyecsca/misc/cfg.py
@@ -110,6 +110,7 @@ class ECConfig:
         One of:
 
          - ``"gmp"``: Requires the GMP library and `gmpy2` package.
+         - ``"flint"``: Requires the flint library and `python-flint` package.
          - ``"python"``: Doesn't require anything.
          - ``"symbolic"``: Requires sympy.
         """
@@ -117,8 +118,8 @@ class ECConfig:
 
     @mod_implementation.setter
     def mod_implementation(self, value: str):
-        if value not in ("python", "gmp", "symbolic"):
-            raise ValueError("Bad Mod implementaiton, can be one of 'python', 'gmp' or 'symbolic'.")
+        if value not in ("python", "gmp", "flint", "symbolic"):
+            raise ValueError("Bad Mod implementaiton, can be one of 'python', 'gmp', 'flint' or 'symbolic'.")
         self._mod_implementation = value
 
 
diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index a697388..8005382 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -180,7 +180,7 @@ def compute_factor_set(
         # Go over all the factors of the intermediate, forget the power
         for factor, power in factor_list:
             # Remove unnecessary variables from the Poly
-            reduced = factor.exclude()
+            reduced = factor.exclude() if not factor.is_univariate else factor
             # If there are only curve parameters, we do not care about the polynomial
             if set(reduced.gens).issubset(curve_params):  # type: ignore[attr-defined]
                 continue