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

Fix sympy compat and enhance formula eval

21 files changed, 495 insertions, 210 deletions

diff --git a/.github/workflows/perf.yml b/.github/workflows/perf.yml
index 6211ea1..69ea332 100644
--- a/.github/workflows/perf.yml
+++ b/.github/workflows/perf.yml
@@ -6,6 +6,7 @@ env:
   LLVM_CONFIG: /usr/bin/llvm-config-10
   PS_PACKAGES: libps4000 libps5000 libps6000
   GMP_PACKAGES: libgmp-dev libmpfr-dev libmpc-dev
+  FLINT_PACKAGES: libflint-dev
   OTHER_PACKAGES: swig gcc libpcsclite-dev llvm-10 libllvm10 llvm-10-dev libpari-dev pari-gp pari-seadata
 
 jobs:
@@ -14,10 +15,10 @@ jobs:
     strategy:
       matrix:
         python-version: ["3.9", "3.10", "3.11"]
-        gmp: [0, 1]
+        mod: ["python", "gmp", "flint"]
     env:
       PYTHON: ${{ matrix.python-version }}
-      USE_GMP: ${{ matrix.gmp }}
+      MOD_IMPL: ${{ matrix.mod }}
     steps:
       - uses: actions/checkout@v4
         with:
@@ -25,10 +26,10 @@ jobs:
       - uses: actions/cache@v4
         with:
           path: ~/.cache/pip
-          key: pip-${{ runner.os }}-${{ matrix.gmp }}-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
+          key: pip-${{ runner.os }}-${{ matrix.mod }}-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
           restore-keys: |
-            pip-${{ runner.os }}-${{ matrix.gmp }}-${{ matrix.python-version }}-
-            pip-${{ runner.os }}-${{ matrix.gmp }}-
+            pip-${{ runner.os }}-${{ matrix.mod }}-${{ matrix.python-version }}-
+            pip-${{ runner.os }}-${{ matrix.mod }}-
             pip-${{ runner.os }}-
       - name: Setup Python ${{ matrix.python-version }}
         uses: actions/setup-python@v5
@@ -42,22 +43,45 @@ jobs:
       - name: Install system dependencies
         run: |
           sudo apt-get install -y $PS_PACKAGES $OTHER_PACKAGES
-          if [ $USE_GMP == 1 ]; then sudo apt-get install -y $GMP_PACKAGES; fi
+          if [ $MOD_IMPL == "gmp" ]; then sudo apt-get install -y $GMP_PACKAGES; fi
+          if [ $MOD_IMPL == "flint" ]; then sudo apt-get install -y $FLINT_PACKAGES; fi
       - name: Install picoscope bindings
         run: |
           python -m pip install -U pip setuptools wheel
           git clone https://github.com/colinoflynn/pico-python && cd pico-python && pip install . && cd ..
           git clone https://github.com/picotech/picosdk-python-wrappers && cd picosdk-python-wrappers && pip install . && cd ..
-      - name: Install dependencies
+      - name: Install
         run: |
-          if [ $USE_GMP == 1 ]; then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, gmp, test, dev]"; fi
-          if [ $USE_GMP == 0 ]; then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, test, dev]"; fi
+          if [ $MOD_IMPL == "gmp" ];    then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, gmp, test, dev]"; fi
+          if [ $MOD_IMPL == "flint" ];  then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, flint, test, dev]"; fi
+          if [ $MOD_IMPL == "python" ]; then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, test, dev]"; fi
       - name: Perf
         run: |
           make perf
       - name: Archive perf results
         uses: actions/upload-artifact@v4
         with:
-          name: perf-results-${{ matrix.gmp }}-${{ matrix.python-version }}
+          name: perf-results-${{ matrix.mod }}-${{ matrix.python-version }}
           path:
             .perf
+  merge:
+    runs-on: ubuntu-20.04
+    needs: perf
+    steps:
+      - name: Download perf results
+        uses: actions/download-artifact@v4
+      - name: Merge
+        run: |
+          mkdir out
+          for dir in */; do if [ "$dir" != "out/" ]; then echo $dir; for f in "$dir"*; do fname=$(basename $f); echo $fname; cat $f >> out/$fname; done; fi; done
+      - name: Upload merged
+        uses: actions/upload-artifact@v4
+        with:
+          name: perf-results
+          path:
+            out
+      - name: Delete old
+        uses: geekyeggo/delete-artifact@v5
+        with:
+          name: |
+              perf-results-*
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 8a31b88..1134792 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -6,6 +6,7 @@ env:
   LLVM_CONFIG: /usr/bin/llvm-config-10
   PS_PACKAGES: libps4000 libps5000 libps6000
   GMP_PACKAGES: libgmp-dev libmpfr-dev libmpc-dev
+  FLINT_PACKAGES: libflint-dev
   OTHER_PACKAGES: swig gcc libpcsclite-dev llvm-10 libllvm10 llvm-10-dev libpari-dev pari-gp pari-seadata
 
 jobs:
@@ -14,10 +15,10 @@ jobs:
     strategy:
       matrix:
         python-version: ["3.9", "3.10", "3.11"]
-        gmp: [0, 1]
+        mod: ["python", "gmp", "flint"]
     env:
       PYTHON: ${{ matrix.python-version }}
-      USE_GMP: ${{ matrix.gmp }}
+      MOD_IMPL: ${{ matrix.mod }}
     steps:
       - uses: actions/checkout@v4
         with:
@@ -25,10 +26,10 @@ jobs:
       - uses: actions/cache@v4
         with:
           path: ~/.cache/pip
-          key: pip-${{ runner.os }}-${{ matrix.gmp }}-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
+          key: pip-${{ runner.os }}-${{ matrix.mod }}-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
           restore-keys: |
-            pip-${{ runner.os }}-${{ matrix.gmp }}-${{ matrix.python-version }}-
-            pip-${{ runner.os }}-${{ matrix.gmp }}-
+            pip-${{ runner.os }}-${{ matrix.mod }}-${{ matrix.python-version }}-
+            pip-${{ runner.os }}-${{ matrix.mod }}-
             pip-${{ runner.os }}-
       - name: Setup Python ${{ matrix.python-version }}
         uses: actions/setup-python@v5
@@ -42,16 +43,18 @@ jobs:
       - name: Install system dependencies
         run: |
           sudo apt-get install -y $PS_PACKAGES $OTHER_PACKAGES
-          if [ $USE_GMP == 1 ]; then sudo apt-get install -y $GMP_PACKAGES; fi
+          if [ $MOD_IMPL == "gmp" ]; then sudo apt-get install -y $GMP_PACKAGES; fi
+          if [ $MOD_IMPL == "flint" ]; then sudo apt-get install -y $FLINT_PACKAGES; fi
       - name: Install picoscope bindings
         run: |
           python -m pip install -U pip setuptools wheel
           git clone https://github.com/colinoflynn/pico-python && cd pico-python && pip install . && cd ..
           git clone https://github.com/picotech/picosdk-python-wrappers && cd picosdk-python-wrappers && pip install . && cd ..
-      - name: Install dependencies
+      - name: Install
         run: |
-          if [ $USE_GMP == 1 ]; then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, gmp, test, dev]"; fi
-          if [ $USE_GMP == 0 ]; then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, test, dev]"; fi
+          if [ $MOD_IMPL == "gmp" ];    then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, gmp, test, dev]"; fi
+          if [ $MOD_IMPL == "flint" ];  then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, flint, test, dev]"; fi
+          if [ $MOD_IMPL == "python" ]; then pip install -e ".[picoscope_sdk, picoscope_alt, chipwhisperer, smartcard, pari, leia, test, dev]"; fi
       - name: Test
         run: |
           make test
diff --git a/docs/installation.rst b/docs/installation.rst
index d973a02..2396171 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -40,6 +40,7 @@ Requirements
        - **Faster arithmetic:**
 
          - gmpy2_ (and also GMP library)
+         - python-flint_ (and also Flint library)
          - cypari2_ (and also PARI library)
 
     *pyecsca* contains data from the `Explicit-Formulas Database`_ by Daniel J. Bernstein and Tanja Lange.
@@ -92,6 +93,7 @@ Requirements
 .. _pyscard: https://pyscard.sourceforge.io/
 .. _leia: https://pypi.org/project/smartleia/
 .. _gmpy2: https://gmpy2.readthedocs.io/
+.. _python-flint: https://fredrikj.net/python-flint/
 .. _cypari2: https://cypari2.readthedocs.io/
 .. _pytest: https://pytest.org
 .. _mypy: http://mypy-lang.org/
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
diff --git a/pyproject.toml b/pyproject.toml
index c6081ab..0c0ef51 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,7 +31,7 @@
   dependencies = [
     "numpy",
     "scipy",
-    "sympy>=1.7.1,<1.13", # <=1.13 required because of https://github.com/sympy/sympy/issues/26789 and https://github.com/sympy/sympy/issues/26791.
+    "sympy>=1.7.1",
     "pandas",
     "atpublic",
     "cython",
@@ -65,6 +65,7 @@
   "smartcard" = ["pyscard"]
   "leia" = ["smartleia"]
   "gmp" = ["gmpy2"]
+  "flint" = ["python-flint>=0.5.0"]
   "pari" = ["cysignals", "cypari2"]
   "dev" = ["mypy", "flake8", "interrogate", "pyinstrument", "black", "types-setuptools", "pydocstyle"]
   "test" = ["pytest>=7.0.0", "coverage", "pytest-cov", "pytest-sugar", "pytest-mock", "nbmake"]
diff --git a/test/ec/perf_formula.py b/test/ec/perf_formula.py
index 0a5f2a0..60b2b7d 100755
--- a/test/ec/perf_formula.py
+++ b/test/ec/perf_formula.py
@@ -1,19 +1,19 @@
 #!/usr/bin/env python
 import click
 
-from pyecsca.ec.mod import has_gmp
+from pyecsca.ec.mod import has_gmp, has_flint
 from pyecsca.ec.params import get_params
 from pyecsca.misc.cfg import TemporaryConfig
 from test.utils import Profiler
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c")), default="py")
+@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
 @click.option(
     "-m",
     "--mod",
-    type=click.Choice(("python", "gmp")),
-    default="gmp" if has_gmp else "python",
+    type=click.Choice(("python", "gmp", "flint")),
+    default="flint" if has_flint else "gmp" if has_gmp else "python",
 )
 @click.option("-o", "--operations", type=click.INT, default=5000)
 @click.option(
diff --git a/test/ec/perf_mod.py b/test/ec/perf_mod.py
index 6a83314..925360e 100755
--- a/test/ec/perf_mod.py
+++ b/test/ec/perf_mod.py
@@ -1,18 +1,18 @@
 #!/usr/bin/env python
 import click
 
-from pyecsca.ec.mod import Mod, has_gmp
+from pyecsca.ec.mod import Mod, has_gmp, has_flint
 from pyecsca.misc.cfg import TemporaryConfig
 from test.utils import Profiler
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c")), default="py")
+@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
 @click.option(
     "-m",
     "--mod",
-    type=click.Choice(("python", "gmp")),
-    default="gmp" if has_gmp else "python",
+    type=click.Choice(("python", "gmp", "flint")),
+    default="flint" if has_flint else "gmp" if has_gmp else "python",
 )
 @click.option("-o", "--operations", type=click.INT, default=100000)
 @click.option(
diff --git a/test/ec/perf_mult.py b/test/ec/perf_mult.py
index d5e6a83..50a7191 100755
--- a/test/ec/perf_mult.py
+++ b/test/ec/perf_mult.py
@@ -3,8 +3,9 @@ from typing import cast
 
 import click
 
+from pyecsca.ec.context import local, DefaultContext
 from pyecsca.ec.formula import AdditionFormula, DoublingFormula
-from pyecsca.ec.mod import has_gmp
+from pyecsca.ec.mod import has_gmp, has_flint
 from pyecsca.ec.mult import LTRMultiplier
 from pyecsca.ec.params import get_params
 from pyecsca.misc.cfg import TemporaryConfig
@@ -12,12 +13,12 @@ from test.utils import Profiler
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c")), default="py")
+@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
 @click.option(
     "-m",
     "--mod",
-    type=click.Choice(("python", "gmp")),
-    default="gmp" if has_gmp else "python",
+    type=click.Choice(("python", "gmp", "flint")),
+    default="flint" if has_flint else "gmp" if has_gmp else "python",
 )
 @click.option("-o", "--operations", type=click.INT, default=50)
 @click.option(
@@ -45,6 +46,17 @@ def main(profiler, mod, operations, directory):
                 one_point = mult.multiply(
                     0x71A55E0C1ABB3A0E069419E0F837BC195F1B9545E69FC51E53C4D48D7FEA3B1A
                 )
+        click.echo(
+            f"Profiling {operations} {p256.curve.prime.bit_length()}-bit scalar multiplication executions (with tracing)..."
+        )
+        with local(DefaultContext()):
+            one_point = p256.generator
+            with Profiler(profiler, directory, f"mult_ltr_rcb_p256_wtrace_{operations}_{mod}"):
+                for _ in range(operations):
+                    mult.init(p256, one_point)
+                    one_point = mult.multiply(
+                        0x71A55E0C1ABB3A0E069419E0F837BC195F1B9545E69FC51E53C4D48D7FEA3B1A
+                    )
 
 
 if __name__ == "__main__":
diff --git a/test/ec/test_divpoly.py b/test/ec/test_divpoly.py
index 63b3b2b..76b5f29 100644
--- a/test/ec/test_divpoly.py
+++ b/test/ec/test_divpoly.py
@@ -111,7 +111,8 @@ def test_divpoly(secp128r1):
         (4,): K(340282366762482138434845932244680310753),
         (6,): K(2),
     }
-    assert divpoly(secp128r1.curve, 4, 0).as_dict() == coeffs_0
+    d_0 = {i: K.from_sympy(d) for i, d in divpoly(secp128r1.curve, 4, 0).as_dict().items()}
+    assert d_0 == coeffs_0
     coeffs_1 = {
         (6, 1): K(4),
         (4, 1): K(340282366762482138434845932244680310723),
@@ -120,7 +121,8 @@ def test_divpoly(secp128r1):
         (1, 1): K(199419663881059632785909763469900629947),
         (0, 1): K(32040881350774765434229461361345098032),
     }
-    assert divpoly(secp128r1.curve, 4, 1).as_dict() == coeffs_1
+    d_1 = {i: K.from_sympy(d) for i, d in divpoly(secp128r1.curve, 4, 1).as_dict().items()}
+    assert d_1 == coeffs_1
     coeffs_2 = {
         (9,): K(8),
         (7,): K(340282366762482138434845932244680310639),
@@ -131,7 +133,8 @@ def test_divpoly(secp128r1):
         (1,): K(51914434605509249526780779992574428819),
         (0,): K(60581150995923875019702403440670701629),
     }
-    assert divpoly(secp128r1.curve, 4, 2).as_dict() == coeffs_2
+    d_2 = {i: K.from_sympy(d) for i, d in divpoly(secp128r1.curve, 4, 2).as_dict().items()}
+    assert d_2 == coeffs_2
 
 
 def test_mult_by_n(secp128r1):
@@ -162,11 +165,11 @@ def test_mult_by_n(secp128r1):
     }
     mx, my = mult_by_n(secp128r1.curve, 2)
     mx_num, mx_denom = mx
-    assert coeffs_mx_num == list(map(lambda x: K(int(x)), mx_num.all_coeffs()))
-    assert coeffs_mx_denom == list(map(lambda x: K(int(x)), mx_denom.all_coeffs()))
+    assert coeffs_mx_num == list(map(K.from_sympy, mx_num.all_coeffs()))
+    assert coeffs_mx_denom == list(map(K.from_sympy, mx_denom.all_coeffs()))
     my_num, my_denom = my
-    assert my_num.as_dict() == coeffs_my_num
-    assert my_denom.as_dict() == coeffs_my_denom
+    assert {i: K.from_sympy(d) for i, d in my_num.as_dict().items()} == coeffs_my_num
+    assert {i: K.from_sympy(d) for i, d in my_denom.as_dict().items()} == coeffs_my_denom
 
 
 def test_mult_by_n_large(secp128r1):
@@ -187,10 +190,10 @@ def test_mult_by_n_large(secp128r1):
             eval(key): K(val) for key, val in sage_data["my"][1].items()  # eval is OK here, skipcq: PYL-W0123
         }
 
-        assert mx[0].as_dict() == sage_data["mx"][0]
-        assert mx[1].as_dict() == sage_data["mx"][1]
-        assert my[0].as_dict() == sage_data["my"][0]
-        assert my[1].as_dict() == sage_data["my"][1]
+        assert {i: K.from_sympy(d) for i, d in mx[0].as_dict().items()} == sage_data["mx"][0]
+        assert {i: K.from_sympy(d) for i, d in mx[1].as_dict().items()} == sage_data["mx"][1]
+        assert {i: K.from_sympy(d) for i, d in my[0].as_dict().items()} == sage_data["my"][0]
+        assert {i: K.from_sympy(d) for i, d in my[1].as_dict().items()} == sage_data["my"][1]
 
 
 def test_mult_by_n_pari(secp128r1):
diff --git a/test/ec/test_formula.py b/test/ec/test_formula.py
index 24cdc7c..7b0bee1 100644
--- a/test/ec/test_formula.py
+++ b/test/ec/test_formula.py
@@ -109,16 +109,32 @@ def test_assumptions(secp128r1, mdbl):
         assert pt is not None
 
 
-def test_parameters():
-    jac_secp128r1 = get_params("secg", "secp128r1", "jacobian")
-    jac_dbl = jac_secp128r1.curve.coordinate_model.formulas["dbl-1998-hnm"]
+@pytest.mark.parametrize(
+    "formula,category,curve,coords",
+    [("dbl-1998-hnm", "secg", "secp128r1", "jacobian"),
+     ("add-2015-rcb", "secg", "secp128r1", "projective"),
+     ("dbl-1987-m-2", "other", "Curve25519", "xz"),
+     ("add-20090311-hwcd", "other", "E-222", "projective")]
+)
+def test_eval(formula, category, curve, coords):
+    params = get_params(category, curve, coords)
+    f = params.curve.coordinate_model.formulas[formula]
+
+    points_aff = [params.curve.affine_random() for _ in range(f.num_inputs)]
+    points = [point.to_model(params.curve.coordinate_model, params.curve) for point in points_aff]
+    expected = params.curve.affine_double(*points_aff) if f.shortname == "dbl" else params.curve.affine_add(*points_aff)
 
-    res = jac_dbl(
-        jac_secp128r1.curve.prime,
-        jac_secp128r1.generator,
-        **jac_secp128r1.curve.parameters,
+    res = f(
+        params.curve.prime,
+        *points,
+        **params.curve.parameters,
     )
     assert res is not None
+    try:
+        res_aff = res[0].to_affine()
+        assert expected == res_aff
+    except NotImplementedError:
+        pass
 
 
 def test_symbolic(secp128r1, dbl):
@@ -139,7 +155,7 @@ def test_symbolic(secp128r1, dbl):
         generator_val = getattr(generator_double, outer_var).x
         for inner_var in coords.variables:
             symbolic_val = symbolic_val.subs(
-                inner_var, k(getattr(secp128r1.generator, inner_var).x)
+                inner_var, int(getattr(secp128r1.generator, inner_var).x)
             )
         assert Mod(int(symbolic_val), p) == Mod(generator_val, p)
 
diff --git a/test/ec/test_mod.py b/test/ec/test_mod.py
index 1b9a83e..9a6cef0 100644
--- a/test/ec/test_mod.py
+++ b/test/ec/test_mod.py
@@ -157,7 +157,6 @@ def test_other():
     assert 5 // b == Mod(4, 7)
     assert a / 3 == Mod(4, 7)
     assert a // 3 == Mod(4, 7)
-    assert divmod(a, b) == (Mod(1, 7), Mod(2, 7))
     assert a + b == Mod(1, 7)
     assert 5 + b == Mod(1, 7)
     assert a + 3 == Mod(1, 7)
@@ -212,21 +211,3 @@ def test_symbolic():
     r = sx * a + b
     assert isinstance(r, SymbolicMod)
     assert r.n == p
-    sa = SymbolicMod(a, p)
-    sb = SymbolicMod(b, p)
-    assert sa == 3
-    assert sa.inverse() == SymbolicMod(k(9), p)
-    assert 1 / sa == SymbolicMod(k(9), p)
-    assert sa + sb == 8
-    assert 1 + sa == 4
-    assert sa - 1 == 2
-    assert 1 - sa == 11
-    assert sa + 1 == 4
-    assert -sa == 10
-    assert sa / 2 == 8
-    assert 2 / sa == 5
-    assert sa // 2 == 8
-    assert 2 // sa == 5
-    assert int(sa) == 3
-    assert sa != sb
-    assert hash(sa) is not None
diff --git a/test/sca/perf_combine.py b/test/sca/perf_combine.py
index b76acb1..bd5db59 100644
--- a/test/sca/perf_combine.py
+++ b/test/sca/perf_combine.py
@@ -16,7 +16,7 @@ from pyecsca.sca import (
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c")), default="py")
+@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
 @click.option("-o", "--operations", type=click.INT, default=100)
 @click.option(
     "-d",
diff --git a/test/sca/perf_zvp.py b/test/sca/perf_zvp.py
index 88bfe74..2175d49 100644
--- a/test/sca/perf_zvp.py
+++ b/test/sca/perf_zvp.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 import click
 
-from pyecsca.ec.mod import has_gmp
+from pyecsca.ec.mod import has_gmp, has_flint
 from pyecsca.misc.cfg import TemporaryConfig
 from pyecsca.sca.re.zvp import zvp_points, map_to_affine
 from pyecsca.ec.formula.unroll import unroll_formula
@@ -10,12 +10,12 @@ from test.utils import Profiler
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c")), default="py")
+@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
 @click.option(
     "-m",
     "--mod",
-    type=click.Choice(("python", "gmp")),
-    default="gmp" if has_gmp else "python",
+    type=click.Choice(("python", "gmp", "flint")),
+    default="flint" if has_flint else "gmp" if has_gmp else "python",
 )
 @click.option("-o", "--operations", type=click.INT, default=1)
 @click.option(
diff --git a/test/utils.py b/test/utils.py
index d4893eb..5d1e80f 100644
--- a/test/utils.py
+++ b/test/utils.py
@@ -1,19 +1,43 @@
 import pstats
 import sys
+import time
 
 from pathlib import Path
 from subprocess import run, PIPE, DEVNULL
+from typing import Union, Literal
 
 from pyinstrument import Profiler as PyProfiler
 from cProfile import Profile as cProfiler
 
 
+class RawTimer:
+    start: int
+    end: int
+    duration: float
+
+    def __enter__(self):
+        self.start = time.perf_counter_ns()
+
+    def __exit__(self, exc_type, exc_val, exc_tb):
+        self.end = time.perf_counter_ns()
+        self.duration = (self.end - self.start) / 1e9
+
+
 class Profiler:
-    def __init__(self, prof_type, output_directory, benchmark_name):
-        self._prof = PyProfiler() if prof_type == "py" else cProfiler()
-        self._prof_type = prof_type
+    def __init__(
+        self,
+        prof_type: Union[Literal["py"], Literal["c"], Literal["raw"]],
+        output_directory: str,
+        benchmark_name: str,
+    ):
+        self._prof: Union[PyProfiler, cProfiler, RawTimer] = {
+            "py": PyProfiler,
+            "c": cProfiler,
+            "raw": RawTimer,
+        }[prof_type]()
+        self._prof_type: Union[Literal["py"], Literal["c"], Literal["raw"]] = prof_type
         self._root_frame = None
-        self._state = None
+        self._state = "out"
         self._output_directory = output_directory
         self._benchmark_name = benchmark_name
 
@@ -25,7 +49,7 @@ class Profiler:
     def __exit__(self, exc_type, exc_val, exc_tb):
         self._prof.__exit__(exc_type, exc_val, exc_tb)
         if self._prof_type == "py":
-            self._root_frame = self._prof.last_session.root_frame()
+            self._root_frame = self._prof.last_session.root_frame()  # type: ignore
         self._state = "out"
         self.output()
         self.save()
@@ -62,14 +86,18 @@ class Profiler:
         if self._state != "out":
             raise ValueError
         if self._prof_type == "py":
-            print(self._prof.output_text(unicode=True, color=True))
-        else:
-            self._prof.print_stats("cumtime")
+            print(self._prof.output_text(unicode=True, color=True))  # type: ignore
+        elif self._prof_type == "c":
+            self._prof.print_stats("cumtime")  # type: ignore
+        elif self._prof_type == "raw":
+            print(f"{self._prof.duration:.4} s")  # type: ignore
 
     def get_time(self) -> float:
         if self._state != "out":
             raise ValueError
         if self._prof_type == "py":
-            return self._root_frame.time
-        else:
+            return self._root_frame.time  # type: ignore
+        elif self._prof_type == "c":
             return pstats.Stats(self._prof).total_tt  # type: ignore
+        elif self._prof_type == "raw":
+            return self._prof.duration  # type: ignore