29 files changed, 251 insertions, 256 deletions

diff --git a/pyecsca/ec/curve.py b/pyecsca/ec/curve.py
index ff89819..fc48215 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)
@@ -94,7 +94,7 @@ class EllipticCurve:
                 if value.n != prime:
                     raise ValueError(f"Parameter {name} has wrong modulus.")
             else:
-                value = Mod(value, prime)
+                value = mod(value, prime)
             self.parameters[name] = value
         self.neutral = neutral
         self.__validate_coord_assumptions()
@@ -147,9 +147,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 +234,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 +314,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 +322,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..d27b881 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, SymbolicMod, mod
 from pyecsca.ec.op import CodeOp, OpType
 from pyecsca.misc.cfg import getconfig
 from pyecsca.misc.cache import sympify
@@ -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
 
@@ -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/base.py b/pyecsca/ec/mod/base.py
index d059120..ed5bbc2 100644
--- a/pyecsca/ec/mod/base.py
+++ b/pyecsca/ec/mod/base.py
@@ -1,5 +1,6 @@
 import random
 import secrets
+from abc import ABC
 from functools import lru_cache, wraps
 
 from public import public
@@ -123,7 +124,7 @@ _mod_order = ["gmp", "flint", "python"]
 
 
 @public
-class Mod:
+class Mod(ABC):
     """
     An element x of ℤₙ.
 
@@ -134,8 +135,8 @@ class Mod:
 
     Has all the usual special methods that upcast integers automatically:
 
-    >>> a = Mod(3, 5)
-    >>> b = Mod(2, 5)
+    >>> a = mod(3, 5)
+    >>> b = mod(2, 5)
     >>> a + b
     0
     >>> a * 2
@@ -163,21 +164,9 @@ class Mod:
     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
-        )
+    def __init__(self, x, n):
+        self.x = x
+        self.n = n
 
     @_check
     def __add__(self, other) -> "Mod":
@@ -269,7 +258,7 @@ class Mod:
         :return: The random :py:class:`Mod`.
         """
         with RandomModAction(n) as action:
-            return action.exit(cls(secrets.randbelow(n), n))
+            return action.exit(mod(secrets.randbelow(n), n))
 
     def __pow__(self, n) -> "Mod":
         return NotImplemented
@@ -288,8 +277,7 @@ class Undefined(Mod):
         return object.__new__(cls)
 
     def __init__(self):
-        self.x = None
-        self.n = None
+        super().__init__(None, None)
 
     def __add__(self, other):
         return NotImplemented
@@ -359,3 +347,17 @@ class Undefined(Mod):
 
     def __pow__(self, n):
         return NotImplemented
+
+
+@public
+def mod(x, n) -> Mod:
+    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
index f3d8e66..47c3c04 100644
--- a/pyecsca/ec/mod/flint.py
+++ b/pyecsca/ec/mod/flint.py
@@ -56,9 +56,6 @@ if has_flint:
         _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],
diff --git a/pyecsca/ec/mod/gmp.py b/pyecsca/ec/mod/gmp.py
index 1d95189..cdeb3b0 100644
--- a/pyecsca/ec/mod/gmp.py
+++ b/pyecsca/ec/mod/gmp.py
@@ -44,9 +44,6 @@ if has_gmp:
         n: gmpy2.mpz
         __slots__ = ("x", "n")
 
-        def __new__(cls, *args, **kwargs):
-            return object.__new__(cls)
-
         def __init__(
             self,
             x: Union[int, gmpy2.mpz],
diff --git a/pyecsca/ec/mod/raw.py b/pyecsca/ec/mod/raw.py
index 06eda1f..5f4c78c 100644
--- a/pyecsca/ec/mod/raw.py
+++ b/pyecsca/ec/mod/raw.py
@@ -15,9 +15,6 @@ class RawMod(Mod):
     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
@@ -115,6 +112,6 @@ class RawMod(Mod):
         return RawMod(pow(self.x, n, self.n), self.n)
 
 
-from pyecsca.ec.mod.base import _mod_classes
+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
index 25ce91d..0c8a4ba 100644
--- a/pyecsca/ec/mod/symbolic.py
+++ b/pyecsca/ec/mod/symbolic.py
@@ -26,9 +26,6 @@ class SymbolicMod(Mod):
     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
@@ -117,6 +114,6 @@ class SymbolicMod(Mod):
         return self.__class__(pow(self.x, n), self.n)
 
 
-from pyecsca.ec.mod.base import _mod_classes
+from pyecsca.ec.mod.base import _mod_classes  # noqa
 
 _mod_classes["symbolic"] = SymbolicMod
diff --git a/pyecsca/ec/mod/test.pyx b/pyecsca/ec/mod/test.pyx
index 080d950..b8c88d3 100644
--- a/pyecsca/ec/mod/test.pyx
+++ b/pyecsca/ec/mod/test.pyx
@@ -1,7 +1,11 @@
-import cython
+cdef class Test:
+    def __init__(self):
+        print("here")
 
+cdef class SubTest(Test):
+    def __init__(self):
+        print("sub init")
 
-@cython.cclass
-class Test:
-    def __new__(cls, *args, **kwargs):
-        pass
+cdef class OtherTest(Test):
+    def __init__(self):
+        print("other init")
diff --git a/pyecsca/ec/params.py b/pyecsca/ec/params.py
index 6693581..7dfaab8 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,7 +258,7 @@ def _create_params(curve, coords, infty):
                 )
             value = ilocals[coordinate]
             if isinstance(value, int):
-                value = Mod(value, field)
+                value = mod(value, field)
             infinity_coords[coordinate] = value
         infinity = Point(coord_model, **infinity_coords)
     elliptic_curve = EllipticCurve(model, coord_model, field, infinity, params)  # type: ignore[arg-type]
@@ -266,8 +266,8 @@ def _create_params(curve, coords, infty):
         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..ba79ace 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)
                         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/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..f83cc14 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, mod
 from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel
 from pyecsca.ec.point import Point, InfinityPoint
 import pytest
@@ -28,9 +28,9 @@ def test_construction(secp128r1_coords):
 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/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))