Merge pull request #67 from J08nY/feat/cython-mod

Prepare Mod for cython

40 files changed, 1425 insertions, 1224 deletions

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index de79e4a..77da1c1 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -11,12 +11,14 @@ repos:
     -   id: check-yaml
     -   id: check-added-large-files
 -   repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v0.991
+    rev: v1.10.1
     hooks:
     -   id: mypy
         additional_dependencies:
             - "types-setuptools"
             - "numpy"
+            - "gmpy2"
+            - "python-flint"
         args: [--ignore-missing-imports, --show-error-codes, --namespace-packages, --explicit-package-bases, --check-untyped-defs]
 -   repo: https://github.com/PyCQA/flake8
     rev: 6.0.0
diff --git a/MANIFEST.in b/MANIFEST.in
index f3de4f2..9cca32a 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,4 +1,5 @@
 include README.md
+include docs/_static/logo_black_full.png
 graft pyecsca/ec/efd/
 graft pyecsca/ec/std/
 global-exclude .git*
diff --git a/docs/index.rst b/docs/index.rst
index 67fd929..34abac6 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -19,9 +19,7 @@ pyecsca [pɪɛtska]
 
 **pyecsca** aims to fill a gap in SCA tooling for Elliptic Curve Cryptography, it focuses on
 black-box implementations of ECC and presents a way to extract implementation information
-about a black-box implementation of ECC through side-channels. The main goal of **pyecsca**
-is to be able to reverse engineer the curve model, coordinate system, addition formulas, scalar
-multiplier and even finite-field implementation details.
+about a black-box implementation of ECC through side-channels.
 
 It currently provides:
 
@@ -64,12 +62,12 @@ It currently provides:
 
     .. grid-item-card::  Core
 
-        The `core <https://github.com/J08nY/pyecsca>`_ package contains the core of the
+        The `core <https://github.com/J08nY/pyecsca>`_ repository contains the core of the
         functionality, except the code generation and notebooks.
 
     .. grid-item-card::  Codegen
 
-        The `codegen <https://github.com/J08nY/pyecsca-codegen>`_ package contains
+        The `codegen <https://github.com/J08nY/pyecsca-codegen>`_ repository contains
         the code generation functionality.
 
     .. grid-item-card::  Notebook
diff --git a/docs/installation.rst b/docs/installation.rst
index 2396171..551fd5f 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -2,8 +2,64 @@
 :fas:`screwdriver-wrench;fa-fw` Installation
 ============================================
 
+**pyecsca** consists of three repositories:
+
+.. grid:: 1 1 3 3
+
+    .. grid-item-card::  Core
+
+        The `core <https://github.com/J08nY/pyecsca>`_ repository contains the core of the
+        functionality, except the code generation and notebooks.
+
+    .. grid-item-card::  Codegen
+
+        The `codegen <https://github.com/J08nY/pyecsca-codegen>`_ repository contains
+        the code generation functionality.
+
+    .. grid-item-card::  Notebook
+
+        The `notebook <https://github.com/J08nY/pyecsca-notebook>`_ repository contains
+        example notebooks that showcase functionality of the toolkit.
+
+
+Core
+====
+
+The core package can be installed either from `pypi <https://pypi.org/project/pyecsca/>`__ or from the
+`source repository <https://github.com/J08nY/pyecsca>`__. There are several extras that can be installed:
+
+- `picoscope_sdk` to enable support for PicoScope oscilloscopes using the picosdk_ package.
+- `picoscope_alt` to enable support for PicoScope oscilloscopes using the picoscope_ package.
+- `chipwhisperer` to enable support for ChipWhisperer_ targets and oscilloscopes.
+- `smartcard` to enable support for smartcard targets using the pyscard_ package.
+- `leia` to enable support for smartcard targets using the leia_ (smartleia) package.
+- `gmp` to enable arithmetic via gmpy2_ (which may or may not be faster).
+- `flint` to enable arithmetic via python-flint_ (which may or may not be faster).
+- `pari` to enable faster division polynomial computation using cypari2_.
+- `dev` to install several packages used in development.
+- `test` to install several packages used for testing.
+- `doc` to install several packages used for building documentation.
+
+You can install these extras like this:
+
+.. code-block:: shell
+
+    pip install pyecsca[smartcard,gmp]
+
+.. note::
+
+    The core repository uses git submodules, make sure to check them out after cloning with: ``git submodule update --init``.
+
+
+The core package contains data from the `Explicit-Formulas Database`_ by Daniel J. Bernstein and Tanja Lange.
+The data was partially changed, to make working with it easier. It is available on Github at `crocs-muni/efd`_.
+
+It uses `ChipWhisperer`_ as one of its targets. It also supports working with Riscure_ Inspector trace sets, which are of a proprietary format.
+
+Optionally, you can Cythonize the ``pyecsca/ec/mod`` subpackage and sometimes gain a performance benefit, YMMV.
+
 Requirements
-============
+------------
 
 .. dropdown:: General
    :open:
@@ -43,12 +99,6 @@ Requirements
          - 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.
-    The data was partially changed, to make working with it easier. It is available on Github at `crocs-muni/efd`_.
-
-    It uses `ChipWhisperer`_ as one of its targets. It also supports working with Riscure_ Inspector trace sets, which are of a proprietary format.
-
-
 .. dropdown:: Testing & Development
 
     See the Makefile for tests, performance measurement, codestyle and type checking commands.
@@ -72,6 +122,31 @@ Requirements
      - sphinx-paramlinks_
      - sphinx-design_
 
+Codegen
+=======
+
+The codegen package requires (and bundles in ``ext/libtommath`` as a git submodule) a version
+of the libtommath library. The package can be either installed from `pypi <https://pypi.org/project/pyecsca-codegen/>`__ or from the
+`source repository <https://github.com/J08nY/pyecsca-codegen>`__. Note that currently, the pypi project
+contains the built package for x86_64 Linux only. Thus, installation from source is preferable.
+
+Assuming you have ``make``, a C compiler and a C cross-compiler for ```arm-none-eabi`` you can just run:
+
+.. code-block:: shell
+
+    pip install .
+
+inside the codegen repository and it should be built and installed automatically.
+
+.. note::
+
+    The codegen repository uses git submodules, make sure to check them out after cloning with: ``git submodule update --init``.
+
+Notebooks
+=========
+
+The notebook repository is included as a submodule in the core repository.
+However, this version can get outdated during active development.
 
 .. _Numpy: https://www.numpy.org
 .. _Scipy: https://www.scipy.org
diff --git a/pyecsca/ec/curve.py b/pyecsca/ec/curve.py
index ff89819..8a3c726 100644
--- a/pyecsca/ec/curve.py
+++ b/pyecsca/ec/curve.py
@@ -12,7 +12,7 @@ from pyecsca.misc.cfg import getconfig
 
 from pyecsca.ec.coordinates import CoordinateModel, AffineCoordinateModel
 from pyecsca.ec.error import raise_unsatisified_assumption
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.model import CurveModel
 from pyecsca.ec.point import Point, InfinityPoint
 
@@ -37,12 +37,12 @@ class EllipticCurve:
 
     >>> from pyecsca.ec.coordinates import AffineCoordinateModel
     >>> affine = AffineCoordinateModel(curve.model)
-    >>> points_P = sorted(curve.affine_lift_x(Mod(5, curve.prime)), key=lambda p: int(p.y))
+    >>> points_P = sorted(curve.affine_lift_x(mod(5, curve.prime)), key=lambda p: int(p.y))
     >>> points_P  # doctest: +NORMALIZE_WHITESPACE
     [Point([x=5, y=31468013646237722594854082025316614106172411895747863909393730389177298123724] in shortw/affine),
      Point([x=5, y=84324075564118526167843364924090959423913731519542450286139900919689799730227] in shortw/affine)]
     >>> P = points_P[0]
-    >>> Q = Point(affine, x=Mod(106156966968002564385990772707119429362097710917623193504777452220576981858057, curve.prime), y=Mod(89283496902772247016522581906930535517715184283144143693965440110672128480043, curve.prime))
+    >>> Q = Point(affine, x=mod(106156966968002564385990772707119429362097710917623193504777452220576981858057, curve.prime), y=mod(89283496902772247016522581906930535517715184283144143693965440110672128480043, curve.prime))
     >>> curve.affine_add(P, Q)
     Point([x=110884201872336783252492544257507655322265785208411447156687491781308462893723, y=17851997459724035659875545393642578516937407971293368958749928013979790074156] in shortw/affine)
     >>> curve.affine_multiply(P, 10)
@@ -93,9 +93,10 @@ class EllipticCurve:
             if isinstance(value, Mod):
                 if value.n != prime:
                     raise ValueError(f"Parameter {name} has wrong modulus.")
+                val = value
             else:
-                value = Mod(value, prime)
-            self.parameters[name] = value
+                val = mod(value, prime)
+            self.parameters[name] = val
         self.neutral = neutral
         self.__validate_coord_assumptions()
 
@@ -147,9 +148,9 @@ class EllipticCurve:
         for line in formulas:
             exec(compile(line, "", mode="exec"), None, locls)  # exec is OK here, skipcq: PYL-W0122
         if not isinstance(locls["x"], Mod):
-            locls["x"] = Mod(locls["x"], self.prime)
+            locls["x"] = mod(locls["x"], self.prime)
         if not isinstance(locls["y"], Mod):
-            locls["y"] = Mod(locls["y"], self.prime)
+            locls["y"] = mod(locls["y"], self.prime)
         return Point(AffineCoordinateModel(self.model), x=locls["x"], y=locls["y"])
 
     def affine_add(self, one: Point, other: Point) -> Point:
@@ -234,9 +235,9 @@ class EllipticCurve:
         for line in self.model.base_neutral:
             exec(compile(line, "", mode="exec"), None, locls)  # exec is OK here, skipcq: PYL-W0122
         if not isinstance(locls["x"], Mod):
-            locls["x"] = Mod(locls["x"], self.prime)
+            locls["x"] = mod(locls["x"], self.prime)
         if not isinstance(locls["y"], Mod):
-            locls["y"] = Mod(locls["y"], self.prime)
+            locls["y"] = mod(locls["y"], self.prime)
         return Point(AffineCoordinateModel(self.model), x=locls["x"], y=locls["y"])
 
     @property
@@ -314,7 +315,7 @@ class EllipticCurve:
                 raise ValueError("Encoded point has bad length")
             coords = {}
             for var in sorted(self.coordinate_model.variables):
-                coords[var] = Mod(int.from_bytes(data[:coord_len], "big"), self.prime)
+                coords[var] = mod(int.from_bytes(data[:coord_len], "big"), self.prime)
                 data = data[coord_len:]
             return Point(self.coordinate_model, **coords)
         elif encoded[0] in (0x02, 0x03):
@@ -322,7 +323,7 @@ class EllipticCurve:
                 data = encoded[1:]
                 if len(data) != coord_len:
                     raise ValueError("Encoded point has bad length")
-                x = Mod(int.from_bytes(data, "big"), self.prime)
+                x = mod(int.from_bytes(data, "big"), self.prime)
                 loc = {**self.parameters, "x": x}
                 rhs = eval(compile(self.model.ysquared, "", mode="eval"), loc)  # eval is OK here, skipcq: PYL-W0123
                 if not rhs.is_residue():
diff --git a/pyecsca/ec/divpoly.py b/pyecsca/ec/divpoly.py
index 722ec97..912fb06 100644
--- a/pyecsca/ec/divpoly.py
+++ b/pyecsca/ec/divpoly.py
@@ -8,7 +8,7 @@ from sympy import symbols, FF, Poly
 import networkx as nx
 
 from pyecsca.ec.curve import EllipticCurve
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.model import ShortWeierstrassModel
 
 has_pari = False
@@ -90,9 +90,9 @@ def a_invariants(curve: EllipticCurve) -> Tuple[Mod, ...]:
     :return: A tuple of 5 a-invariants (a1, a2, a3, a4, a6).
     """
     if isinstance(curve.model, ShortWeierstrassModel):
-        a1 = Mod(0, curve.prime)
-        a2 = Mod(0, curve.prime)
-        a3 = Mod(0, curve.prime)
+        a1 = mod(0, curve.prime)
+        a2 = mod(0, curve.prime)
+        a3 = mod(0, curve.prime)
         a4 = curve.parameters["a"]
         a6 = curve.parameters["b"]
         return a1, a2, a3, a4, a6
diff --git a/pyecsca/ec/formula/base.py b/pyecsca/ec/formula/base.py
index fee39aa..04d6c19 100644
--- a/pyecsca/ec/formula/base.py
+++ b/pyecsca/ec/formula/base.py
@@ -13,7 +13,7 @@ from sympy import FF, symbols, Poly
 from pyecsca.ec.context import ResultAction
 from pyecsca.ec import context
 from pyecsca.ec.error import UnsatisfiedAssumptionError, raise_unsatisified_assumption
-from pyecsca.ec.mod import Mod, SymbolicMod
+from pyecsca.ec.mod import Mod, mod, SymbolicMod
 from pyecsca.ec.op import CodeOp, OpType
 from pyecsca.misc.cfg import getconfig
 from pyecsca.misc.cache import sympify
@@ -173,7 +173,7 @@ class Formula(ABC):
         used_symbols = sorted(expr.free_symbols)
         used_params = []
         for symbol in used_symbols:
-            if (value := params.get(str(symbol), None)) is not None:
+            if (value := params.get(symbol.name, None)) is not None:
                 used_params.append(value)
                 if isinstance(value, SymbolicMod):
                     expr = expr.xreplace({symbol: value.x})
@@ -191,7 +191,7 @@ class Formula(ABC):
                 domain = FF(field)
                 numerator, denominator = expr.as_numer_denom()
                 val = int(domain.from_sympy(numerator) / domain.from_sympy(denominator))
-                params[lhs] = Mod(val, field)
+                params[lhs] = mod(val, field)
             _assumption_cache[cache_key] = params[lhs]
         return True
 
@@ -200,14 +200,14 @@ class Formula(ABC):
         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 (value := params.get(symbol.name, 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:
+        if len(remaining) > 1 or (param := remaining[0].name) not in self.parameters:
             raise ValueError(
                 f"This formula couldn't be executed due to an unsupported assumption ({assumption_string})."
             )
@@ -216,7 +216,7 @@ class Formula(ABC):
         poly = Poly(numerator, symbols(param), domain=domain)
         roots = poly.ground_roots()
         for root in roots:
-            params[param] = Mod(int(domain.from_sympy(root)), field)
+            params[param] = mod(int(domain.from_sympy(root)), field)
             return
         raise UnsatisfiedAssumptionError(
             f"Unsatisfied assumption in the formula ({assumption_string}).\n"
@@ -274,7 +274,7 @@ class Formula(ABC):
                         f"Bad stuff happened in op {op}, floats will pollute the results."
                     )
                 if not isinstance(op_result, Mod):
-                    op_result = Mod(op_result, field)
+                    op_result = mod(op_result, field)
                 if context.current is not None:
                     action.add_operation(op, op_result)
                 params[op.result] = op_result
diff --git a/pyecsca/ec/formula/switch_sign.py b/pyecsca/ec/formula/switch_sign.py
index b5bc146..e6397a0 100644
--- a/pyecsca/ec/formula/switch_sign.py
+++ b/pyecsca/ec/formula/switch_sign.py
@@ -7,7 +7,7 @@ from pyecsca.ec.op import OpType, CodeOp
 from pyecsca.ec.formula.base import Formula
 from pyecsca.ec.formula.graph import FormulaGraph, ConstantNode, CodeOpNode, CodeFormula
 from pyecsca.ec.point import Point
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 
 
 @public
@@ -59,7 +59,7 @@ def sign_test(output_signs: Dict[str, int], coordinate_model: Any):
             out_var = out[: out.index(ind)]
             if not out_var.isalpha():
                 continue
-            point_dict[out_var] = Mod(sign, p)
+            point_dict[out_var] = mod(sign, p)
         point = Point(coordinate_model, **point_dict)
         try:
             apoint = point.to_affine()
@@ -68,7 +68,7 @@ def sign_test(output_signs: Dict[str, int], coordinate_model: Any):
             if scale is None:
                 raise BadSignSwitch
             apoint = scale(p, point)[0]
-        if set(apoint.coords.values()) != {Mod(1, p)}:
+        if set(apoint.coords.values()) != {mod(1, p)}:
             raise BadSignSwitch
 
 
diff --git a/pyecsca/ec/mod.py b/pyecsca/ec/mod.py
deleted file mode 100644
index ecdf156..0000000
--- a/pyecsca/ec/mod.py
+++ /dev/null
@@ -1,963 +0,0 @@
-"""
-Provides several implementations of an element of ℤₙ.
-
-The base class :py:class:`Mod` dynamically
-dispatches to the implementation chosen by the runtime configuration of the library
-(see :py:class:`pyecsca.misc.cfg.Config`). A Python integer based implementation is available under
-:py:class:`RawMod`. A symbolic implementation based on sympy is available under :py:class:`SymbolicMod`. If
-`gmpy2` is installed, a GMP based implementation is available under :py:class:`GMPMod`. If `python-flint` is
-installed, a flint based implementation is available under :py:class:`FlintMod`.
-"""
-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
-
-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
-
-has_gmp = False
-try:
-    import gmpy2
-
-    has_gmp = True
-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: int, b: int) -> int:
-    """Euclid's greatest common denominator algorithm."""
-    if abs(a) < abs(b):
-        return gcd(b, a)
-
-    while abs(b) > 0:
-        _, r = divmod(a, b)
-        a, b = b, r
-
-    return a
-
-
-@public
-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)
-        return y, x, d
-
-    if abs(b) == 0:
-        return 1, 0, a
-
-    x1, x2, y1, y2 = 0, 1, 1, 0
-    while abs(b) > 0:
-        q, r = divmod(a, b)
-        x = x2 - q * x1
-        y = y2 - q * y1
-        a, b, x2, x1, y2, y1 = b, r, x1, x, y1, y
-
-    return x2, y2, a
-
-
-@public
-def jacobi(x: int, n: int) -> int:
-    """Jacobi symbol."""
-    if n <= 0:
-        raise ValueError("'n' must be a positive integer.")
-    if n % 2 == 0:
-        raise ValueError("'n' must be odd.")
-    x %= n
-    r = 1
-    while x != 0:
-        while x % 2 == 0:
-            x //= 2
-            nm8 = n % 8
-            if nm8 in (3, 5):
-                r = -r
-        x, n = n, x
-        if x % 4 == 3 and n % 4 == 3:
-            r = -r
-        x %= n
-    return r if n == 1 else 0
-
-
-@public
-@lru_cache
-def miller_rabin(n: int, rounds: int = 50) -> bool:
-    """Miller-Rabin probabilistic primality test."""
-    if n in (2, 3):
-        return True
-
-    if n % 2 == 0:
-        return False
-
-    r, s = 0, n - 1
-    while s % 2 == 0:
-        r += 1
-        s //= 2
-    for _ in range(rounds):
-        a = random.randrange(2, n - 1)
-        x = pow(a, s, n)
-        if x in (1, n - 1):
-            continue
-        for _ in range(r - 1):
-            x = pow(x, 2, n)
-            if x == n - 1:
-                break
-        else:
-            return False
-    return True
-
-
-def _check(func):
-    @wraps(func)
-    def method(self, other):
-        if self.__class__ is not type(other):
-            other = self.__class__(other, self.n)
-        elif self.n != other.n:
-            raise ValueError
-        return func(self, other)
-
-    return method
-
-
-@public
-class RandomModAction(ResultAction):
-    """A random sampling from Z_n."""
-
-    order: int
-
-    def __init__(self, order: int):
-        super().__init__()
-        self.order = order
-
-    def __repr__(self):
-        return f"{self.__class__.__name__}({self.order:x})"
-
-
-_mod_classes: Dict[str, Type] = {}
-_mod_order = ["gmp", "flint", "python"]
-
-
-@public
-class Mod:
-    """
-    An element x of ℤₙ.
-
-    .. note::
-        This class dispatches to one of :py:class:`RawMod`, :py:class:`GMPMod` or :py:class:`FlintMod`
-        based on what packages are installed and what implementation is configured (see
-        :py:mod:`pyecsca.misc.cfg`).
-
-    Has all the usual special methods that upcast integers automatically:
-
-    >>> a = Mod(3, 5)
-    >>> b = Mod(2, 5)
-    >>> a + b
-    0
-    >>> a * 2
-    1
-    >>> a == 3
-    True
-    >>> a == -2
-    True
-    >>> -a
-    2
-
-    Plus some additional useful things:
-
-    >>> a.inverse()
-    2
-    >>> a.is_residue()
-    False
-    >>> (a**2).is_residue()
-    True
-    >>> (a**2).sqrt() in (a, -a)
-    True
-    """
-
-    x: Any
-    n: Any
-    __slots__ = ("x", "n")
-
-    def __new__(cls, *args, **kwargs) -> "Mod":
-        if cls != Mod:
-            return cls.__new__(cls, *args, **kwargs)
-        if not _mod_classes:
-            raise ValueError("Cannot find any working Mod class.")
-        selected_class = getconfig().ec.mod_implementation
-        if selected_class not in _mod_classes:
-            # Fallback to something
-            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
-        )
-
-    @_check
-    def __add__(self, other) -> "Mod":
-        return self.__class__((self.x + other.x) % self.n, self.n)
-
-    @_check
-    def __radd__(self, other) -> "Mod":
-        return self + other
-
-    @_check
-    def __sub__(self, other) -> "Mod":
-        return self.__class__((self.x - other.x) % self.n, self.n)
-
-    @_check
-    def __rsub__(self, other) -> "Mod":
-        return -self + other
-
-    def __neg__(self) -> "Mod":
-        return self.__class__(self.n - self.x, self.n)
-
-    def bit_length(self):
-        """
-        Compute the bit length of this element (in its positive integer representation).
-
-        :return: The bit-length.
-        """
-        raise NotImplementedError
-
-    def inverse(self) -> "Mod":
-        """
-        Invert the element.
-
-        :return: The inverse.
-        :raises: :py:class:`NonInvertibleError` if the element is not invertible.
-        """
-        raise NotImplementedError
-
-    def __invert__(self) -> "Mod":
-        return self.inverse()
-
-    def is_residue(self) -> bool:
-        """Whether this element is a quadratic residue (only implemented for prime modulus)."""
-        raise NotImplementedError
-
-    def sqrt(self) -> "Mod":
-        """
-        Compute the modular square root of this element (only implemented for prime modulus).
-
-        Uses the `Tonelli-Shanks <https://en.wikipedia.org/wiki/Tonelli–Shanks_algorithm>`_ algorithm.
-        """
-        raise NotImplementedError
-
-    @_check
-    def __mul__(self, other) -> "Mod":
-        return self.__class__((self.x * other.x) % self.n, self.n)
-
-    @_check
-    def __rmul__(self, other) -> "Mod":
-        return self * other
-
-    @_check
-    def __truediv__(self, other) -> "Mod":
-        return self * ~other
-
-    @_check
-    def __rtruediv__(self, other) -> "Mod":
-        return ~self * other
-
-    @_check
-    def __floordiv__(self, other) -> "Mod":
-        return self * ~other
-
-    @_check
-    def __rfloordiv__(self, other) -> "Mod":
-        return ~self * other
-
-    def __bytes__(self) -> bytes:
-        raise NotImplementedError
-
-    def __int__(self) -> int:
-        raise NotImplementedError
-
-    @classmethod
-    def random(cls, n: int) -> "Mod":
-        """
-        Generate a random :py:class:`Mod` in ℤₙ.
-
-        :param n: The order.
-        :return: The random :py:class:`Mod`.
-        """
-        with RandomModAction(n) as action:
-            return action.exit(cls(secrets.randbelow(n), n))
-
-    def __pow__(self, n) -> "Mod":
-        return NotImplemented
-
-    def __str__(self):
-        return str(self.x)
-
-
-@public
-class RawMod(Mod):
-    """An element x of ℤₙ (implemented using Python integers)."""
-
-    x: int
-    n: int
-    __slots__ = ("x", "n")
-
-    def __new__(cls, *args, **kwargs):
-        return object.__new__(cls)
-
-    def __init__(self, x: int, n: int):
-        self.x = x % n
-        self.n = n
-
-    def bit_length(self):
-        return self.x.bit_length()
-
-    def inverse(self) -> "RawMod":
-        if self.x == 0:
-            raise_non_invertible()
-        x, _, d = extgcd(self.x, self.n)
-        if d != 1:
-            raise_non_invertible()
-        return RawMod(x, self.n)
-
-    def is_residue(self):
-        if not miller_rabin(self.n):
-            raise NotImplementedError
-        if self.x == 0:
-            return True
-        if self.n == 2:
-            return self.x in (0, 1)
-        legendre_symbol = jacobi(self.x, self.n)
-        return legendre_symbol == 1
-
-    def sqrt(self) -> "RawMod":
-        if not miller_rabin(self.n):
-            raise NotImplementedError
-        if self.x == 0:
-            return RawMod(0, self.n)
-        if not self.is_residue():
-            raise_non_residue()
-        if self.n % 4 == 3:
-            return self ** int((self.n + 1) // 4)
-        q = self.n - 1
-        s = 0
-        while q % 2 == 0:
-            q //= 2
-            s += 1
-
-        z = 2
-        while RawMod(z, self.n).is_residue():
-            z += 1
-
-        m = s
-        c = RawMod(z, self.n) ** q
-        t = self**q
-        r_exp = (q + 1) // 2
-        r = self**r_exp
-
-        while t != 1:
-            i = 1
-            while not (t ** (2**i)) == 1:
-                i += 1
-            two_exp = m - (i + 1)
-            b = c ** int(RawMod(2, self.n) ** two_exp)
-            m = int(RawMod(i, self.n))
-            c = b**2
-            t *= c
-            r *= b
-        return r
-
-    def __bytes__(self):
-        return self.x.to_bytes((self.n.bit_length() + 7) // 8, byteorder="big")
-
-    def __int__(self):
-        return self.x
-
-    def __eq__(self, other):
-        if type(other) is int:
-            return self.x == (other % self.n)
-        if type(other) is not RawMod:
-            return False
-        return self.x == other.x and self.n == other.n
-
-    def __ne__(self, other):
-        return not self == other
-
-    def __repr__(self):
-        return str(self.x)
-
-    def __hash__(self):
-        return hash(("RawMod", self.x, self.n))
-
-    def __pow__(self, n) -> "RawMod":
-        if type(n) is not int:
-            raise TypeError
-        if n == 0:
-            return RawMod(1, self.n)
-        if n < 0:
-            return self.inverse() ** (-n)
-        if n == 1:
-            return RawMod(self.x, self.n)
-
-        return RawMod(pow(self.x, n, self.n), self.n)
-
-
-_mod_classes["python"] = RawMod
-
-
-@public
-class Undefined(Mod):
-    """A special undefined element."""
-
-    __slots__ = ("x", "n")
-
-    def __new__(cls, *args, **kwargs):
-        return object.__new__(cls)
-
-    def __init__(self):
-        self.x = None
-        self.n = None
-
-    def __add__(self, other):
-        return NotImplemented
-
-    def __radd__(self, other):
-        return NotImplemented
-
-    def __sub__(self, other):
-        return NotImplemented
-
-    def __rsub__(self, other):
-        return NotImplemented
-
-    def __neg__(self):
-        raise NotImplementedError
-
-    def bit_length(self):
-        raise NotImplementedError
-
-    def inverse(self):
-        raise NotImplementedError
-
-    def sqrt(self):
-        raise NotImplementedError
-
-    def is_residue(self):
-        raise NotImplementedError
-
-    def __invert__(self):
-        raise NotImplementedError
-
-    def __mul__(self, other):
-        return NotImplemented
-
-    def __rmul__(self, other):
-        return NotImplemented
-
-    def __truediv__(self, other):
-        return NotImplemented
-
-    def __rtruediv__(self, other):
-        return NotImplemented
-
-    def __floordiv__(self, other):
-        return NotImplemented
-
-    def __rfloordiv__(self, other):
-        return NotImplemented
-
-    def __bytes__(self):
-        raise NotImplementedError
-
-    def __int__(self):
-        raise NotImplementedError
-
-    def __eq__(self, other):
-        return False
-
-    def __ne__(self, other):
-        return False
-
-    def __repr__(self):
-        return "Undefined"
-
-    def __hash__(self):
-        return hash("Undefined") + 1
-
-    def __pow__(self, n):
-        return NotImplemented
-
-
-@public
-class SymbolicMod(Mod):
-    """A symbolic element x of ℤₙ (implemented using sympy)."""
-
-    x: Expr
-    n: int
-    __slots__ = ("x", "n")
-
-    def __new__(cls, *args, **kwargs):
-        return object.__new__(cls)
-
-    def __init__(self, x: Expr, n: int):
-        self.x = x
-        self.n = n
-
-    @_check
-    def __add__(self, other) -> "SymbolicMod":
-        return self.__class__((self.x + other.x), self.n)
-
-    @_check
-    def __radd__(self, other) -> "SymbolicMod":
-        return self + other
-
-    @_check
-    def __sub__(self, other) -> "SymbolicMod":
-        return self.__class__((self.x - other.x), self.n)
-
-    @_check
-    def __rsub__(self, other) -> "SymbolicMod":
-        return -self + other
-
-    def __neg__(self) -> "SymbolicMod":
-        return self.__class__(-self.x, self.n)
-
-    def bit_length(self):
-        raise NotImplementedError
-
-    def inverse(self) -> "SymbolicMod":
-        return self.__class__(self.x ** (-1), self.n)
-
-    def sqrt(self) -> "SymbolicMod":
-        raise NotImplementedError
-
-    def is_residue(self):
-        raise NotImplementedError
-
-    def __invert__(self) -> "SymbolicMod":
-        return self.inverse()
-
-    @_check
-    def __mul__(self, other) -> "SymbolicMod":
-        return self.__class__(self.x * other.x, self.n)
-
-    @_check
-    def __rmul__(self, other) -> "SymbolicMod":
-        return self * other
-
-    @_check
-    def __truediv__(self, other) -> "SymbolicMod":
-        return self * ~other
-
-    @_check
-    def __rtruediv__(self, other) -> "SymbolicMod":
-        return ~self * other
-
-    @_check
-    def __floordiv__(self, other) -> "SymbolicMod":
-        return self * ~other
-
-    @_check
-    def __rfloordiv__(self, other) -> "SymbolicMod":
-        return ~self * other
-
-    def __bytes__(self):
-        return int(self.x).to_bytes((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 % self.n
-        if type(other) is not SymbolicMod:
-            return False
-        return self.x == other.x and self.n == other.n
-
-    def __ne__(self, other):
-        return not self == other
-
-    def __repr__(self):
-        return str(self.x)
-
-    def __hash__(self):
-        return hash(("SymbolicMod", self.x, self.n))
-
-    def __pow__(self, n) -> "SymbolicMod":
-        return self.__class__(pow(self.x, n), self.n)
-
-
-_mod_classes["symbolic"] = SymbolicMod
-
-if has_gmp:
-
-    @lru_cache
-    def _gmpy_is_prime(x) -> bool:
-        return gmpy2.is_prime(x)
-
-    @public
-    class GMPMod(Mod):
-        """An element x of ℤₙ. Implemented by GMP."""
-
-        x: gmpy2.mpz
-        n: gmpy2.mpz
-        __slots__ = ("x", "n")
-
-        def __new__(cls, *args, **kwargs):
-            return object.__new__(cls)
-
-        def __init__(
-            self,
-            x: Union[int, gmpy2.mpz],
-            n: Union[int, gmpy2.mpz],
-            ensure: bool = True,
-        ):
-            if ensure:
-                self.n = gmpy2.mpz(n)
-                self.x = gmpy2.mpz(x % self.n)
-            else:
-                self.n = n
-                self.x = x
-
-        def bit_length(self):
-            return self.x.bit_length()
-
-        def inverse(self) -> "GMPMod":
-            if self.x == 0:
-                raise_non_invertible()
-            if self.x == 1:
-                return GMPMod(gmpy2.mpz(1), self.n, ensure=False)
-            try:
-                res = gmpy2.invert(self.x, self.n)
-            except ZeroDivisionError:
-                raise_non_invertible()
-                res = gmpy2.mpz(0)
-            return GMPMod(res, self.n, ensure=False)
-
-        def is_residue(self) -> bool:
-            if not _gmpy_is_prime(self.n):
-                raise NotImplementedError
-            if self.x == 0:
-                return True
-            if self.n == 2:
-                return self.x in (0, 1)
-            return gmpy2.legendre(self.x, self.n) == 1
-
-        def sqrt(self) -> "GMPMod":
-            if not _gmpy_is_prime(self.n):
-                raise NotImplementedError
-            if self.x == 0:
-                return GMPMod(gmpy2.mpz(0), self.n, ensure=False)
-            if not self.is_residue():
-                raise_non_residue()
-            if self.n % 4 == 3:
-                return self ** int((self.n + 1) // 4)
-            q = self.n - 1
-            s = 0
-            while q % 2 == 0:
-                q //= 2
-                s += 1
-
-            z = gmpy2.mpz(2)
-            while GMPMod(z, self.n, ensure=False).is_residue():
-                z += 1
-
-            m = s
-            c = GMPMod(z, self.n, 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(GMPMod(gmpy2.mpz(2), self.n, ensure=False) ** two_exp)
-                m = int(GMPMod(gmpy2.mpz(i), self.n, ensure=False))
-                c = b**2
-                t *= c
-                r *= b
-            return r
-
-        @_check
-        def __add__(self, other) -> "GMPMod":
-            return GMPMod((self.x + other.x) % self.n, self.n, ensure=False)
-
-        @_check
-        def __sub__(self, other) -> "GMPMod":
-            return GMPMod((self.x - other.x) % self.n, self.n, ensure=False)
-
-        def __neg__(self) -> "GMPMod":
-            return GMPMod(self.n - self.x, self.n, ensure=False)
-
-        @_check
-        def __mul__(self, other) -> "GMPMod":
-            return GMPMod((self.x * other.x) % self.n, self.n, ensure=False)
-
-        def __bytes__(self):
-            return int(self.x).to_bytes((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 == (gmpy2.mpz(other) % self.n)
-            if type(other) is not GMPMod:
-                return False
-            return self.x == other.x and self.n == other.n
-
-        def __ne__(self, other):
-            return not self == other
-
-        def __repr__(self):
-            return str(int(self.x))
-
-        def __hash__(self):
-            return hash(("GMPMod", self.x, self.n))
-
-        def __pow__(self, n) -> "GMPMod":
-            if type(n) not in (int, gmpy2.mpz):
-                raise TypeError
-            if n == 0:
-                return GMPMod(gmpy2.mpz(1), self.n, ensure=False)
-            if n < 0:
-                return self.inverse() ** (-n)
-            if n == 1:
-                return GMPMod(self.x, self.n, ensure=False)
-            return GMPMod(
-                gmpy2.powmod(self.x, gmpy2.mpz(n), self.n), self.n, ensure=False
-            )
-
-    _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 flint."""
-
-        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/mod/__init__.py b/pyecsca/ec/mod/__init__.py
new file mode 100644
index 0000000..a5b1521
--- /dev/null
+++ b/pyecsca/ec/mod/__init__.py
@@ -0,0 +1,16 @@
+"""
+Provides several implementations of an element of ℤₙ.
+
+The base class :py:class:`Mod` dynamically
+dispatches to the implementation chosen by the runtime configuration of the library
+(see :py:class:`pyecsca.misc.cfg.Config`). A Python integer based implementation is available under
+:py:class:`RawMod`. A symbolic implementation based on sympy is available under :py:class:`SymbolicMod`. If
+`gmpy2` is installed, a GMP based implementation is available under :py:class:`GMPMod`. If `python-flint` is
+installed, a flint based implementation is available under :py:class:`FlintMod`.
+"""
+
+from .base import *
+from .raw import *
+from .symbolic import *
+from .gmp import *
+from .flint import *
diff --git a/pyecsca/ec/mod/base.py b/pyecsca/ec/mod/base.py
new file mode 100644
index 0000000..67e2b97
--- /dev/null
+++ b/pyecsca/ec/mod/base.py
@@ -0,0 +1,367 @@
+import random
+import secrets
+from functools import lru_cache, wraps
+
+from public import public
+from typing import Tuple, Any, Dict, Type
+
+from pyecsca.ec.context import ResultAction
+from pyecsca.misc.cfg import getconfig
+
+
+@public
+def gcd(a: int, b: int) -> int:
+    """Euclid's greatest common denominator algorithm."""
+    if abs(a) < abs(b):
+        return gcd(b, a)
+
+    while abs(b) > 0:
+        _, r = divmod(a, b)
+        a, b = b, r
+
+    return a
+
+
+@public
+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)
+        return y, x, d
+
+    if abs(b) == 0:
+        return 1, 0, a
+
+    x1, x2, y1, y2 = 0, 1, 1, 0
+    while abs(b) > 0:
+        q, r = divmod(a, b)
+        x = x2 - q * x1
+        y = y2 - q * y1
+        a, b, x2, x1, y2, y1 = b, r, x1, x, y1, y
+
+    return x2, y2, a
+
+
+@public
+def jacobi(x: int, n: int) -> int:
+    """Jacobi symbol."""
+    if n <= 0:
+        raise ValueError("'n' must be a positive integer.")
+    if n % 2 == 0:
+        raise ValueError("'n' must be odd.")
+    x %= n
+    r = 1
+    while x != 0:
+        while x % 2 == 0:
+            x //= 2
+            nm8 = n % 8
+            if nm8 in (3, 5):
+                r = -r
+        x, n = n, x
+        if x % 4 == 3 and n % 4 == 3:
+            r = -r
+        x %= n
+    return r if n == 1 else 0
+
+
+@public
+@lru_cache
+def miller_rabin(n: int, rounds: int = 50) -> bool:
+    """Miller-Rabin probabilistic primality test."""
+    if n in (2, 3):
+        return True
+
+    if n % 2 == 0:
+        return False
+
+    r, s = 0, n - 1
+    while s % 2 == 0:
+        r += 1
+        s //= 2
+    for _ in range(rounds):
+        a = random.randrange(2, n - 1)
+        x = pow(a, s, n)
+        if x in (1, n - 1):
+            continue
+        for _ in range(r - 1):
+            x = pow(x, 2, n)
+            if x == n - 1:
+                break
+        else:
+            return False
+    return True
+
+
+def _check(func):
+    @wraps(func)
+    def method(self, other):
+        if self.__class__ is not type(other):
+            other = self.__class__(other, self.n)
+        elif self.n != other.n:
+            raise ValueError
+        return func(self, other)
+
+    return method
+
+
+@public
+class RandomModAction(ResultAction):
+    """A random sampling from Z_n."""
+
+    order: int
+
+    def __init__(self, order: int):
+        super().__init__()
+        self.order = order
+
+    def __repr__(self):
+        return f"{self.__class__.__name__}({self.order:x})"
+
+
+_mod_classes: Dict[str, Type] = {}
+_mod_order = ["gmp", "flint", "python"]
+
+
+@public
+class Mod:
+    """
+    An element x of ℤₙ.
+
+    .. note::
+        This class dispatches to one of :py:class:`RawMod`, :py:class:`GMPMod` or :py:class:`FlintMod`
+        based on what packages are installed and what implementation is configured (see
+        :py:mod:`pyecsca.misc.cfg`).
+
+    Has all the usual special methods that upcast integers automatically:
+
+    >>> a = mod(3, 5)
+    >>> b = mod(2, 5)
+    >>> a + b
+    0
+    >>> a * 2
+    1
+    >>> a == 3
+    True
+    >>> a == -2
+    True
+    >>> -a
+    2
+
+    Plus some additional useful things:
+
+    >>> a.inverse()
+    2
+    >>> a.is_residue()
+    False
+    >>> (a**2).is_residue()
+    True
+    >>> (a**2).sqrt() in (a, -a)
+    True
+    """
+
+    x: Any
+    n: Any
+    __slots__ = ("x", "n")
+
+    def __init__(self, x, n):
+        raise TypeError("Abstract")
+
+    @_check
+    def __add__(self, other) -> "Mod":
+        return self.__class__((self.x + other.x) % self.n, self.n)
+
+    @_check
+    def __radd__(self, other) -> "Mod":
+        return self + other
+
+    @_check
+    def __sub__(self, other) -> "Mod":
+        return self.__class__((self.x - other.x) % self.n, self.n)
+
+    @_check
+    def __rsub__(self, other) -> "Mod":
+        return -self + other
+
+    def __neg__(self) -> "Mod":
+        return self.__class__(self.n - self.x, self.n)
+
+    def bit_length(self):
+        """
+        Compute the bit length of this element (in its positive integer representation).
+
+        :return: The bit-length.
+        """
+        raise NotImplementedError
+
+    def inverse(self) -> "Mod":
+        """
+        Invert the element.
+
+        :return: The inverse.
+        :raises: :py:class:`NonInvertibleError` if the element is not invertible.
+        """
+        raise NotImplementedError
+
+    def __invert__(self) -> "Mod":
+        return self.inverse()
+
+    def is_residue(self) -> bool:
+        """Whether this element is a quadratic residue (only implemented for prime modulus)."""
+        raise NotImplementedError
+
+    def sqrt(self) -> "Mod":
+        """
+        Compute the modular square root of this element (only implemented for prime modulus).
+
+        Uses the `Tonelli-Shanks <https://en.wikipedia.org/wiki/Tonelli–Shanks_algorithm>`_ algorithm.
+        """
+        raise NotImplementedError
+
+    @_check
+    def __mul__(self, other) -> "Mod":
+        return self.__class__((self.x * other.x) % self.n, self.n)
+
+    @_check
+    def __rmul__(self, other) -> "Mod":
+        return self * other
+
+    @_check
+    def __truediv__(self, other) -> "Mod":
+        return self * ~other
+
+    @_check
+    def __rtruediv__(self, other) -> "Mod":
+        return ~self * other
+
+    @_check
+    def __floordiv__(self, other) -> "Mod":
+        return self * ~other
+
+    @_check
+    def __rfloordiv__(self, other) -> "Mod":
+        return ~self * other
+
+    def __bytes__(self) -> bytes:
+        raise NotImplementedError
+
+    def __int__(self) -> int:
+        raise NotImplementedError
+
+    @classmethod
+    def random(cls, n: int) -> "Mod":
+        """
+        Generate a random :py:class:`Mod` in ℤₙ.
+
+        :param n: The order.
+        :return: The random :py:class:`Mod`.
+        """
+        with RandomModAction(n) as action:
+            return action.exit(mod(secrets.randbelow(n), n))
+
+    def __pow__(self, n, _=None) -> "Mod":
+        return NotImplemented
+
+    def __str__(self):
+        return str(self.x)
+
+
+@public
+class Undefined(Mod):
+    """A special undefined element."""
+
+    __slots__ = ("x", "n")
+
+    def __init__(self):
+        self.x = None
+        self.n = None
+
+    def __add__(self, other):
+        return NotImplemented
+
+    def __radd__(self, other):
+        return NotImplemented
+
+    def __sub__(self, other):
+        return NotImplemented
+
+    def __rsub__(self, other):
+        return NotImplemented
+
+    def __neg__(self):
+        raise NotImplementedError
+
+    def bit_length(self):
+        raise NotImplementedError
+
+    def inverse(self):
+        raise NotImplementedError
+
+    def sqrt(self):
+        raise NotImplementedError
+
+    def is_residue(self):
+        raise NotImplementedError
+
+    def __invert__(self):
+        raise NotImplementedError
+
+    def __mul__(self, other):
+        return NotImplemented
+
+    def __rmul__(self, other):
+        return NotImplemented
+
+    def __truediv__(self, other):
+        return NotImplemented
+
+    def __rtruediv__(self, other):
+        return NotImplemented
+
+    def __floordiv__(self, other):
+        return NotImplemented
+
+    def __rfloordiv__(self, other):
+        return NotImplemented
+
+    def __bytes__(self):
+        raise NotImplementedError
+
+    def __int__(self):
+        raise NotImplementedError
+
+    def __eq__(self, other):
+        return False
+
+    def __ne__(self, other):
+        return False
+
+    def __repr__(self):
+        return "Undefined"
+
+    def __hash__(self):
+        return hash("Undefined") + 1
+
+    def __pow__(self, n, _=None):
+        return NotImplemented
+
+
+@public
+def mod(x: int, n: int) -> Mod:
+    """
+    Construct a :py:class:`Mod`.
+
+    :param x: The value.
+    :param n: The modulus.
+    :return: A selected Mod implementation object.
+    :raises: ValueError in case a working Mod implementation cannot be found.
+    """
+    if not _mod_classes:
+        raise ValueError("Cannot find any working Mod class.")
+    selected_class = getconfig().ec.mod_implementation
+    if selected_class not in _mod_classes:
+        # Fallback to something
+        for fallback in _mod_order:
+            if fallback in _mod_classes:
+                selected_class = fallback
+                break
+    return _mod_classes[selected_class](x, n)
diff --git a/pyecsca/ec/mod/flint.py b/pyecsca/ec/mod/flint.py
new file mode 100644
index 0000000..d5fd432
--- /dev/null
+++ b/pyecsca/ec/mod/flint.py
@@ -0,0 +1,231 @@
+import warnings
+from functools import lru_cache, wraps
+from typing import Union
+
+from public import public
+
+from pyecsca.ec.error import (
+    raise_non_invertible,
+    raise_non_residue,
+    NonResidueError,
+    NonResidueWarning,
+)
+from pyecsca.ec.mod.base import Mod
+
+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
+
+
+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 flint."""
+
+        x: flint.fmpz_mod
+        _ctx: flint.fmpz_mod_ctx
+        __slots__ = ("x", "_ctx")
+
+        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, _=None) -> "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"])
+
+    from pyecsca.ec.mod.base import _mod_classes
+
+    _mod_classes["flint"] = FlintMod
diff --git a/pyecsca/ec/mod/gmp.py b/pyecsca/ec/mod/gmp.py
new file mode 100644
index 0000000..87301c0
--- /dev/null
+++ b/pyecsca/ec/mod/gmp.py
@@ -0,0 +1,173 @@
+from functools import lru_cache, wraps
+from typing import Union
+
+from public import public
+
+from pyecsca.ec.mod.base import Mod
+from pyecsca.ec.error import (
+    raise_non_invertible,
+    raise_non_residue,
+)
+
+has_gmp = False
+try:
+    import gmpy2
+
+    has_gmp = True
+except ImportError:
+    gmpy2 = None
+
+
+def _check(func):
+    @wraps(func)
+    def method(self, other):
+        if self.__class__ is not type(other):
+            other = self.__class__(other, self.n)
+        elif self.n != other.n:
+            raise ValueError
+        return func(self, other)
+
+    return method
+
+
+if has_gmp:
+
+    @lru_cache
+    def _gmpy_is_prime(x) -> bool:
+        return gmpy2.is_prime(x)
+
+    @public
+    class GMPMod(Mod):
+        """An element x of ℤₙ. Implemented by GMP."""
+
+        x: gmpy2.mpz
+        n: gmpy2.mpz
+        __slots__ = ("x", "n")
+
+        def __init__(
+            self,
+            x: Union[int, gmpy2.mpz],
+            n: Union[int, gmpy2.mpz],
+            ensure: bool = True,
+        ):
+            if ensure:
+                self.n = gmpy2.mpz(n)
+                self.x = gmpy2.mpz(x % self.n)
+            else:
+                self.n = n
+                self.x = x
+
+        def bit_length(self):
+            return self.x.bit_length()
+
+        def inverse(self) -> "GMPMod":
+            if self.x == 0:
+                raise_non_invertible()
+            if self.x == 1:
+                return GMPMod(gmpy2.mpz(1), self.n, ensure=False)
+            try:
+                res = gmpy2.invert(self.x, self.n)
+            except ZeroDivisionError:
+                raise_non_invertible()
+                res = gmpy2.mpz(0)
+            return GMPMod(res, self.n, ensure=False)
+
+        def is_residue(self) -> bool:
+            if not _gmpy_is_prime(self.n):
+                raise NotImplementedError
+            if self.x == 0:
+                return True
+            if self.n == 2:
+                return self.x in (0, 1)
+            return gmpy2.legendre(self.x, self.n) == 1
+
+        def sqrt(self) -> "GMPMod":
+            if not _gmpy_is_prime(self.n):
+                raise NotImplementedError
+            if self.x == 0:
+                return GMPMod(gmpy2.mpz(0), self.n, ensure=False)
+            if not self.is_residue():
+                raise_non_residue()
+            if self.n % 4 == 3:
+                return self ** int((self.n + 1) // 4)
+            q = self.n - 1
+            s = 0
+            while q % 2 == 0:
+                q //= 2
+                s += 1
+
+            z = gmpy2.mpz(2)
+            while GMPMod(z, self.n, ensure=False).is_residue():
+                z += 1
+
+            m = s
+            c = GMPMod(z, self.n, 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(GMPMod(gmpy2.mpz(2), self.n, ensure=False) ** two_exp)
+                m = int(GMPMod(gmpy2.mpz(i), self.n, ensure=False))
+                c = b**2
+                t *= c
+                r *= b
+            return r
+
+        @_check
+        def __add__(self, other) -> "GMPMod":
+            return GMPMod((self.x + other.x) % self.n, self.n, ensure=False)
+
+        @_check
+        def __sub__(self, other) -> "GMPMod":
+            return GMPMod((self.x - other.x) % self.n, self.n, ensure=False)
+
+        def __neg__(self) -> "GMPMod":
+            return GMPMod(self.n - self.x, self.n, ensure=False)
+
+        @_check
+        def __mul__(self, other) -> "GMPMod":
+            return GMPMod((self.x * other.x) % self.n, self.n, ensure=False)
+
+        def __bytes__(self):
+            return int(self.x).to_bytes((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 == (gmpy2.mpz(other) % self.n)
+            if type(other) is not GMPMod:
+                return False
+            return self.x == other.x and self.n == other.n
+
+        def __ne__(self, other):
+            return not self == other
+
+        def __repr__(self):
+            return str(int(self.x))
+
+        def __hash__(self):
+            return hash(("GMPMod", self.x, self.n))
+
+        def __pow__(self, n, _=None) -> "GMPMod":
+            if type(n) not in (int, gmpy2.mpz):
+                raise TypeError
+            if n == 0:
+                return GMPMod(gmpy2.mpz(1), self.n, ensure=False)
+            if n < 0:
+                return self.inverse() ** (-n)
+            if n == 1:
+                return GMPMod(self.x, self.n, ensure=False)
+            return GMPMod(
+                gmpy2.powmod(self.x, gmpy2.mpz(n), self.n), self.n, ensure=False
+            )
+
+    from pyecsca.ec.mod.base import _mod_classes
+
+    _mod_classes["gmp"] = GMPMod
diff --git a/pyecsca/ec/mod/raw.py b/pyecsca/ec/mod/raw.py
new file mode 100644
index 0000000..90e0f9d
--- /dev/null
+++ b/pyecsca/ec/mod/raw.py
@@ -0,0 +1,117 @@
+from public import public
+from pyecsca.ec.error import (
+    raise_non_invertible,
+    raise_non_residue,
+)
+
+from pyecsca.ec.mod.base import Mod, extgcd, miller_rabin, jacobi
+
+
+@public
+class RawMod(Mod):
+    """An element x of ℤₙ (implemented using Python integers)."""
+
+    x: int
+    n: int
+    __slots__ = ("x", "n")
+
+    def __init__(self, x: int, n: int):
+        self.x = x % n
+        self.n = n
+
+    def bit_length(self):
+        return self.x.bit_length()
+
+    def inverse(self) -> "RawMod":
+        if self.x == 0:
+            raise_non_invertible()
+        x, _, d = extgcd(self.x, self.n)
+        if d != 1:
+            raise_non_invertible()
+        return RawMod(x, self.n)
+
+    def is_residue(self):
+        if not miller_rabin(self.n):
+            raise NotImplementedError
+        if self.x == 0:
+            return True
+        if self.n == 2:
+            return self.x in (0, 1)
+        legendre_symbol = jacobi(self.x, self.n)
+        return legendre_symbol == 1
+
+    def sqrt(self) -> "RawMod":
+        if not miller_rabin(self.n):
+            raise NotImplementedError
+        if self.x == 0:
+            return RawMod(0, self.n)
+        if not self.is_residue():
+            raise_non_residue()
+        if self.n % 4 == 3:
+            return self ** int((self.n + 1) // 4)
+        q = self.n - 1
+        s = 0
+        while q % 2 == 0:
+            q //= 2
+            s += 1
+
+        z = 2
+        while RawMod(z, self.n).is_residue():
+            z += 1
+
+        m = s
+        c = RawMod(z, self.n) ** q
+        t = self**q
+        r_exp = (q + 1) // 2
+        r = self**r_exp
+
+        while t != 1:
+            i = 1
+            while not (t ** (2**i)) == 1:
+                i += 1
+            two_exp = m - (i + 1)
+            b = c ** int(RawMod(2, self.n) ** two_exp)
+            m = int(RawMod(i, self.n))
+            c = b**2
+            t *= c
+            r *= b
+        return r
+
+    def __bytes__(self):
+        return self.x.to_bytes((self.n.bit_length() + 7) // 8, byteorder="big")
+
+    def __int__(self):
+        return self.x
+
+    def __eq__(self, other):
+        if type(other) is int:
+            return self.x == (other % self.n)
+        if type(other) is not RawMod:
+            return False
+        return self.x == other.x and self.n == other.n
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __repr__(self):
+        return str(self.x)
+
+    def __hash__(self):
+        return hash(("RawMod", self.x, self.n))
+
+    def __pow__(self, n, _=None) -> "RawMod":
+        if type(n) is not int:
+            raise TypeError
+        if n == 0:
+            return RawMod(1, self.n)
+        if n < 0:
+            return self.inverse() ** (-n)
+        if n == 1:
+            return RawMod(self.x, self.n)
+
+        return RawMod(pow(self.x, n, self.n), self.n)
+
+
+from pyecsca.ec.mod.base import _mod_classes  # noqa
+
+_mod_classes["python"] = RawMod
diff --git a/pyecsca/ec/mod/symbolic.py b/pyecsca/ec/mod/symbolic.py
new file mode 100644
index 0000000..57e7233
--- /dev/null
+++ b/pyecsca/ec/mod/symbolic.py
@@ -0,0 +1,119 @@
+from functools import wraps
+
+from public import public
+from sympy import Expr
+
+from pyecsca.ec.mod.base import Mod
+
+
+def _check(func):
+    @wraps(func)
+    def method(self, other):
+        if self.__class__ is not type(other):
+            other = self.__class__(other, self.n)
+        elif self.n != other.n:
+            raise ValueError
+        return func(self, other)
+
+    return method
+
+
+@public
+class SymbolicMod(Mod):
+    """A symbolic element x of ℤₙ (implemented using sympy)."""
+
+    x: Expr
+    n: int
+    __slots__ = ("x", "n")
+
+    def __init__(self, x: Expr, n: int):
+        self.x = x
+        self.n = n
+
+    @_check
+    def __add__(self, other) -> "SymbolicMod":
+        return self.__class__((self.x + other.x), self.n)
+
+    @_check
+    def __radd__(self, other) -> "SymbolicMod":
+        return self + other
+
+    @_check
+    def __sub__(self, other) -> "SymbolicMod":
+        return self.__class__((self.x - other.x), self.n)
+
+    @_check
+    def __rsub__(self, other) -> "SymbolicMod":
+        return -self + other
+
+    def __neg__(self) -> "SymbolicMod":
+        return self.__class__(-self.x, self.n)
+
+    def bit_length(self):
+        raise NotImplementedError
+
+    def inverse(self) -> "SymbolicMod":
+        return self.__class__(self.x ** (-1), self.n)
+
+    def sqrt(self) -> "SymbolicMod":
+        raise NotImplementedError
+
+    def is_residue(self):
+        raise NotImplementedError
+
+    def __invert__(self) -> "SymbolicMod":
+        return self.inverse()
+
+    @_check
+    def __mul__(self, other) -> "SymbolicMod":
+        return self.__class__(self.x * other.x, self.n)
+
+    @_check
+    def __rmul__(self, other) -> "SymbolicMod":
+        return self * other
+
+    @_check
+    def __truediv__(self, other) -> "SymbolicMod":
+        return self * ~other
+
+    @_check
+    def __rtruediv__(self, other) -> "SymbolicMod":
+        return ~self * other
+
+    @_check
+    def __floordiv__(self, other) -> "SymbolicMod":
+        return self * ~other
+
+    @_check
+    def __rfloordiv__(self, other) -> "SymbolicMod":
+        return ~self * other
+
+    def __bytes__(self):
+        return int(self.x).to_bytes((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 % self.n
+        if type(other) is not SymbolicMod:
+            return False
+        return self.x == other.x and self.n == other.n
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __repr__(self):
+        return str(self.x)
+
+    def __hash__(self):
+        return hash(("SymbolicMod", self.x, self.n))
+
+    def __pow__(self, n, _=None) -> "SymbolicMod":
+        return self.__class__(pow(self.x, n), self.n)
+
+
+from pyecsca.ec.mod.base import _mod_classes  # noqa
+
+_mod_classes["symbolic"] = SymbolicMod
diff --git a/pyecsca/ec/params.py b/pyecsca/ec/params.py
index 6693581..2e437a2 100644
--- a/pyecsca/ec/params.py
+++ b/pyecsca/ec/params.py
@@ -18,7 +18,7 @@ from pyecsca.misc.cache import sympify
 from pyecsca.ec.coordinates import AffineCoordinateModel, CoordinateModel
 from pyecsca.ec.curve import EllipticCurve
 from pyecsca.ec.error import raise_unsatisified_assumption
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.model import (
     CurveModel,
     ShortWeierstrassModel,
@@ -189,7 +189,7 @@ def _create_params(curve, coords, infty):
     else:
         raise ValueError("Unknown curve model.")
     params = {
-        name: Mod(int(curve["params"][name]["raw"], 16), field) for name in param_names
+        name: mod(int(curve["params"][name]["raw"], 16), field) for name in param_names
     }
 
     # Check coordinate model name and assumptions
@@ -233,7 +233,7 @@ def _create_params(curve, coords, infty):
                 poly = Poly(numerator, symbols(param), domain=k)
                 roots = poly.ground_roots()
                 for root in roots:
-                    params[param] = Mod(int(k.from_sympy(root)), field)
+                    params[param] = mod(int(k.from_sympy(root)), field)
                     break
                 else:
                     raise_unsatisified_assumption(
@@ -258,16 +258,17 @@ def _create_params(curve, coords, infty):
                 )
             value = ilocals[coordinate]
             if isinstance(value, int):
-                value = Mod(value, field)
-            infinity_coords[coordinate] = value
+                infinity_coords[coordinate] = mod(value, field)
+            else:
+                infinity_coords[coordinate] = value
         infinity = Point(coord_model, **infinity_coords)
     elliptic_curve = EllipticCurve(model, coord_model, field, infinity, params)  # type: ignore[arg-type]
     if "generator" not in curve:
         raise ValueError("Cannot construct curve, missing generator.")
     affine = Point(
         AffineCoordinateModel(model),
-        x=Mod(int(curve["generator"]["x"]["raw"], 16), field),
-        y=Mod(int(curve["generator"]["y"]["raw"], 16), field),
+        x=mod(int(curve["generator"]["x"]["raw"], 16), field),
+        y=mod(int(curve["generator"]["y"]["raw"], 16), field),
     )
     if not isinstance(coord_model, AffineCoordinateModel):
         generator = affine.to_model(coord_model, elliptic_curve)
diff --git a/pyecsca/ec/point.py b/pyecsca/ec/point.py
index e8d5d5f..1b11e3f 100644
--- a/pyecsca/ec/point.py
+++ b/pyecsca/ec/point.py
@@ -6,7 +6,7 @@ from public import public
 
 from pyecsca.ec.context import ResultAction
 from pyecsca.ec.coordinates import AffineCoordinateModel, CoordinateModel
-from pyecsca.ec.mod import Mod, Undefined
+from pyecsca.ec.mod import Mod, Undefined, mod
 from pyecsca.ec.op import CodeOp
 
 
@@ -131,7 +131,7 @@ class Point:
             for var in coordinate_model.variables:
                 if var in locls:
                     result[var] = (
-                        Mod(locls[var], curve.prime)
+                        mod(locls[var], curve.prime)  # type: ignore
                         if not isinstance(locls[var], Mod)
                         else locls[var]
                     )
diff --git a/pyecsca/ec/signature.py b/pyecsca/ec/signature.py
index 331af6e..088c931 100644
--- a/pyecsca/ec/signature.py
+++ b/pyecsca/ec/signature.py
@@ -7,7 +7,7 @@ from public import public
 
 from pyecsca.ec.context import Action
 from pyecsca.ec.formula import AdditionFormula
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import ScalarMultiplier
 from pyecsca.ec.params import DomainParameters
 from pyecsca.ec.point import Point
@@ -165,7 +165,7 @@ class Signature:
         if nonce is None:
             return Mod.random(self.params.order)
         else:
-            return Mod(nonce, self.params.order)
+            return mod(nonce, self.params.order)
 
     def _do_sign(self, nonce: Mod, digest: bytes) -> SignatureResult:
         z = int.from_bytes(digest, byteorder="big")
@@ -174,8 +174,8 @@ class Signature:
         self.mult.init(self.params, self.params.generator)
         point = self.mult.multiply(int(nonce))
         affine_point = point.to_affine()
-        r = Mod(int(affine_point.x), self.params.order)
-        s = nonce.inverse() * (Mod(z, self.params.order) + r * self.privkey)
+        r = mod(int(affine_point.x), self.params.order)
+        s = nonce.inverse() * (mod(z, self.params.order) + r * self.privkey)
         return SignatureResult(int(r), int(s))
 
     def sign_hash(self, digest: bytes, nonce: Optional[int] = None) -> SignatureResult:
@@ -203,16 +203,16 @@ class Signature:
         z = int.from_bytes(digest, byteorder="big")
         if len(digest) * 8 > self.params.order.bit_length():
             z >>= len(digest) * 8 - self.params.order.bit_length()
-        c = Mod(signature.s, self.params.order).inverse()
-        u1 = Mod(z, self.params.order) * c
-        u2 = Mod(signature.r, self.params.order) * c
+        c = mod(signature.s, self.params.order).inverse()
+        u1 = mod(z, self.params.order) * c
+        u2 = mod(signature.r, self.params.order) * c
         self.mult.init(self.params, self.params.generator)
         p1 = self.mult.multiply(int(u1))
         self.mult.init(self.params, self.pubkey)
         p2 = self.mult.multiply(int(u2))
         p = self.add(self.params.curve.prime, p1, p2, **self.params.curve.parameters)[0]
         affine = p.to_affine()
-        v = Mod(int(affine.x), self.params.order)
+        v = mod(int(affine.x), self.params.order)
         return signature.r == int(v)
 
     def verify_hash(self, signature: SignatureResult, digest: bytes) -> bool:
diff --git a/pyecsca/ec/transformations.py b/pyecsca/ec/transformations.py
index 20ecde7..818a6dd 100644
--- a/pyecsca/ec/transformations.py
+++ b/pyecsca/ec/transformations.py
@@ -6,7 +6,7 @@ from sympy import FF, symbols, Poly
 
 from pyecsca.ec.coordinates import AffineCoordinateModel
 from pyecsca.ec.curve import EllipticCurve
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.model import (
     ShortWeierstrassModel,
     MontgomeryModel,
@@ -210,7 +210,7 @@ def __sw_ab(params: DomainParameters) -> Generator[Tuple[Mod, Mod], None, None]:
     if not roots:
         raise ValueError("Curve cannot be transformed (x^3 + ax + b has no root).")
     for root in roots:
-        alpha = Mod(int(root), params.curve.prime)
+        alpha = mod(int(root), params.curve.prime)
         beta = (3 * alpha**2 + params.curve.parameters["a"]).sqrt()
         yield alpha, beta
 
@@ -255,11 +255,11 @@ def SW2TE(params: DomainParameters) -> DomainParameters:
 
     def map_point(a, b, pt, aff):
         if params.curve.is_neutral(pt):
-            u = Mod(0, params.curve.prime)
-            v = Mod(1, params.curve.prime)
+            u = mod(0, params.curve.prime)
+            v = mod(1, params.curve.prime)
         elif pt.x == alpha and pt.y == 0:
-            u = Mod(0, params.curve.prime)
-            v = Mod(-1, params.curve.prime)
+            u = mod(0, params.curve.prime)
+            v = mod(-1, params.curve.prime)
         else:
             u = (pt.x - alpha) / pt.y
             v = (pt.x - alpha - beta) / (pt.x - alpha + beta)
diff --git a/pyecsca/sca/re/rpa.py b/pyecsca/sca/re/rpa.py
index c055274..0456eac 100644
--- a/pyecsca/sca/re/rpa.py
+++ b/pyecsca/sca/re/rpa.py
@@ -20,7 +20,7 @@ from pyecsca.ec.formula import (
     DifferentialAdditionFormula,
     LadderFormula,
 )
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import (
     ScalarMultiplicationAction,
     PrecomputationAction,
@@ -129,13 +129,13 @@ def rpa_point_0y(params: DomainParameters) -> Optional[Point]:
         y = params.curve.parameters["b"].sqrt()
         # TODO: We can take the negative as well.
         return Point(
-            AffineCoordinateModel(params.curve.model), x=Mod(0, params.curve.prime), y=y
+            AffineCoordinateModel(params.curve.model), x=mod(0, params.curve.prime), y=y
         )
     elif isinstance(params.curve.model, MontgomeryModel):
         return Point(
             AffineCoordinateModel(params.curve.model),
-            x=Mod(0, params.curve.prime),
-            y=Mod(0, params.curve.prime),
+            x=mod(0, params.curve.prime),
+            y=mod(0, params.curve.prime),
         )
     else:
         raise NotImplementedError
@@ -155,15 +155,15 @@ def rpa_point_x0(params: DomainParameters) -> Optional[Point]:
         roots = poly.ground_roots()
         if not roots:
             return None
-        x = Mod(int(next(iter(roots.keys()))), params.curve.prime)
+        x = mod(int(next(iter(roots.keys()))), params.curve.prime)
         return Point(
-            AffineCoordinateModel(params.curve.model), x=x, y=Mod(0, params.curve.prime)
+            AffineCoordinateModel(params.curve.model), x=x, y=mod(0, params.curve.prime)
         )
     elif isinstance(params.curve.model, MontgomeryModel):
         return Point(
             AffineCoordinateModel(params.curve.model),
-            x=Mod(0, params.curve.prime),
-            y=Mod(0, params.curve.prime),
+            x=mod(0, params.curve.prime),
+            y=mod(0, params.curve.prime),
         )
     else:
         raise NotImplementedError
@@ -269,7 +269,7 @@ class RPA(RE):
                 # Go over the parents and map them to multiples of the base (plus-minus sign)
                 init_multiples = set(
                     map(
-                        lambda v: Mod(v, params.order),
+                        lambda v: mod(v, params.order),
                         (init_context.points[parent] for parent in init_parents),
                     )
                 )
@@ -283,7 +283,7 @@ class RPA(RE):
                 )
                 multiply_multiples = set(
                     map(
-                        lambda v: Mod(v, params.order),
+                        lambda v: mod(v, params.order),
                         (ctx.points[parent] for parent in multiply_parents),
                     )
                 )
diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index 8005382..4954baa 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -23,7 +23,7 @@ from pyecsca.ec.formula import (
 )
 from pyecsca.ec.formula.fake import FakePoint, FakeFormula
 from pyecsca.ec.formula.unroll import unroll_formula
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import ScalarMultiplier
 from pyecsca.ec.params import DomainParameters
 from pyecsca.ec.point import Point
@@ -444,7 +444,7 @@ def zvp_points(poly: Poly, curve: EllipticCurve, k: int, n: int) -> Set[Point]:
                 points.add(point)
     elif only_2:
         # if only_2, dlog sub is not necessary, then multiply with k_inverse to obtain target point
-        k_inv = Mod(k, n).inverse()
+        k_inv = mod(k, n).inverse()
         for point in solve_easy_dcp(eliminated, curve):
             inputs = {"x2": point.x, "y2": point.y, **curve.parameters}
             res = poly.eval([inputs[str(gen)] for gen in poly.gens])  # type: ignore[attr-defined]
@@ -471,7 +471,7 @@ def zvp_points(poly: Poly, curve: EllipticCurve, k: int, n: int) -> Set[Point]:
 
 def _deterministic_point_x(curve: EllipticCurve) -> int:
     """Obtain a "random" coordinate `x` on  given curve."""
-    x = Mod(1, curve.prime)
+    x = mod(1, curve.prime)
     while True:
         points = curve.affine_lift_x(x)
         if points:
@@ -502,7 +502,7 @@ def solve_easy_dcp(xonly_polynomial: Poly, curve: EllipticCurve) -> Set[Point]:
         roots = final.ground_roots().keys()
 
     for root in roots:
-        points.update(curve.affine_lift_x(Mod(int(root), curve.prime)))
+        points.update(curve.affine_lift_x(mod(int(root), curve.prime)))
     return points
 
 
@@ -533,7 +533,7 @@ def solve_hard_dcp(xonly_polynomial: Poly, curve: EllipticCurve, k: int) -> Set[
 
     # Finally lift the roots to find the points (if any)
     for root in roots:
-        points.update(curve.affine_lift_x(Mod(int(root), curve.prime)))
+        points.update(curve.affine_lift_x(mod(int(root), curve.prime)))
     return points
 
 
diff --git a/test/ec/perf_formula.py b/test/ec/perf_formula.py
index 60b2b7d..1092ab5 100755
--- a/test/ec/perf_formula.py
+++ b/test/ec/perf_formula.py
@@ -1,19 +1,27 @@
 #!/usr/bin/env python
 import click
 
-from pyecsca.ec.mod import has_gmp, has_flint
+from pyecsca.ec.mod.gmp import has_gmp
+from pyecsca.ec.mod.flint import 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", "raw")), default="py")
+@click.option(
+    "-p",
+    "--profiler",
+    type=click.Choice(("py", "c", "raw")),
+    default="py",
+    envvar="PROF",
+)
 @click.option(
     "-m",
     "--mod",
     type=click.Choice(("python", "gmp", "flint")),
     default="flint" if has_flint else "gmp" if has_gmp else "python",
+    envvar="MOD",
 )
 @click.option("-o", "--operations", type=click.INT, default=5000)
 @click.option(
@@ -35,7 +43,10 @@ def main(profiler, mod, operations, directory):
         )
         one_point = p256.generator
         with Profiler(
-            profiler, directory, f"formula_dbl2016rcb_p256_{operations}_{mod}"
+            profiler,
+            directory,
+            f"formula_dbl2016rcb_p256_{operations}_{mod}",
+            operations,
         ):
             for _ in range(operations):
                 one_point = dbl(p256.curve.prime, one_point, **p256.curve.parameters)[0]
@@ -44,7 +55,10 @@ def main(profiler, mod, operations, directory):
         )
         other_point = p256.generator
         with Profiler(
-            profiler, directory, f"formula_add2016rcb_p256_{operations}_{mod}"
+            profiler,
+            directory,
+            f"formula_add2016rcb_p256_{operations}_{mod}",
+            operations,
         ):
             for _ in range(operations):
                 one_point = add(
@@ -58,7 +72,10 @@ def main(profiler, mod, operations, directory):
         )
         eone_point = ed25519.generator
         with Profiler(
-            profiler, directory, f"formula_mdbl2008hwcd_ed25519_{operations}_{mod}"
+            profiler,
+            directory,
+            f"formula_mdbl2008hwcd_ed25519_{operations}_{mod}",
+            operations,
         ):
             for _ in range(operations):
                 dblg(ed25519.curve.prime, eone_point, **ed25519.curve.parameters)
diff --git a/test/ec/perf_mod.py b/test/ec/perf_mod.py
index 925360e..7783f09 100755
--- a/test/ec/perf_mod.py
+++ b/test/ec/perf_mod.py
@@ -1,18 +1,27 @@
 #!/usr/bin/env python
 import click
 
-from pyecsca.ec.mod import Mod, has_gmp, has_flint
+from pyecsca.ec.mod.base import Mod, mod as make_mod
+from pyecsca.ec.mod.gmp import has_gmp
+from pyecsca.ec.mod.flint import has_flint
 from pyecsca.misc.cfg import TemporaryConfig
 from test.utils import Profiler
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
+@click.option(
+    "-p",
+    "--profiler",
+    type=click.Choice(("py", "c", "raw")),
+    default="py",
+    envvar="PROF",
+)
 @click.option(
     "-m",
     "--mod",
     type=click.Choice(("python", "gmp", "flint")),
     default="flint" if has_flint else "gmp" if has_gmp else "python",
+    envvar="MOD",
 )
 @click.option("-o", "--operations", type=click.INT, default=100000)
 @click.option(
@@ -26,21 +35,27 @@ def main(profiler, mod, operations, directory):
     with TemporaryConfig() as cfg:
         cfg.ec.mod_implementation = mod
         n = 0xFFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
-        a = Mod(0x11111111111111111111111111111111, n)
-        b = Mod(0xBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB, n)
+        a = make_mod(0x11111111111111111111111111111111, n)
+        b = make_mod(0xBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB, n)
         click.echo(f"Profiling {operations} {n.bit_length()}-bit modular inverse...")
-        with Profiler(profiler, directory, f"mod_256b_inverse_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_inverse_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 a.inverse()
         click.echo(
             f"Profiling {operations} {n.bit_length()}-bit modular square root..."
         )
-        with Profiler(profiler, directory, f"mod_256b_sqrt_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_sqrt_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 a.sqrt()
         click.echo(f"Profiling {operations} {n.bit_length()}-bit modular multiply...")
         c = a
-        with Profiler(profiler, directory, f"mod_256b_multiply_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_multiply_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 c = c * b
         click.echo(
@@ -48,33 +63,46 @@ def main(profiler, mod, operations, directory):
         )
         c = a
         with Profiler(
-            profiler, directory, f"mod_256b_constmultiply_{operations}_{mod}"
+            profiler,
+            directory,
+            f"mod_256b_constmultiply_{operations}_{mod}",
+            operations,
         ):
             for _ in range(operations):
                 c = c * 48006
         click.echo(f"Profiling {operations} {n.bit_length()}-bit modular square...")
         c = a
-        with Profiler(profiler, directory, f"mod_256b_square_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_square_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
-                c = c ** 2
+                c = c**2
         click.echo(f"Profiling {operations} {n.bit_length()}-bit modular add...")
         c = a
-        with Profiler(profiler, directory, f"mod_256b_add_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_add_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 c = c + b
         click.echo(f"Profiling {operations} {n.bit_length()}-bit modular subtract...")
         c = a
-        with Profiler(profiler, directory, f"mod_256b_subtract_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_subtract_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 c = c - b
         click.echo(
             f"Profiling {operations} {n.bit_length()}-bit modular quadratic residue checks..."
         )
-        with Profiler(profiler, directory, f"mod_256b_isresidue_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_isresidue_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 a.is_residue()
         click.echo(f"Profiling {operations} {n.bit_length()}-bit modular random...")
-        with Profiler(profiler, directory, f"mod_256b_random_{operations}_{mod}"):
+        with Profiler(
+            profiler, directory, f"mod_256b_random_{operations}_{mod}", operations
+        ):
             for _ in range(operations):
                 Mod.random(n)
 
diff --git a/test/ec/perf_mult.py b/test/ec/perf_mult.py
index 50a7191..402aa26 100755
--- a/test/ec/perf_mult.py
+++ b/test/ec/perf_mult.py
@@ -5,7 +5,8 @@ import click
 
 from pyecsca.ec.context import local, DefaultContext
 from pyecsca.ec.formula import AdditionFormula, DoublingFormula
-from pyecsca.ec.mod import has_gmp, has_flint
+from pyecsca.ec.mod.gmp import has_gmp
+from pyecsca.ec.mod.flint import has_flint
 from pyecsca.ec.mult import LTRMultiplier
 from pyecsca.ec.params import get_params
 from pyecsca.misc.cfg import TemporaryConfig
@@ -40,7 +41,7 @@ def main(profiler, mod, operations, directory):
             f"Profiling {operations} {p256.curve.prime.bit_length()}-bit scalar multiplication executions..."
         )
         one_point = p256.generator
-        with Profiler(profiler, directory, f"mult_ltr_rcb_p256_{operations}_{mod}"):
+        with Profiler(profiler, directory, f"mult_ltr_rcb_p256_{operations}_{mod}", operations):
             for _ in range(operations):
                 mult.init(p256, one_point)
                 one_point = mult.multiply(
@@ -51,7 +52,7 @@ def main(profiler, mod, operations, directory):
         )
         with local(DefaultContext()):
             one_point = p256.generator
-            with Profiler(profiler, directory, f"mult_ltr_rcb_p256_wtrace_{operations}_{mod}"):
+            with Profiler(profiler, directory, f"mult_ltr_rcb_p256_wtrace_{operations}_{mod}", operations):
                 for _ in range(operations):
                     mult.init(p256, one_point)
                     one_point = mult.multiply(
diff --git a/test/ec/test_curve.py b/test/ec/test_curve.py
index cd0639a..42bbcf6 100644
--- a/test/ec/test_curve.py
+++ b/test/ec/test_curve.py
@@ -4,7 +4,7 @@ import pytest
 from pyecsca.ec.coordinates import AffineCoordinateModel
 from pyecsca.ec.curve import EllipticCurve
 from pyecsca.ec.error import UnsatisfiedAssumptionError
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.model import MontgomeryModel
 from pyecsca.ec.point import Point, InfinityPoint
 
@@ -34,7 +34,7 @@ def test_init(secp128r1):
             secp128r1.curve.coordinate_model,
             15,
             InfinityPoint(secp128r1.curve.coordinate_model),
-            parameters={"a": Mod(1, 5), "b": Mod(2, 5)},
+            parameters={"a": mod(1, 5), "b": mod(2, 5)},
         )
 
 
@@ -64,17 +64,17 @@ def test_is_on_curve(secp128r1, curve25519):
     assert secp128r1.curve.is_on_curve(secp128r1.curve.neutral)
     pt = Point(
         secp128r1.curve.coordinate_model,
-        X=Mod(0x161FF7528B899B2D0C28607CA52C5B86, secp128r1.curve.prime),
-        Y=Mod(0xCF5AC8395BAFEB13C02DA292DDED7A83, secp128r1.curve.prime),
-        Z=Mod(1, secp128r1.curve.prime),
+        X=mod(0x161FF7528B899B2D0C28607CA52C5B86, secp128r1.curve.prime),
+        Y=mod(0xCF5AC8395BAFEB13C02DA292DDED7A83, secp128r1.curve.prime),
+        Z=mod(1, secp128r1.curve.prime),
     )
     assert secp128r1.curve.is_on_curve(pt)
     assert secp128r1.curve.is_on_curve(pt.to_affine())
     other = Point(
         secp128r1.curve.coordinate_model,
-        X=Mod(0x161FF7528B899B2D0C28607CA52C5B86, secp128r1.curve.prime),
-        Y=Mod(0xAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA, secp128r1.curve.prime),
-        Z=Mod(1, secp128r1.curve.prime),
+        X=mod(0x161FF7528B899B2D0C28607CA52C5B86, secp128r1.curve.prime),
+        Y=mod(0xAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA, secp128r1.curve.prime),
+        Z=mod(1, secp128r1.curve.prime),
     )
     assert not secp128r1.curve.is_on_curve(other)
     assert not secp128r1.curve.is_on_curve(curve25519.generator)
@@ -83,8 +83,8 @@ def test_is_on_curve(secp128r1, curve25519):
 def test_affine_add(secp128r1):
     pt = Point(
         AffineCoordinateModel(secp128r1.curve.model),
-        x=Mod(0xEB916224EDA4FB356421773573297C15, secp128r1.curve.prime),
-        y=Mod(0xBCDAF32A2C08FD4271228FEF35070848, secp128r1.curve.prime),
+        x=mod(0xEB916224EDA4FB356421773573297C15, secp128r1.curve.prime),
+        y=mod(0xBCDAF32A2C08FD4271228FEF35070848, secp128r1.curve.prime),
     )
     affine_base = secp128r1.generator.to_affine()
     assert secp128r1.curve.affine_add(affine_base, pt) is not None
diff --git a/test/ec/test_formula.py b/test/ec/test_formula.py
index 7b0bee1..4f64d36 100644
--- a/test/ec/test_formula.py
+++ b/test/ec/test_formula.py
@@ -21,7 +21,7 @@ from pyecsca.ec.formula.partitions import (
     generate_partitioned_formulas,
 )
 from pyecsca.ec.formula.switch_sign import generate_switched_formulas
-from pyecsca.ec.mod import SymbolicMod, Mod
+from pyecsca.ec.mod import SymbolicMod, Mod, mod
 from pyecsca.misc.cfg import TemporaryConfig
 from pyecsca.ec.error import UnsatisfiedAssumptionError
 from pyecsca.ec.params import get_params, DomainParameters
@@ -157,7 +157,7 @@ def test_symbolic(secp128r1, dbl):
             symbolic_val = symbolic_val.subs(
                 inner_var, int(getattr(secp128r1.generator, inner_var).x)
             )
-        assert Mod(int(symbolic_val), p) == Mod(generator_val, p)
+        assert mod(int(symbolic_val), p) == mod(generator_val, p)
 
 
 def test_pickle(add, dbl):
diff --git a/test/ec/test_key_agreement.py b/test/ec/test_key_agreement.py
index 248d172..60e48f8 100644
--- a/test/ec/test_key_agreement.py
+++ b/test/ec/test_key_agreement.py
@@ -13,7 +13,7 @@ from pyecsca.ec.key_agreement import (
     ECDH_SHA384,
     ECDH_SHA512,
 )
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import LTRMultiplier
 import test.data.ec
 from pyecsca.ec.params import get_params
@@ -29,7 +29,7 @@ def mult(secp128r1):
 
 @pytest.fixture()
 def keypair_a(secp128r1, mult):
-    priv_a = Mod(0xDEADBEEF, secp128r1.order)
+    priv_a = mod(0xDEADBEEF, secp128r1.order)
     mult.init(secp128r1, secp128r1.generator)
     pub_a = mult.multiply(int(priv_a))
     return priv_a, pub_a
@@ -37,7 +37,7 @@ def keypair_a(secp128r1, mult):
 
 @pytest.fixture()
 def keypair_b(secp128r1, mult):
-    priv_b = Mod(0xCAFEBABE, secp128r1.order)
+    priv_b = mod(0xCAFEBABE, secp128r1.order)
     mult.init(secp128r1, secp128r1.generator)
     pub_b = mult.multiply(int(priv_b))
     return priv_b, pub_b
@@ -58,15 +58,15 @@ def test_ka_secg():
     add = secp160r1.curve.coordinate_model.formulas["add-2015-rcb"]
     dbl = secp160r1.curve.coordinate_model.formulas["dbl-2015-rcb"]
     mult = LTRMultiplier(add, dbl)
-    privA = Mod(int(secg_data["keyA"]["priv"], 16), secp160r1.order)
+    privA = mod(int(secg_data["keyA"]["priv"], 16), secp160r1.order)
     pubA_affine = Point(affine_model,
-                        x=Mod(int(secg_data["keyA"]["pub"]["x"], 16), secp160r1.curve.prime),
-                        y=Mod(int(secg_data["keyA"]["pub"]["y"], 16), secp160r1.curve.prime))
+                        x=mod(int(secg_data["keyA"]["pub"]["x"], 16), secp160r1.curve.prime),
+                        y=mod(int(secg_data["keyA"]["pub"]["y"], 16), secp160r1.curve.prime))
     pubA = pubA_affine.to_model(secp160r1.curve.coordinate_model, secp160r1.curve)
-    privB = Mod(int(secg_data["keyB"]["priv"], 16), secp160r1.order)
+    privB = mod(int(secg_data["keyB"]["priv"], 16), secp160r1.order)
     pubB_affine = Point(affine_model,
-                        x=Mod(int(secg_data["keyB"]["pub"]["x"], 16), secp160r1.curve.prime),
-                        y=Mod(int(secg_data["keyB"]["pub"]["y"], 16), secp160r1.curve.prime))
+                        x=mod(int(secg_data["keyB"]["pub"]["x"], 16), secp160r1.curve.prime),
+                        y=mod(int(secg_data["keyB"]["pub"]["y"], 16), secp160r1.curve.prime))
     pubB = pubB_affine.to_model(secp160r1.curve.coordinate_model, secp160r1.curve)
 
     algoAB = ECDH_SHA1(copy(mult), secp160r1, pubA, privB)
diff --git a/test/ec/test_mod.py b/test/ec/test_mod.py
index 9a6cef0..13f8850 100644
--- a/test/ec/test_mod.py
+++ b/test/ec/test_mod.py
@@ -4,16 +4,17 @@ import pytest
 from sympy import FF, symbols
 
 from pyecsca.ec.mod import (
+    mod,
     Mod,
     gcd,
     extgcd,
     Undefined,
     miller_rabin,
-    has_gmp,
     RawMod,
     SymbolicMod,
     jacobi,
 )
+from pyecsca.ec.mod.gmp import has_gmp
 from pyecsca.ec.error import (
     NonInvertibleError,
     NonResidueError,
@@ -46,120 +47,120 @@ def test_miller_rabin():
 
 def test_inverse():
     p = 0xFFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
-    assert Mod(
+    assert mod(
         0x702BDAFD3C1C837B23A1CB196ED7F9FADB333C5CFE4A462BE32ADCD67BFB6AC1, p
-    ).inverse() == Mod(0x1CB2E5274BBA085C4CA88EEDE75AE77949E7A410C80368376E97AB22EB590F9D, p)
+    ).inverse() == mod(0x1CB2E5274BBA085C4CA88EEDE75AE77949E7A410C80368376E97AB22EB590F9D, p)
     with pytest.raises(NonInvertibleError):
-        Mod(0, p).inverse()
+        mod(0, p).inverse()
     with pytest.raises(NonInvertibleError):
-        Mod(5, 10).inverse()
+        mod(5, 10).inverse()
     getconfig().ec.no_inverse_action = "warning"
     with warnings.catch_warnings(record=True) as w:
-        Mod(0, p).inverse()
+        mod(0, p).inverse()
         assert issubclass(w[0].category, NonInvertibleWarning)
     with warnings.catch_warnings(record=True) as w:
-        Mod(5, 10).inverse()
+        mod(5, 10).inverse()
         assert issubclass(w[0].category, NonInvertibleWarning)
     getconfig().ec.no_inverse_action = "ignore"
-    Mod(0, p).inverse()
-    Mod(5, 10).inverse()
+    mod(0, p).inverse()
+    mod(5, 10).inverse()
     getconfig().ec.no_inverse_action = "error"
 
 
 def test_is_residue():
-    assert Mod(4, 11).is_residue()
-    assert not Mod(11, 31).is_residue()
-    assert Mod(0, 7).is_residue()
-    assert Mod(1, 2).is_residue()
+    assert mod(4, 11).is_residue()
+    assert not mod(11, 31).is_residue()
+    assert mod(0, 7).is_residue()
+    assert mod(1, 2).is_residue()
 
 
 def test_bit_length():
-    x = Mod(3, 5)
+    x = mod(3, 5)
     assert x.bit_length() == 2
 
 
 def test_sqrt():
     p = 0xFFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
-    assert Mod(
+    assert mod(
         0xFFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC, p
     ).sqrt() in (
                0x9ADD512515B70D9EC471151C1DEC46625CD18B37BDE7CA7FB2C8B31D7033599D,
                0x6522AED9EA48F2623B8EEAE3E213B99DA32E74C9421835804D374CE28FCCA662,
            )
     with pytest.raises(NonResidueError):
-        Mod(
+        mod(
             0x702BDAFD3C1C837B23A1CB196ED7F9FADB333C5CFE4A462BE32ADCD67BFB6AC1, p
         ).sqrt()
     getconfig().ec.non_residue_action = "warning"
     with warnings.catch_warnings(record=True) as w:
-        Mod(
+        mod(
             0x702BDAFD3C1C837B23A1CB196ED7F9FADB333C5CFE4A462BE32ADCD67BFB6AC1, p
         ).sqrt()
         assert issubclass(w[0].category, NonResidueWarning)
     getconfig().ec.non_residue_action = "ignore"
-    Mod(
+    mod(
         0x702BDAFD3C1C837B23A1CB196ED7F9FADB333C5CFE4A462BE32ADCD67BFB6AC1, p
     ).sqrt()
     with TemporaryConfig() as cfg:
         cfg.ec.non_residue_action = "warning"
         with warnings.catch_warnings(record=True) as w:
-            Mod(
+            mod(
                 0x702BDAFD3C1C837B23A1CB196ED7F9FADB333C5CFE4A462BE32ADCD67BFB6AC1,
                 p,
             ).sqrt()
             assert issubclass(w[0].category, NonResidueWarning)
-    assert Mod(0, p).sqrt() == Mod(0, p)
+    assert mod(0, p).sqrt() == mod(0, p)
     q = 0x75D44FEE9A71841AE8403C0C251FBAD
-    assert Mod(0x591E0DB18CF1BD81A11B2985A821EB3, q).sqrt() in \
+    assert mod(0x591E0DB18CF1BD81A11B2985A821EB3, q).sqrt() in \
            (0x113B41A1A2B73F636E73BE3F9A3716E, 0x64990E4CF7BA44B779CC7DCC8AE8A3F)
     getconfig().ec.non_residue_action = "error"
 
 
 def test_eq():
-    assert Mod(1, 7) == 1
-    assert Mod(1, 7) != "1"
-    assert Mod(1, 7) == Mod(1, 7)
-    assert Mod(1, 7) != Mod(5, 7)
-    assert Mod(1, 7) != Mod(1, 5)
+    assert mod(1, 7) == 1
+    assert mod(1, 7) != "1"
+    assert mod(1, 7) == mod(1, 7)
+    assert mod(1, 7) != mod(5, 7)
+    assert mod(1, 7) != mod(1, 5)
 
 
 def test_pow():
-    a = Mod(5, 7)
+    a = mod(5, 7)
     assert a ** (-1) == a.inverse()
-    assert a ** 0 == Mod(1, 7)
+    assert a ** 0 == mod(1, 7)
     assert a ** (-2) == a.inverse() ** 2
 
 
 def test_wrong_mod():
-    a = Mod(5, 7)
-    b = Mod(4, 11)
+    a = mod(5, 7)
+    b = mod(4, 11)
     with pytest.raises(ValueError):
         a + b
 
 
 def test_wrong_pow():
-    a = Mod(5, 7)
-    c = Mod(4, 11)
+    a = mod(5, 7)
+    c = mod(4, 11)
     with pytest.raises(TypeError):
         a ** c
 
 
 def test_other():
-    a = Mod(5, 7)
-    b = Mod(3, 7)
+    a = mod(5, 7)
+    b = mod(3, 7)
     assert int(-a) == 2
     assert str(a) == "5"
-    assert 6 - a == Mod(1, 7)
+    assert 6 - a == mod(1, 7)
     assert a != b
-    assert a / b == Mod(4, 7)
-    assert a // b == Mod(4, 7)
-    assert 5 / b == Mod(4, 7)
-    assert 5 // b == Mod(4, 7)
-    assert a / 3 == Mod(4, 7)
-    assert a // 3 == Mod(4, 7)
-    assert a + b == Mod(1, 7)
-    assert 5 + b == Mod(1, 7)
-    assert a + 3 == Mod(1, 7)
+    assert a / b == mod(4, 7)
+    assert a // b == mod(4, 7)
+    assert 5 / b == mod(4, 7)
+    assert 5 // b == mod(4, 7)
+    assert a / 3 == mod(4, 7)
+    assert a // 3 == mod(4, 7)
+    assert a + b == mod(1, 7)
+    assert 5 + b == mod(1, 7)
+    assert a + 3 == mod(1, 7)
     assert a != 6
     assert hash(a) is not None
 
@@ -198,7 +199,7 @@ def test_implementation():
         pytest.skip("Only makes sense if more Mod implementations are available.")
     with TemporaryConfig() as cfg:
         cfg.ec.mod_implementation = "python"
-        assert isinstance(Mod(5, 7), RawMod)
+        assert isinstance(mod(5, 7), RawMod)
 
 
 def test_symbolic():
diff --git a/test/ec/test_mult.py b/test/ec/test_mult.py
index fb10b13..e477b07 100644
--- a/test/ec/test_mult.py
+++ b/test/ec/test_mult.py
@@ -1,9 +1,9 @@
 from itertools import product
-from typing import Sequence
+from typing import Sequence, List
 
 import pytest
 
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import (
     DoubleAndAddMultiplier,
     LTRMultiplier,
@@ -191,21 +191,21 @@ def test_ladder(curve25519):
 )
 def test_ladder_full(curve25519, scalar, x, res):
     p = curve25519.curve.prime
-    point = Point(curve25519.curve.coordinate_model, X=Mod(x, p), Z=Mod(1, p))
-    result = Point(curve25519.curve.coordinate_model, X=Mod(res, p), Z=Mod(1, p))
+    point = Point(curve25519.curve.coordinate_model, X=mod(x, p), Z=mod(1, p))
+    result = Point(curve25519.curve.coordinate_model, X=mod(res, p), Z=mod(1, p))
 
     mult = LadderMultiplier(
         curve25519.curve.coordinate_model.formulas["ladd-1987-m"],
         curve25519.curve.coordinate_model.formulas["dbl-1987-m"],
         # complete=False
     )
-    fixed = int(Mod(scalar, curve25519.order))
+    fixed = int(mod(scalar, curve25519.order))
 
     mult.init(curve25519, point)
     computed = mult.multiply(fixed)
 
-    point_aff = list(curve25519.curve.affine_lift_x(Mod(x, p)))[0]
-    result_aff = list(curve25519.curve.affine_lift_x(Mod(res, p)))[0]
+    point_aff = list(curve25519.curve.affine_lift_x(mod(x, p)))[0]
+    result_aff = list(curve25519.curve.affine_lift_x(mod(res, p)))[0]
     computed_aff = curve25519.curve.affine_multiply(point_aff, scalar)
 
     scale = curve25519.curve.coordinate_model.formulas["scale"]
@@ -480,7 +480,7 @@ def test_basic_multipliers(secp128r1, num, add, dbl):
         + bgmws
         + combs
     )
-    results = []
+    results: List[Point] = []
     for mult in mults:
         mult.init(secp128r1, secp128r1.generator)
         res = mult.multiply(num)
diff --git a/test/ec/test_op.py b/test/ec/test_op.py
index ece8cd9..8d6996b 100644
--- a/test/ec/test_op.py
+++ b/test/ec/test_op.py
@@ -3,7 +3,7 @@ from ast import parse
 import pytest
 
 from pyecsca.ec.formula import OpResult
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.op import CodeOp, OpType
 
 
@@ -21,8 +21,8 @@ def test_str(name, module, result, op_type):
 
 
 @pytest.mark.parametrize("name,module,locals,result",
-                         [("add", "x = a+b", {"a": Mod(5, 21), "b": Mod(7, 21)}, Mod(12, 21)),
-                          ("sub", "x = a-b", {"a": Mod(7, 21), "b": Mod(5, 21)}, Mod(2, 21)), ])
+                         [("add", "x = a+b", {"a": mod(5, 21), "b": mod(7, 21)}, mod(12, 21)),
+                          ("sub", "x = a-b", {"a": mod(7, 21), "b": mod(5, 21)}, mod(2, 21)), ])
 def test_call(name, module, locals, result):
     code = parse(module, mode="exec")
     op = CodeOp(code)
@@ -31,9 +31,9 @@ def test_call(name, module, locals, result):
 
 
 def test_opresult_repr():
-    res = OpResult("a", Mod(7, 11), OpType.Neg, "b")
+    res = OpResult("a", mod(7, 11), OpType.Neg, "b")
     assert repr(res) == "a = -b"
-    res = OpResult("a", Mod(5, 7), OpType.Add, "c", 3)
+    res = OpResult("a", mod(5, 7), OpType.Add, "c", 3)
     assert repr(res) == "a = c+3"
-    res = OpResult("a", Mod(3, 11), OpType.Inv, "d")
+    res = OpResult("a", mod(3, 11), OpType.Inv, "d")
     assert repr(res) == "a = 1/d"
diff --git a/test/ec/test_params.py b/test/ec/test_params.py
index 8ca7252..891a7ca 100644
--- a/test/ec/test_params.py
+++ b/test/ec/test_params.py
@@ -5,7 +5,7 @@ from importlib_resources import files, as_file
 import pytest
 
 import test.data.ec
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.point import Point, InfinityPoint
 from pyecsca.misc.cfg import TemporaryConfig
 from pyecsca.ec.coordinates import AffineCoordinateModel
@@ -131,11 +131,11 @@ def test_custom_params():
     model = ShortWeierstrassModel()
     coords = model.coordinates["projective"]
     p = 0xd7d1247f
-    a = Mod(0xa4a44016, p)
-    b = Mod(0x73f76716, p)
+    a = mod(0xa4a44016, p)
+    b = mod(0x73f76716, p)
     n = 0xd7d2a475
     h = 1
-    gx, gy, gz = Mod(0x54eed6d7, p), Mod(0x6f1e55ac, p), Mod(1, p)
+    gx, gy, gz = mod(0x54eed6d7, p), mod(0x6f1e55ac, p), mod(1, p)
     generator = Point(coords, X=gx, Y=gy, Z=gz)
     neutral = InfinityPoint(coords)
 
diff --git a/test/ec/test_point.py b/test/ec/test_point.py
index 9513840..52c2a71 100644
--- a/test/ec/test_point.py
+++ b/test/ec/test_point.py
@@ -2,7 +2,7 @@ import pickle
 from contextlib import nullcontext as does_not_raise
 from pyecsca.ec.coordinates import AffineCoordinateModel
 from pyecsca.ec.params import get_params
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import mod
 from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel
 from pyecsca.ec.point import Point, InfinityPoint
 import pytest
@@ -22,15 +22,15 @@ def test_construction(secp128r1_coords):
     with pytest.raises(ValueError):
         Point(secp128r1_coords)
     with pytest.raises(ValueError):
-        Point(secp128r1_coords, X=Mod(1, 3), Y=Mod(2, 7), Z=Mod(1, 3))
+        Point(secp128r1_coords, X=mod(1, 3), Y=mod(2, 7), Z=mod(1, 3))
 
 
 def test_to_affine(secp128r1, secp128r1_coords, affine_model):
     pt = Point(
         secp128r1_coords,
-        X=Mod(0x161FF7528B899B2D0C28607CA52C5B86, secp128r1.curve.prime),
-        Y=Mod(0xCF5AC8395BAFEB13C02DA292DDED7A83, secp128r1.curve.prime),
-        Z=Mod(1, secp128r1.curve.prime),
+        X=mod(0x161FF7528B899B2D0C28607CA52C5B86, secp128r1.curve.prime),
+        Y=mod(0xCF5AC8395BAFEB13C02DA292DDED7A83, secp128r1.curve.prime),
+        Z=mod(1, secp128r1.curve.prime),
     )
     affine = pt.to_affine()
 
@@ -55,8 +55,8 @@ def test_to_affine(secp128r1, secp128r1_coords, affine_model):
 def test_to_model(secp128r1, secp128r1_coords, affine_model):
     affine = Point(
         affine_model,
-        x=Mod(0xABCD, secp128r1.curve.prime),
-        y=Mod(0xEF, secp128r1.curve.prime),
+        x=mod(0xABCD, secp128r1.curve.prime),
+        y=mod(0xEF, secp128r1.curve.prime),
     )
     other = affine.to_model(secp128r1_coords, secp128r1.curve)
 
@@ -64,7 +64,7 @@ def test_to_model(secp128r1, secp128r1_coords, affine_model):
     assert set(other.coords.keys()) == set(secp128r1_coords.variables)
     assert other.coords["X"] == affine.coords["x"]
     assert other.coords["Y"] == affine.coords["y"]
-    assert other.coords["Z"] == Mod(1, secp128r1.curve.prime)
+    assert other.coords["Z"] == mod(1, secp128r1.curve.prime)
 
     infty = InfinityPoint(AffineCoordinateModel(secp128r1.curve.model))
     other_infty = infty.to_model(secp128r1_coords, secp128r1.curve)
@@ -101,21 +101,21 @@ def test_to_from_affine(category, curve, coords, raises):
 def test_equals(secp128r1, secp128r1_coords):
     pt = Point(
         secp128r1_coords,
-        X=Mod(0x4, secp128r1.curve.prime),
-        Y=Mod(0x6, secp128r1.curve.prime),
-        Z=Mod(2, secp128r1.curve.prime),
+        X=mod(0x4, secp128r1.curve.prime),
+        Y=mod(0x6, secp128r1.curve.prime),
+        Z=mod(2, secp128r1.curve.prime),
     )
     other = Point(
         secp128r1_coords,
-        X=Mod(0x2, secp128r1.curve.prime),
-        Y=Mod(0x3, secp128r1.curve.prime),
-        Z=Mod(1, secp128r1.curve.prime),
+        X=mod(0x2, secp128r1.curve.prime),
+        Y=mod(0x3, secp128r1.curve.prime),
+        Z=mod(1, secp128r1.curve.prime),
     )
     third = Point(
         secp128r1_coords,
-        X=Mod(0x5, secp128r1.curve.prime),
-        Y=Mod(0x3, secp128r1.curve.prime),
-        Z=Mod(1, secp128r1.curve.prime),
+        X=mod(0x5, secp128r1.curve.prime),
+        Y=mod(0x3, secp128r1.curve.prime),
+        Z=mod(1, secp128r1.curve.prime),
     )
     assert pt.equals(other)
     assert pt != other
@@ -140,11 +140,11 @@ def test_equals(secp128r1, secp128r1_coords):
     mont = MontgomeryModel()
     different = Point(
         mont.coordinates["xz"],
-        X=Mod(
+        X=mod(
             0x64DACCD2656420216545E5F65221EB,
             0xAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA,
         ),
-        Z=Mod(1, 0xAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA),
+        Z=mod(1, 0xAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA),
     )
     assert not pt.equals(different)
     assert pt != different
@@ -153,9 +153,9 @@ def test_equals(secp128r1, secp128r1_coords):
 def test_bytes(secp128r1, secp128r1_coords):
     pt = Point(
         secp128r1_coords,
-        X=Mod(0x4, secp128r1.curve.prime),
-        Y=Mod(0x6, secp128r1.curve.prime),
-        Z=Mod(2, secp128r1.curve.prime),
+        X=mod(0x4, secp128r1.curve.prime),
+        Y=mod(0x6, secp128r1.curve.prime),
+        Z=mod(2, secp128r1.curve.prime),
     )
     assert bytes(pt) == \
            b"\x04\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x04\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x06\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x02"
@@ -165,9 +165,9 @@ def test_bytes(secp128r1, secp128r1_coords):
 def test_iter(secp128r1, secp128r1_coords):
     pt = Point(
         secp128r1_coords,
-        X=Mod(0x4, secp128r1.curve.prime),
-        Y=Mod(0x6, secp128r1.curve.prime),
-        Z=Mod(2, secp128r1.curve.prime),
+        X=mod(0x4, secp128r1.curve.prime),
+        Y=mod(0x6, secp128r1.curve.prime),
+        Z=mod(2, secp128r1.curve.prime),
     )
     t = tuple(pt)
     assert len(t) == 3
@@ -180,9 +180,9 @@ def test_iter(secp128r1, secp128r1_coords):
 def test_pickle(secp128r1, secp128r1_coords):
     pt = Point(
         secp128r1_coords,
-        X=Mod(0x4, secp128r1.curve.prime),
-        Y=Mod(0x6, secp128r1.curve.prime),
-        Z=Mod(2, secp128r1.curve.prime),
+        X=mod(0x4, secp128r1.curve.prime),
+        Y=mod(0x6, secp128r1.curve.prime),
+        Z=mod(2, secp128r1.curve.prime),
     )
     pickle.dumps(secp128r1_coords)
     assert pt == pickle.loads(pickle.dumps(pt))
diff --git a/test/ec/test_regress.py b/test/ec/test_regress.py
index 8d54e98..cbfb08b 100644
--- a/test/ec/test_regress.py
+++ b/test/ec/test_regress.py
@@ -7,7 +7,7 @@ from pyecsca.ec.coordinates import AffineCoordinateModel
 from pyecsca.ec.curve import EllipticCurve
 from pyecsca.ec.error import UnsatisfiedAssumptionError
 from pyecsca.ec.formula import AdditionFormula, DoublingFormula, ScalingFormula
-from pyecsca.ec.mod import Mod, SymbolicMod
+from pyecsca.ec.mod import Mod, SymbolicMod, mod
 from pyecsca.ec.model import MontgomeryModel, EdwardsModel
 from pyecsca.ec.params import get_params
 from pyecsca.ec.mult import LTRMultiplier
@@ -48,13 +48,13 @@ def test_issue_9():
     model = MontgomeryModel()
     coords = model.coordinates["xz"]
     p = 19
-    neutral = Point(coords, X=Mod(1, p), Z=Mod(0, p))
-    curve = EllipticCurve(model, coords, p, neutral, {"a": Mod(8, p), "b": Mod(1, p)})
-    base = Point(coords, X=Mod(12, p), Z=Mod(1, p))
+    neutral = Point(coords, X=mod(1, p), Z=mod(0, p))
+    curve = EllipticCurve(model, coords, p, neutral, {"a": mod(8, p), "b": mod(1, p)})
+    base = Point(coords, X=mod(12, p), Z=mod(1, p))
     formula = coords.formulas["dbl-1987-m-2"]
     res = formula(p, base, **curve.parameters)[0]
     assert res is not None
-    affine_base = Point(AffineCoordinateModel(model), x=Mod(12, p), y=Mod(2, p))
+    affine_base = Point(AffineCoordinateModel(model), x=mod(12, p), y=mod(2, p))
     dbase = curve.affine_double(affine_base).to_model(coords, curve)
     ladder = coords.formulas["ladd-1987-m-3"]
     one, other = ladder(p, base, dbase, base, **curve.parameters)
@@ -67,14 +67,14 @@ def test_issue_10():
     coords = model.coordinates["yz"]
     coords_sqr = model.coordinates["yzsquared"]
     p = 0x1D
-    c = Mod(1, p)
-    d = Mod(0x1C, p)
+    c = mod(1, p)
+    d = mod(0x1C, p)
     r = d.sqrt()
-    neutral = Point(coords, Y=c * r, Z=Mod(1, p))
+    neutral = Point(coords, Y=c * r, Z=mod(1, p))
     curve = EllipticCurve(model, coords, p, neutral, {"c": c, "d": d, "r": r})
-    neutral_affine = Point(AffineCoordinateModel(model), x=Mod(0, p), y=c)
+    neutral_affine = Point(AffineCoordinateModel(model), x=mod(0, p), y=c)
     assert neutral == neutral_affine.to_model(coords, curve)
-    neutral_sqr = Point(coords_sqr, Y=c ** 2 * r, Z=Mod(1, p))
+    neutral_sqr = Point(coords_sqr, Y=c ** 2 * r, Z=mod(1, p))
     assert neutral_sqr == neutral_affine.to_model(coords_sqr, curve)
 
 
@@ -103,23 +103,23 @@ def test_issue_14():
     with pytest.raises(UnsatisfiedAssumptionError):
         # p is 3 mod 4, so there is no square root of -1
         p = 19
-        c = Mod(2, p)
-        d = Mod(10, p)
+        c = mod(2, p)
+        d = mod(10, p)
         curve = EllipticCurve(model, coords, p, InfinityPoint(coords), {"c": c, "d": d})
-        Paff = Point(affine, x=Mod(0xD, p), y=Mod(0x9, p))
+        Paff = Point(affine, x=mod(0xD, p), y=mod(0x9, p))
         P = Paff.to_model(coords, curve)
-        Qaff = Point(affine, x=Mod(0x4, p), y=Mod(0x12, p))
+        Qaff = Point(affine, x=mod(0x4, p), y=mod(0x12, p))
         Q = Qaff.to_model(coords, curve)
         formula(p, P, Q, **curve.parameters)[0]
 
     # p is 1 mod 4, so there is a square root of -1
     p = 29
-    c = Mod(2, p)
-    d = Mod(10, p)
+    c = mod(2, p)
+    d = mod(10, p)
     curve = EllipticCurve(model, coords, p, InfinityPoint(coords), {"c": c, "d": d})
-    Paff = Point(affine, x=Mod(0xD, p), y=Mod(0x9, p))
+    Paff = Point(affine, x=mod(0xD, p), y=mod(0x9, p))
     P = Paff.to_model(coords, curve)
-    Qaff = Point(affine, x=Mod(0x4, p), y=Mod(0x12, p))
+    Qaff = Point(affine, x=mod(0x4, p), y=mod(0x12, p))
     Q = Qaff.to_model(coords, curve)
     PQaff = curve.affine_add(Paff, Qaff)
     R = formula(p, P, Q, **curve.parameters)[0]
diff --git a/test/ec/test_signature.py b/test/ec/test_signature.py
index 1244d6e..df00aa6 100644
--- a/test/ec/test_signature.py
+++ b/test/ec/test_signature.py
@@ -4,7 +4,7 @@ import pytest
 from importlib_resources import files
 import test.data.ec
 from pyecsca.ec.coordinates import AffineCoordinateModel
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import LTRMultiplier
 from pyecsca.ec.params import get_params
 from pyecsca.ec.point import Point
@@ -33,7 +33,7 @@ def mult(secp128r1, add):
 
 @pytest.fixture()
 def keypair(secp128r1, mult):
-    priv = Mod(0xDEADBEEF, secp128r1.order)
+    priv = mod(0xDEADBEEF, secp128r1.order)
     mult.init(secp128r1, secp128r1.generator)
     pub = mult.multiply(int(priv))
     return priv, pub
@@ -126,11 +126,11 @@ def test_ecdsa_nist():
     add = P192.curve.coordinate_model.formulas["add-2015-rcb"]
     dbl = P192.curve.coordinate_model.formulas["dbl-2015-rcb"]
     mult = LTRMultiplier(add, dbl)
-    priv = Mod(int(nist_data["priv"], 16), P192.order)
+    priv = mod(int(nist_data["priv"], 16), P192.order)
 
     pub_affine = Point(affine_model,
-                       x=Mod(int(nist_data["pub"]["x"], 16), P192.curve.prime),
-                       y=Mod(int(nist_data["pub"]["y"], 16), P192.curve.prime))
+                       x=mod(int(nist_data["pub"]["x"], 16), P192.curve.prime),
+                       y=mod(int(nist_data["pub"]["y"], 16), P192.curve.prime))
     pub = pub_affine.to_model(P192.curve.coordinate_model, P192.curve)
 
     signer = ECDSA_SHA1(mult, P192, add, pub, priv)
diff --git a/test/sca/perf_combine.py b/test/sca/perf_combine.py
index bd5db59..8ac274c 100644
--- a/test/sca/perf_combine.py
+++ b/test/sca/perf_combine.py
@@ -16,7 +16,13 @@ from pyecsca.sca import (
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
+@click.option(
+    "-p",
+    "--profiler",
+    type=click.Choice(("py", "c", "raw")),
+    default="py",
+    envvar="PROF",
+)
 @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 2175d49..34090df 100644
--- a/test/sca/perf_zvp.py
+++ b/test/sca/perf_zvp.py
@@ -1,7 +1,8 @@
 #!/usr/bin/env python
 import click
 
-from pyecsca.ec.mod import has_gmp, has_flint
+from pyecsca.ec.mod.gmp import has_gmp
+from pyecsca.ec.mod.flint import 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 +11,19 @@ from test.utils import Profiler
 
 
 @click.command()
-@click.option("-p", "--profiler", type=click.Choice(("py", "c", "raw")), default="py")
+@click.option(
+    "-p",
+    "--profiler",
+    type=click.Choice(("py", "c", "raw")),
+    default="py",
+    envvar="PROF",
+)
 @click.option(
     "-m",
     "--mod",
     type=click.Choice(("python", "gmp", "flint")),
     default="flint" if has_flint else "gmp" if has_gmp else "python",
+    envvar="MOD",
 )
 @click.option("-o", "--operations", type=click.INT, default=1)
 @click.option(
diff --git a/test/sca/test_leakage_models.py b/test/sca/test_leakage_models.py
index 592a100..08a3b07 100644
--- a/test/sca/test_leakage_models.py
+++ b/test/sca/test_leakage_models.py
@@ -1,6 +1,6 @@
 from pyecsca.ec.context import local, DefaultContext
 from pyecsca.ec.formula import FormulaAction, OpResult
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import LTRMultiplier
 from pyecsca.ec.op import OpType
 from pyecsca.sca.attack.leakage_model import Identity, Bit, Slice, HammingWeight, HammingDistance, BitLength
@@ -8,13 +8,13 @@ import pytest
 
 
 def test_identity():
-    val = Mod(3, 7)
+    val = mod(3, 7)
     lm = Identity()
     assert lm(val) == 3
 
 
 def test_bit():
-    val = Mod(3, 7)
+    val = mod(3, 7)
     lm = Bit(0)
     assert lm(val) == 1
     lm = Bit(4)
@@ -24,7 +24,7 @@ def test_bit():
 
 
 def test_slice():
-    val = Mod(0b11110000, 0xf00)
+    val = mod(0b11110000, 0xf00)
     lm = Slice(0, 4)
     assert lm(val) == 0
     lm = Slice(1, 5)
@@ -36,20 +36,20 @@ def test_slice():
 
 
 def test_hamming_weight():
-    val = Mod(0b11110000, 0xf00)
+    val = mod(0b11110000, 0xf00)
     lm = HammingWeight()
     assert lm(val) == 4
 
 
 def test_hamming_distance():
-    a = Mod(0b11110000, 0xf00)
-    b = Mod(0b00010000, 0xf00)
+    a = mod(0b11110000, 0xf00)
+    b = mod(0b00010000, 0xf00)
     lm = HammingDistance()
     assert lm(a, b) == 3
 
 
 def test_bit_length():
-    a = Mod(0b11110000, 0xf00)
+    a = mod(0b11110000, 0xf00)
     lm = BitLength()
     assert lm(a) == 8
 
diff --git a/test/sca/test_rpa.py b/test/sca/test_rpa.py
index fe38c41..2fa533d 100644
--- a/test/sca/test_rpa.py
+++ b/test/sca/test_rpa.py
@@ -3,7 +3,7 @@ import pytest
 from pyecsca.ec.context import local
 from pyecsca.ec.model import ShortWeierstrassModel
 from pyecsca.ec.curve import EllipticCurve
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import (
     LTRMultiplier,
     RTLMultiplier,
@@ -57,15 +57,15 @@ def neg(coords):
 @pytest.fixture()
 def rpa_params(model, coords):
     p = 0x85D265945A4F5681
-    a = Mod(0x7FC57B4110698BC0, p)
-    b = Mod(0x37113EA591B04527, p)
-    gx = Mod(0x80D2D78FDDB97597, p)
-    gy = Mod(0x5586D818B7910930, p)
+    a = mod(0x7FC57B4110698BC0, p)
+    b = mod(0x37113EA591B04527, p)
+    gx = mod(0x80D2D78FDDB97597, p)
+    gy = mod(0x5586D818B7910930, p)
     # (0x4880bcf620852a54, 0) RPA point
     # (0, 0x6bed3155c9ada064) RPA point
 
-    infty = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))
-    g = Point(coords, X=gx, Y=gy, Z=Mod(1, p))
+    infty = Point(coords, X=mod(0, p), Y=mod(1, p), Z=mod(0, p))
+    g = Point(coords, X=gx, Y=gy, Z=mod(1, p))
     curve = EllipticCurve(model, coords, p, infty, dict(a=a, b=b))
     return DomainParameters(curve, g, 0x85D265932D90785C, 1)
 
@@ -85,15 +85,15 @@ def test_0y_point(rpa_params):
 @pytest.fixture()
 def distinguish_params(model, coords):
     p = 0xcb5e1d94a6168511
-    a = Mod(0xb166ca7d2dfbf69f, p)
-    b = Mod(0x855bb40cb6937c4b, p)
-    gx = Mod(0x253b2638bd13d6f4, p)
-    gy = Mod(0x1e91a1a182287e71, p)
+    a = mod(0xb166ca7d2dfbf69f, p)
+    b = mod(0x855bb40cb6937c4b, p)
+    gx = mod(0x253b2638bd13d6f4, p)
+    gy = mod(0x1e91a1a182287e71, p)
     # (0x4880bcf620852a54, 0) RPA point
     # (0, 0x6bed3155c9ada064) RPA point
 
-    infty = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))
-    g = Point(coords, X=gx, Y=gy, Z=Mod(1, p))
+    infty = Point(coords, X=mod(0, p), Y=mod(1, p), Z=mod(0, p))
+    g = Point(coords, X=gx, Y=gy, Z=mod(1, p))
     curve = EllipticCurve(model, coords, p, infty, dict(a=a, b=b))
     return DomainParameters(curve, g, 0xcb5e1d94601a3ac5, 1)
 
diff --git a/test/sca/test_zvp.py b/test/sca/test_zvp.py
index 99b7105..f266e55 100644
--- a/test/sca/test_zvp.py
+++ b/test/sca/test_zvp.py
@@ -1,7 +1,7 @@
 import pytest
 
 from pyecsca.ec.coordinates import AffineCoordinateModel
-from pyecsca.ec.mod import Mod
+from pyecsca.ec.mod import Mod, mod
 from pyecsca.ec.mult import LTRMultiplier, AccumulationOrder
 from pyecsca.ec.point import Point
 from pyecsca.sca.re.zvp import (
@@ -233,8 +233,8 @@ def test_zvp(secp128r1, formula):
 def test_points(secp128r1, poly_str, point, k):
     pt = Point(
         AffineCoordinateModel(secp128r1.curve.model),
-        x=Mod(point[0], secp128r1.curve.prime),
-        y=Mod(point[1], secp128r1.curve.prime),
+        x=mod(point[0], secp128r1.curve.prime),
+        y=mod(point[1], secp128r1.curve.prime),
     )
     poly_expr = sympify(poly_str)
     poly = Poly(poly_expr, domain=FF(secp128r1.curve.prime))
diff --git a/test/utils.py b/test/utils.py
index 5d1e80f..1276ad6 100644
--- a/test/utils.py
+++ b/test/utils.py
@@ -29,6 +29,7 @@ class Profiler:
         prof_type: Union[Literal["py"], Literal["c"], Literal["raw"]],
         output_directory: str,
         benchmark_name: str,
+        operations: int = 0
     ):
         self._prof: Union[PyProfiler, cProfiler, RawTimer] = {
             "py": PyProfiler,
@@ -40,6 +41,7 @@ class Profiler:
         self._state = "out"
         self._output_directory = output_directory
         self._benchmark_name = benchmark_name
+        self._operations = operations
 
     def __enter__(self):
         self._prof.__enter__()
@@ -90,7 +92,7 @@ class Profiler:
         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
+            print(f"{self._prof.duration:.4f}s {(self._operations/self._prof.duration) if self._operations else '-':.1f}op/s")  # type: ignore
 
     def get_time(self) -> float:
         if self._state != "out":