8 files changed, 138 insertions, 41 deletions

diff --git a/Makefile b/Makefile
index 4c2d6c4..9ebe608 100644
--- a/Makefile
+++ b/Makefile
@@ -23,7 +23,7 @@ codestyle:
 	flake8 --ignore=E501,F405,F403,F401,E126 pyecsca
 
 doc-coverage:
-    interrogate -vv -nmps pyecsca
+	interrogate -vv -nmps pyecsca
 
 docs:
 	$(MAKE) -C docs apidoc
diff --git a/pyecsca/ec/error.py b/pyecsca/ec/error.py
new file mode 100644
index 0000000..f350d9e
--- /dev/null
+++ b/pyecsca/ec/error.py
@@ -0,0 +1,31 @@
+from public import public
+
+
+@public
+class NonInvertibleError(ArithmeticError):
+    pass
+
+
+@public
+class NonInvertibleWarning(UserWarning):
+    pass
+
+
+@public
+class NonResidueError(ArithmeticError):
+    pass
+
+
+@public
+class NonResidueWarning(UserWarning):
+    pass
+
+
+@public
+class UnsatisfiedAssumptionError(ValueError):
+    pass
+
+
+@public
+class UnsatisfiedAssumptionWarning(UserWarning):
+    pass
diff --git a/pyecsca/ec/formula.py b/pyecsca/ec/formula.py
index dd30605..169f4cb 100644
--- a/pyecsca/ec/formula.py
+++ b/pyecsca/ec/formula.py
@@ -1,12 +1,15 @@
 from abc import ABC, abstractmethod
 from ast import parse, Expression
+from astunparse import unparse
 from itertools import product
-from typing import List, Set, Any, ClassVar, MutableMapping, Tuple, Union
+from typing import List, Set, Any, ClassVar, MutableMapping, Tuple, Union, Dict
 
 from pkg_resources import resource_stream
 from public import public
+from sympy import sympify, FF, symbols, Poly
 
 from .context import ResultAction, getcontext, NullContext
+from .error import UnsatisfiedAssumptionError
 from .mod import Mod
 from .op import CodeOp, OpType
 
@@ -38,11 +41,17 @@ class OpResult(object):
 class FormulaAction(ResultAction):
     """An execution of a formula, on some input points and parameters, with some outputs."""
     formula: "Formula"
+    """The formula that was executed."""
     inputs: MutableMapping[str, Mod]
+    """The input variables (point coordinates and parameters)."""
     input_points: List[Any]
+    """The input points."""
     intermediates: MutableMapping[str, List[OpResult]]
+    """Intermediates computed during execution."""
     outputs: MutableMapping[str, OpResult]
+    """The output variables."""
     output_points: List[Any]
+    """The output points."""
 
     def __init__(self, formula: "Formula", *points: Any, **inputs: Mod):
         super().__init__()
@@ -62,8 +71,8 @@ class FormulaAction(ResultAction):
                 parents.append(self.intermediates[parent][-1])
             elif parent in self.inputs:
                 parents.append(self.inputs[parent])
-        l = self.intermediates.setdefault(op.result, list())
-        l.append(OpResult(op.result, value, op.operator, *parents))
+        li = self.intermediates.setdefault(op.result, list())
+        li.append(OpResult(op.result, value, op.operator, *parents))
 
     def add_result(self, point: Any, **outputs: Mod):
         if isinstance(getcontext(), NullContext):
@@ -79,18 +88,29 @@ class FormulaAction(ResultAction):
         return f"{self.__class__.__name__}({self.formula}, {self.input_points}) = {self.output_points}"
 
 
+@public
 class Formula(ABC):
     """A formula operating on points."""
     name: str
+    """Name of the formula."""
+    shortname: ClassVar[str]
+    """A shortname for the type of the formula."""
     coordinate_model: Any
+    """Coordinate model of the formula."""
     meta: MutableMapping[str, Any]
+    """Meta information about the formula, such as its source."""
     parameters: List[str]
+    """Formula parameters (i.e. new parameters introduced by the formula, like `half = 1/2`)."""
     assumptions: List[Expression]
+    """Assumptions of the formula (e.g. `Z1 == 1` or `2*half == 1`)."""
     code: List[CodeOp]
-    shortname: ClassVar[str]
+    """The collection of ops that constitute the code of the formula."""
     num_inputs: ClassVar[int]
+    """Number of inputs (points) of the formula."""
     num_outputs: ClassVar[int]
+    """Number of outputs (points) of the formula."""
     unified: bool
+    """Whether the formula is specifies that it is unified."""
 
     def __call__(self, *points: Any, **params: Mod) -> Tuple[Any, ...]:
         """
@@ -101,13 +121,42 @@ class Formula(ABC):
         :return: The resulting point(s).
         """
         from .point import Point
+        # Validate number of inputs.
         if len(points) != self.num_inputs:
             raise ValueError(f"Wrong number of inputs for {self}.")
+        # Validate input points and unroll them into input params.
         for i, point in enumerate(points):
             if point.coordinate_model != self.coordinate_model:
                 raise ValueError(f"Wrong coordinate model of point {point}.")
             for coord, value in point.coords.items():
                 params[coord + str(i + 1)] = value
+        # Validate assumptions and compute formula parameters.
+        for assumption in self.assumptions:
+            assumption_string = unparse(assumption)[1:-2]
+            lhs, rhs = assumption_string.split(" == ")
+            if lhs in params:
+                # Handle an assumption check on value of input points.
+                alocals: Dict[str, Union[Mod, int]] = {**params}
+                compiled = compile(assumption, "", mode="eval")
+                holds = eval(compiled, None, alocals)
+                if not holds:
+                    raise UnsatisfiedAssumptionError(f"Unsatisfied assumption in the formula ({assumption_string}).")
+            else:
+                field = int(params[next(iter(params.keys()))].n)  # This is nasty...
+                k = FF(field)
+                expr = sympify(f"{rhs} - {lhs}")
+                for curve_param, value in params.items():
+                    expr = expr.subs(curve_param, k(value))
+                if len(expr.free_symbols) > 1 or (param := str(expr.free_symbols.pop())) not in self.parameters:
+                    raise ValueError(
+                        f"This formula couldn't be executed due to an unsupported asusmption ({assumption_string}).")
+                poly = Poly(expr, symbols(param), domain=k)
+                roots = poly.ground_roots()
+                for root in roots.keys():
+                    params[param] = Mod(int(root), field)
+                    break
+                else:
+                    raise UnsatisfiedAssumptionError(f"Unsatisfied assumption in the formula ({assumption_string}).")
         with FormulaAction(self, *points, **params) as action:
             for op in self.code:
                 op_result = op(**params)
@@ -219,7 +268,7 @@ class EFDFormula(Formula):
                     self.parameters.append(line[10:])
                 elif line.startswith("assume"):
                     self.assumptions.append(
-                            parse(line[7:].replace("=", "==").replace("^", "**"), mode="eval"))
+                        parse(line[7:].replace("=", "==").replace("^", "**"), mode="eval"))
                 elif line.startswith("unified"):
                     self.unified = True
                 line = f.readline().decode("ascii")
@@ -259,7 +308,7 @@ class EFDFormula(Formula):
 
 
 @public
-class AdditionFormula(Formula):
+class AdditionFormula(Formula, ABC):
     shortname = "add"
     num_inputs = 2
     num_outputs = 1
@@ -271,7 +320,7 @@ class AdditionEFDFormula(AdditionFormula, EFDFormula):
 
 
 @public
-class DoublingFormula(Formula):
+class DoublingFormula(Formula, ABC):
     shortname = "dbl"
     num_inputs = 1
     num_outputs = 1
@@ -283,7 +332,7 @@ class DoublingEFDFormula(DoublingFormula, EFDFormula):
 
 
 @public
-class TriplingFormula(Formula):
+class TriplingFormula(Formula, ABC):
     shortname = "tpl"
     num_inputs = 1
     num_outputs = 1
@@ -295,7 +344,7 @@ class TriplingEFDFormula(TriplingFormula, EFDFormula):
 
 
 @public
-class NegationFormula(Formula):
+class NegationFormula(Formula, ABC):
     shortname = "neg"
     num_inputs = 1
     num_outputs = 1
@@ -307,7 +356,7 @@ class NegationEFDFormula(NegationFormula, EFDFormula):
 
 
 @public
-class ScalingFormula(Formula):
+class ScalingFormula(Formula, ABC):
     shortname = "scl"
     num_inputs = 1
     num_outputs = 1
@@ -319,7 +368,7 @@ class ScalingEFDFormula(ScalingFormula, EFDFormula):
 
 
 @public
-class DifferentialAdditionFormula(Formula):
+class DifferentialAdditionFormula(Formula, ABC):
     shortname = "dadd"
     num_inputs = 3
     num_outputs = 1
@@ -331,7 +380,7 @@ class DifferentialAdditionEFDFormula(DifferentialAdditionFormula, EFDFormula):
 
 
 @public
-class LadderFormula(Formula):
+class LadderFormula(Formula, ABC):
     shortname = "ladd"
     num_inputs = 3
     num_outputs = 2
diff --git a/pyecsca/ec/mod.py b/pyecsca/ec/mod.py
index a28949c..3f97637 100644
--- a/pyecsca/ec/mod.py
+++ b/pyecsca/ec/mod.py
@@ -2,7 +2,11 @@ import random
 import secrets
 from functools import wraps, lru_cache
 from abc import ABC, abstractmethod
-from typing import Type
+from public import public
+
+from .error import NonInvertibleError, NonResidueError
+from .context import ResultAction
+
 
 has_gmp = False
 try:
@@ -12,10 +16,6 @@ try:
 except ImportError:
     pass
 
-from public import public
-
-from .context import ResultAction
-
 
 @public
 def gcd(a, b):
@@ -92,16 +92,6 @@ def check(func):
 
 
 @public
-class NonInvertibleError(ArithmeticError):
-    pass
-
-
-@public
-class NonResidueError(ArithmeticError):
-    pass
-
-
-@public
 class RandomModAction(ResultAction):
     """A random sampling from Z_n."""
     order: int
@@ -471,6 +461,7 @@ if has_gmp:
                 return GMPMod(self.x, self.n)
             return GMPMod(gmpy2.powmod(self.x, gmpy2.mpz(n), self.n), self.n)
 
+
     Mod = GMPMod
 else:
     Mod = RawMod
diff --git a/pyecsca/ec/model.py b/pyecsca/ec/model.py
index 4832e3c..bfb9184 100644
--- a/pyecsca/ec/model.py
+++ b/pyecsca/ec/model.py
@@ -8,6 +8,7 @@ from public import public
 from .coordinates import EFDCoordinateModel, CoordinateModel
 
 
+@public
 class CurveModel(object):
     """A model(form) of an elliptic curve."""
     name: str
diff --git a/pyecsca/ec/params.py b/pyecsca/ec/params.py
index f92b2ca..2ec9706 100644
--- a/pyecsca/ec/params.py
+++ b/pyecsca/ec/params.py
@@ -1,5 +1,5 @@
 import json
-from sympy import Poly, PythonFiniteField, symbols, sympify
+from sympy import Poly, FF, symbols, sympify
 from astunparse import unparse
 from io import RawIOBase, BufferedIOBase
 from os.path import join
@@ -11,6 +11,7 @@ from public import public
 
 from .coordinates import AffineCoordinateModel, CoordinateModel
 from .curve import EllipticCurve
+from .error import UnsatisfiedAssumptionError
 from .mod import Mod
 from .model import (CurveModel, ShortWeierstrassModel, MontgomeryModel, EdwardsModel,
                     TwistedEdwardsModel)
@@ -133,17 +134,17 @@ def _create_params(curve, coords, infty):
                 exec(compiled, None, alocals)
                 for param, value in alocals.items():
                     if params[param] != value:
-                        raise ValueError(
+                        raise UnsatisfiedAssumptionError(
                             f"Coordinate model {coord_model} has an unsatisifed assumption on the {param} parameter (= {value}).")
             except NameError:
-                k = PythonFiniteField(field)
+                k = FF(field)
                 assumption_string = unparse(assumption)
                 lhs, rhs = assumption_string.split(" = ")
                 expr = sympify(f"{rhs} - {lhs}")
                 for curve_param, value in params.items():
                     expr = expr.subs(curve_param, k(value))
                 if len(expr.free_symbols) > 1 or (param := str(expr.free_symbols.pop())) not in coord_model.parameters:
-                    raise ValueError(f"This coordinate model couldn't be loaded due to unsupported asusmption ({assumption_string}).")
+                    raise ValueError(f"This coordinate model couldn't be loaded due to an unsupported assumption ({assumption_string}).")
                 poly = Poly(expr, symbols(param), domain=k)
                 roots = poly.ground_roots()
                 for root in roots.keys():
@@ -151,7 +152,7 @@ def _create_params(curve, coords, infty):
                         params[param] = Mod(int(root), field)
                         break
                 else:
-                    raise ValueError(f"Coordinate model {coord_model} has an unsatisifed assumption on the {param} parameter (0 = {expr}).")
+                    raise UnsatisfiedAssumptionError(f"Coordinate model {coord_model} has an unsatisifed assumption on the {param} parameter (0 = {expr}).")
 
     # Construct the point at infinity
     infinity: Point
@@ -238,20 +239,25 @@ def load_params(file: Union[str, Path, BinaryIO], coords: str, infty: bool = Tru
 
     return _create_params(curve, coords, infty)
 
+
 @public
 def get_category(category: str, coords: Union[str, Callable[[str], str]],
                   infty: Union[bool, Callable[[str], bool]] = True) -> DomainParameterCategory:
     """
+    Retrieve a category from the std-curves database at https://github.com/J08nY/std-curves.
 
-    :param category:
-    :param coords:
-    :param infty:
-    :return:
+    :param category: The category to retrieve.
+    :param coords: The name of the coordinate system to use. Can be a callable that takes
+                   as argument the name of the curve and produces the coordinate system to use for that curve.
+    :param infty: Whether to use the special :py:class:InfinityPoint (`True`) or try to use the
+                  point at infinity of the coordinate system. Can be a callable that takes
+                  as argument the name of the curve and returns the infinity option to use for that curve.
+    :return: The category.
     """
     listing = resource_listdir(__name__, "std")
     categories = list(entry for entry in listing if resource_isdir(__name__, join("std", entry)))
     if category not in categories:
-        raise ValueError("Category {} not found.".format(category))
+        raise ValueError(f"Category {category} not found.")
     json_path = join("std", category, "curves.json")
     with resource_stream(__name__, json_path) as f:
         return load_category(f, coords, infty)
@@ -273,7 +279,7 @@ def get_params(category: str, name: str, coords: str, infty: bool = True) -> Dom
     listing = resource_listdir(__name__, "std")
     categories = list(entry for entry in listing if resource_isdir(__name__, join("std", entry)))
     if category not in categories:
-        raise ValueError("Category {} not found.".format(category))
+        raise ValueError(f"Category {category} not found.")
     json_path = join("std", category, "curves.json")
     with resource_stream(__name__, json_path) as f:
         category_json = json.load(f)
@@ -281,6 +287,6 @@ def get_params(category: str, name: str, coords: str, infty: bool = True) -> Dom
         if curve["name"] == name:
             break
     else:
-        raise ValueError("Curve {} not found in category {}.".format(name, category))
+        raise ValueError(f"Curve {name} not found in category {category}.")
 
     return _create_params(curve, coords, infty)
diff --git a/test/ec/test_formula.py b/test/ec/test_formula.py
index d5f8392..ed4e6de 100644
--- a/test/ec/test_formula.py
+++ b/test/ec/test_formula.py
@@ -1,6 +1,8 @@
 from unittest import TestCase
 
+from pyecsca.ec.error import UnsatisfiedAssumptionError
 from pyecsca.ec.params import get_params
+from pyecsca.ec.point import Point
 
 
 class FormulaTests(TestCase):
@@ -9,6 +11,9 @@ class FormulaTests(TestCase):
         self.secp128r1 = get_params("secg", "secp128r1", "projective")
         self.add = self.secp128r1.curve.coordinate_model.formulas["add-2007-bl"]
         self.dbl = self.secp128r1.curve.coordinate_model.formulas["dbl-2007-bl"]
+        self.mdbl = self.secp128r1.curve.coordinate_model.formulas["mdbl-2007-bl"]
+        self.jac_secp128r1 = get_params("secg", "secp128r1", "jacobian")
+        self.jac_dbl = self.jac_secp128r1.curve.coordinate_model.formulas["dbl-1998-hnm"]
 
     def test_wrong_call(self):
         with self.assertRaises(ValueError):
@@ -36,3 +41,16 @@ class FormulaTests(TestCase):
         self.assertEqual(self.add.num_powers, 0)
         self.assertEqual(self.add.num_squarings, 6)
         self.assertEqual(self.add.num_addsubs, 10)
+
+    def test_assumptions(self):
+        res = self.mdbl(self.secp128r1.generator, **self.secp128r1.curve.parameters)
+        self.assertIsNotNone(res)
+
+        coords = {name: value * 5 for name, value in self.secp128r1.generator.coords.items()}
+        other = Point(self.secp128r1.generator.coordinate_model, **coords)
+        with self.assertRaises(UnsatisfiedAssumptionError):
+            self.mdbl(other, **self.secp128r1.curve.parameters)
+
+    def test_parameters(self):
+        res = self.jac_dbl(self.jac_secp128r1.generator, **self.jac_secp128r1.curve.parameters)
+        self.assertIsNotNone(res)
diff --git a/test/ec/test_mod.py b/test/ec/test_mod.py
index 3a14c91..71a5a74 100644
--- a/test/ec/test_mod.py
+++ b/test/ec/test_mod.py
@@ -1,6 +1,7 @@
 from unittest import TestCase
 
-from pyecsca.ec.mod import Mod, gcd, extgcd, Undefined, miller_rabin, NonResidueError, NonInvertibleError
+from pyecsca.ec.mod import Mod, gcd, extgcd, Undefined, miller_rabin
+from pyecsca.ec.error import NonInvertibleError, NonResidueError
 
 
 class ModTests(TestCase):