Add multiplication-by-n polynomial computation to divpoly.

5 files changed, 282 insertions, 186 deletions

diff --git a/pyecsca/ec/divpoly.py b/pyecsca/ec/divpoly.py
new file mode 100644
index 0000000..2f2b635
--- /dev/null
+++ b/pyecsca/ec/divpoly.py
@@ -0,0 +1,263 @@
+from typing import Tuple, Dict, Union, Set, Mapping
+
+from sympy import symbols, FF, Poly
+import networkx as nx
+
+from .curve import EllipticCurve
+from .mod import Mod
+from .model import ShortWeierstrassModel
+
+
+def values(*ns: int) -> Mapping[int, Tuple[int, ...]]:
+    done: Set[int] = set()
+    vals = {}
+    todo: Set[int] = set()
+    todo.update(ns)
+    while todo:
+        val = todo.pop()
+        if val in done:
+            continue
+        new: Tuple[int, ...] = ()
+        if val == -2:
+            new = (-1,)
+        elif val == -1:
+            pass
+        elif val < 0:
+            raise ValueError(f"bad {val}")
+        elif val in (0, 1, 2, 3):
+            pass
+        elif val == 4:
+            new = (-2, 3)
+        elif val % 2 == 0:
+            m = (val - 2) // 2
+            new = (m + 1, m + 3, m, m - 1, m + 2)
+        else:
+            m = (val - 1) // 2
+            if m % 2 == 0:
+                new = (-2, m + 2, m, m - 1, m + 1)
+            else:
+                new = (m + 2, m, -2, m - 1, m + 1)
+        if new:
+            todo.update(new)
+        vals[val] = new
+        done.add(val)
+    return vals
+
+
+def dep_graph(*ns: int):
+    g = nx.DiGraph()
+    vals = values(*ns)
+    for k, v in vals.items():
+        if v:
+            for e in v:
+                g.add_edge(k, e)
+        else:
+            g.add_node(k)
+    return g, vals
+
+
+def dep_map(*ns: int):
+    g, vals = dep_graph(*ns)
+    current: Set[int] = set()
+    ls = []
+    for vert in nx.lexicographical_topological_sort(g, key=lambda v: -sum(g[v].keys())):
+        if vert in current:
+            current.remove(vert)
+        ls.append((vert, set(current)))
+        current.update(vals[vert])
+    ls.reverse()
+    return ls, vals
+
+
+def a_invariants(curve: EllipticCurve) -> Tuple[Mod, ...]:
+    """
+    Compute the a-invariants of the curve.
+
+    :param curve: The elliptic curve (only ShortWeierstrass model).
+    :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)
+        a4 = curve.parameters["a"]
+        a6 = curve.parameters["b"]
+        return a1, a2, a3, a4, a6
+    else:
+        raise NotImplementedError
+
+
+def b_invariants(curve: EllipticCurve) -> Tuple[Mod, ...]:
+    """
+    Compute the b-invariants of the curve.
+
+    :param curve: The elliptic curve (only ShortWeierstrass model).
+    :return: A tuple of 4 b-invariants (b2, b4, b6, b8).
+    """
+    if isinstance(curve.model, ShortWeierstrassModel):
+        a1, a2, a3, a4, a6 = a_invariants(curve)
+        return (a1 * a1 + 4 * a2,
+                a1 * a3 + 2 * a4,
+                a3 ** 2 + 4 * a6,
+                a1 ** 2 * a6 + 4 * a2 * a6 - a1 * a3 * a4 + a2 * a3 ** 2 - a4 ** 2)
+    else:
+        raise NotImplementedError
+
+
+def divpoly0(curve: EllipticCurve, *ns: int) -> Mapping[int, Poly]:
+    """
+    Basically sagemath's division_polynomial_0 but more clever memory management.
+
+    As sagemath says:
+
+        Return the `n^{th}` torsion (division) polynomial, without
+        the 2-torsion factor if `n` is even, as a polynomial in `x`.
+
+        These are the polynomials `g_n` defined in [MT1991]_, but with
+        the sign flipped for even `n`, so that the leading coefficient is
+        always positive.
+
+    :param curve: The elliptic curve.
+    :param ns: The values to compute the polynomial for.
+    :return:
+    """
+    xs = symbols("x")
+
+    K = FF(curve.prime)
+    Kx = lambda r: Poly(r, xs, domain=K)  # noqa
+
+    x = Kx(xs)
+
+    b2, b4, b6, b8 = map(lambda b: Kx(int(b)), b_invariants(curve))
+    ls, vals = dep_map(*ns)
+
+    mem: Dict[int, Poly] = {}
+    for i, keep in ls:
+        if i == -2:
+            val = mem[-1] ** 2
+        elif i == -1:
+            val = Kx(4) * x ** 3 + b2 * x ** 2 + Kx(2) * b4 * x + b6
+        elif i == 0:
+            val = Kx(0)
+        elif i < 0:
+            raise ValueError("n must be a positive integer (or -1 or -2)")
+        elif i == 1 or i == 2:
+            val = Kx(1)
+        elif i == 3:
+            val = Kx(3) * x ** 4 + b2 * x ** 3 + Kx(3) * b4 * x ** 2 + Kx(3) * b6 * x + b8
+        elif i == 4:
+            val = -mem[-2] + (Kx(6) * x ** 2 + b2 * x + b4) * mem[3]
+        elif i % 2 == 0:
+            m = (i - 2) // 2
+            val = mem[m + 1] * (mem[m + 3] * mem[m] ** 2 - mem[m - 1] * mem[m + 2] ** 2)
+        else:
+            m = (i - 1) // 2
+            if m % 2 == 0:
+                val = mem[-2] * mem[m + 2] * mem[m] ** 3 - mem[m - 1] * mem[m + 1] ** 3
+            else:
+                val = mem[m + 2] * mem[m] ** 3 - mem[-2] * mem[m - 1] * mem[m + 1] ** 3
+        for dl in set(mem.keys()).difference(keep).difference(ns):
+            del mem[dl]
+        mem[i] = val
+
+    return mem
+
+
+def divpoly(curve: EllipticCurve, n: int, two_torsion_multiplicity: int = 2) -> Poly:
+    """
+    Compute the n-th division polynomial.
+
+    :param curve:
+    :param n:
+    :param two_torsion_multiplicity:
+    :return:
+    """
+    f: Poly = divpoly0(curve, n)[n]
+    a1, a2, a3, a4, a6 = a_invariants(curve)
+    xs, ys = symbols("x y")
+    x = Poly(xs, xs, domain=f.domain)
+    y = Poly(ys, ys, domain=f.domain)
+
+    if two_torsion_multiplicity == 0:
+        return f
+    elif two_torsion_multiplicity == 1:
+        if n % 2 == 0:
+            Kxy = lambda r: Poly(r, xs, ys, domain=f.domain)  # noqa
+            return Kxy(f) * (Kxy(2) * y + Kxy(a1) * x + Kxy(a3))
+        else:
+            return f
+    elif two_torsion_multiplicity == 2:
+        if n % 2 == 0:
+            return f * divpoly0(curve, -1)[-1]
+        else:
+            return f
+
+
+def mult_by_n(curve: EllipticCurve, n: int) -> Tuple[Tuple[Poly, Poly], Tuple[Poly, Poly]]:
+    """
+    Compute the multiplication-by-n map on an elliptic curve.
+
+    :param curve: Curve to compute on.
+    :param n: Scalar.
+    :return:
+    """
+    xs, ys = symbols("x y")
+    K = FF(curve.prime)
+    x = Poly(xs, xs, domain=K)
+    y = Poly(ys, ys, domain=K)
+    Kxy = lambda r: Poly(r, xs, ys, domain=K)  # noqa
+
+    if n == 1:
+        return x
+
+    a1, a2, a3, a4, a6 = a_invariants(curve)
+
+    polys = divpoly0(curve, -2, -1, n - 1, n, n + 1, n + 2)
+    mx_denom = polys[n] ** 2
+    if n % 2 == 0:
+        mx_num = x * polys[-1] * polys[n] ** 2 - polys[n - 1] * polys[n + 1]
+        mx_denom *= polys[-1]
+    else:
+        mx_num = x * polys[n] ** 2 - polys[-1] * polys[n - 1] * polys[n + 1]
+
+    # Alternative that makes the denominator monic by dividing the
+    # numerator by the leading coefficient. Sage does this
+    # simplification when asking for multiplication_by_m with the
+    # x-only=True, as then the poly is an univariate object.
+    # lc = K(mx_denom.LC())
+    # mx = (mx_num.quo(lc), mx_denom.monic())
+    mx = (mx_num, mx_denom)
+
+    # The following lines compute
+    # my = ((2*y+a1*x+a3)*mx.derivative(x)/m - a1*mx-a3)/2
+    # just as sage does, but using sympy and step-by-step
+    # tracking the numerator and denominator of the fraction.
+
+    # mx.derivative()
+    mxd_num = mx[1] * mx[0].diff() - mx[0] * mx[1].diff()
+    mxd_denom = mx[1] ** 2
+
+    # mx.derivative()/m
+    mxd_dn_num = mxd_num
+    mxd_dn_denom = mxd_denom * Kxy(n)
+
+    # (2*y+a1*x+a3)*mx.derivative(x)/m
+    mxd_full_num = mxd_dn_num * (Kxy(2) * y + Kxy(a1) * x + Kxy(a3))
+    mxd_full_denom = mxd_dn_denom
+
+    # a1*mx
+    a1mx_num = (Kxy(a1) * mx[0]).quo(Kxy(2))
+    a1mx_denom = mx[1]  # noqa
+
+    # a3
+    a3_num = Kxy(a3) * mx[1]
+    a3_denom = mx[1]  # noqa
+
+    # The mx.derivative part has a different denominator, basically mx[1]^2 * m
+    # so the rest needs to be multiplied by this factor when subtracitng.
+    mxd_fact = mx[1] * n
+
+    my_num = (mxd_full_num - a1mx_num * mxd_fact - a3_num * mxd_fact)
+    my_denom = mxd_full_denom * Kxy(2)
+    my = (my_num, my_denom)
+    return mx, my
diff --git a/pyecsca/sca/attack/leakage_model.py b/pyecsca/sca/attack/leakage_model.py
index 937ec65..f9adcff 100644
--- a/pyecsca/sca/attack/leakage_model.py
+++ b/pyecsca/sca/attack/leakage_model.py
@@ -5,7 +5,7 @@ from typing import Literal, ClassVar
 from numpy.random import default_rng
 from public import public
 
-from pyecsca.sca import Trace
+from ...sca.trace import Trace
 
 if sys.version_info[0] < 3 or sys.version_info[0] == 3 and sys.version_info[1] < 10:
     def hw(i):
diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index 836c18e..22f6f79 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -4,21 +4,25 @@ Provides functionality inspired by the Zero-value point attack.
   Zero-Value Point Attacks on Elliptic Curve Cryptosystem, Toru Akishita & Tsuyoshi Takagi , ISC '03
   `<https://doi.org/10.1007/10958513_17>`_
 """
-from typing import Tuple, Dict, Union, Set
+from typing import List
 
-from sympy import symbols, FF, Poly
-import networkx as nx
+from sympy import symbols
 
-from pyecsca.ec.context import DefaultContext, local
-from pyecsca.ec.curve import EllipticCurve
-from pyecsca.ec.formula import Formula
-from pyecsca.ec.mod import SymbolicMod, Mod
-from pyecsca.ec.model import ShortWeierstrassModel
-from pyecsca.ec.point import Point
-from pyecsca.misc.cfg import TemporaryConfig
+from ...ec.context import DefaultContext, local
+from ...ec.formula import Formula
+from ...ec.mod import SymbolicMod
+from ...ec.point import Point
+from ...misc.cfg import TemporaryConfig
 
 
-def unroll_formula(formula: Formula, prime: int):
+def unroll_formula(formula: Formula, prime: int) -> List[SymbolicMod]:
+    """
+    Unroll a given formula symbolically to obtain symbolic expressions for its intermediate values.
+
+    :param formula: Formula to unroll.
+    :param prime: Field to unroll over.
+    :return: List of symbolic intermediate values.
+    """
     inputs = [Point(formula.coordinate_model,
                     **{var: SymbolicMod(symbols(var + str(i)), prime) for var in formula.coordinate_model.variables})
               for i in
@@ -28,174 +32,3 @@ def unroll_formula(formula: Formula, prime: int):
         cfg.ec.mod_implementation = "symbolic"
         formula(prime, *inputs, **params)
     return [op_result.value for op_result in ctx.actions.get_by_index([0])[0].op_results]
-
-
-def values(*ns: int):
-    done: Set[int] = set()
-    vals = {}
-    todo: Set[int] = set()
-    todo.update(ns)
-    while todo:
-        val = todo.pop()
-        if val in done:
-            continue
-        new: Tuple[int, ...] = ()
-        if val == -2:
-            new = (-1,)
-        elif val == -1:
-            pass
-        elif val <= 0:
-            raise ValueError(f"bad {val}")
-        elif val == 1 or val == 2 or val == 3:
-            pass
-        elif val == 4:
-            new = (-2, 3)
-        elif val % 2 == 0:
-            m = (val - 2) // 2
-            new = (m + 1, m + 3, m, m - 1, m + 2)
-        else:
-            m = (val - 1) // 2
-            if m % 2 == 0:
-                new = (-2, m + 2, m, m - 1, m + 1)
-            else:
-                new = (m + 2, m, -2, m - 1, m + 1)
-        if new:
-            todo.update(new)
-        vals[val] = new
-        done.add(val)
-    return vals
-
-
-def dep_graph(*ns: int):
-    g = nx.DiGraph()
-    vals = values(*ns)
-    for k, v in vals.items():
-        if v:
-            for e in v:
-                g.add_edge(k, e)
-        else:
-            g.add_node(k)
-    return g, vals
-
-
-def dep_map(*ns: int):
-    g, vals = dep_graph(*ns)
-    current: Set[int] = set()
-    ls = []
-    for vert in nx.lexicographical_topological_sort(g, key=lambda v: -sum(g[v].keys())):
-        if vert in current:
-            current.remove(vert)
-        ls.append((vert, set(current)))
-        current.update(vals[vert])
-    ls.reverse()
-    return ls, vals
-
-
-def a_invariants(curve: EllipticCurve) -> Tuple[Mod, ...]:
-    if isinstance(curve.model, ShortWeierstrassModel):
-        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
-    else:
-        raise NotImplementedError
-
-
-def b_invariants(curve: EllipticCurve) -> Tuple[Mod, ...]:
-    if isinstance(curve.model, ShortWeierstrassModel):
-        a1, a2, a3, a4, a6 = a_invariants(curve)
-        return (a1 * a1 + 4 * a2,
-                a1 * a3 + 2 * a4,
-                a3 ** 2 + 4 * a6,
-                a1 ** 2 * a6 + 4 * a2 * a6 - a1 * a3 * a4 + a2 * a3 ** 2 - a4 ** 2)
-    else:
-        raise NotImplementedError
-
-
-def divpoly0(curve: EllipticCurve, *ns: int) -> Union[Poly, Tuple[Poly, ...]]:
-    # Basically sagemath's division_polynomial_0 but more clever memory management
-    # and dependency computation.
-    xs = symbols("x")
-
-    K = FF(curve.prime)
-    Kx = lambda r: Poly(r, xs, domain=K)
-
-    x = Kx(xs)
-
-    b2, b4, b6, b8 = map(lambda b: Kx(int(b)), b_invariants(curve))
-    ls, vals = dep_map(*ns)
-
-    mem: Dict[int, Poly] = {}
-    for i, keep in ls:
-        if i == -2:
-            val = mem[-1] ** 2
-        elif i == -1:
-            val = Kx(4) * x ** 3 + b2 * x ** 2 + Kx(2) * b4 * x + b6
-        elif i <= 0:
-            raise ValueError("n must be a positive integer (or -1 or -2)")
-        elif i == 1 or i == 2:
-            val = Kx(1)
-        elif i == 3:
-            val = Kx(3) * x ** 4 + b2 * x ** 3 + Kx(3) * b4 * x ** 2 + Kx(3) * b6 * x + b8
-        elif i == 4:
-            val = -mem[-2] + (Kx(6) * x ** 2 + b2 * x + b4) * mem[3]
-        elif i % 2 == 0:
-            m = (i - 2) // 2
-            val = mem[m + 1] * (mem[m + 3] * mem[m] ** 2 - mem[m - 1] * mem[m + 2] ** 2)
-        else:
-            m = (i - 1) // 2
-            if m % 2 == 0:
-                val = mem[-2] * mem[m + 2] * mem[m] ** 3 - mem[m - 1] * mem[m + 1] ** 3
-            else:
-                val = mem[m + 2] * mem[m] ** 3 - mem[-2] * mem[m - 1] * mem[m + 1] ** 3
-        for dl in set(mem.keys()).difference(keep).difference(ns):
-            del mem[dl]
-        mem[i] = val
-
-    if len(ns) == 1:
-        return mem[ns[0]]
-    else:
-        return tuple(mem[n] for n in ns)
-
-
-def divpoly(curve: EllipticCurve, n: int, two_torsion_multiplicity: int = 2) -> Poly:
-    f: Poly = divpoly0(curve, n)
-    a1, a2, a3, a4, a6 = a_invariants(curve)
-    xs, ys = symbols("x y")
-    x = Poly(xs, xs, domain=f.domain)
-    y = Poly(ys, ys, domain=f.domain)
-
-    if two_torsion_multiplicity == 0:
-        return f
-    elif two_torsion_multiplicity == 1:
-        if n % 2 == 0:
-            Kxy = lambda r: Poly(r, xs, ys, domain=f.domain)
-            return Kxy(f) * (Kxy(2) * y + Kxy(a1) * x + Kxy(a3))
-        else:
-            return f
-    elif two_torsion_multiplicity == 2:
-        if n % 2 == 0:
-            return f * divpoly0(curve, -1)
-        else:
-            return f
-
-
-def mult_by_n(curve: EllipticCurve, n: int) -> Tuple[Poly, Poly]:
-    xs = symbols("x")
-    K = FF(curve.prime)
-    x = Poly(xs, xs, domain=K)
-
-    if n == 1:
-        return x
-
-    polys = divpoly0(curve, -2, -1, n - 1, n, n + 1)
-    denom = polys[3] ** 2
-    if n % 2 == 0:
-        num = x * polys[1] * polys[3] ** 2 - polys[2] * polys[4]
-        denom *= polys[1]
-    else:
-        num = x * polys[3] ** 2 - polys[1] * polys[2] * polys[4]
-    lc = K(denom.LC())
-    return num.quo(lc), denom.monic()
diff --git a/pyecsca/sca/stacked_traces/combine.py b/pyecsca/sca/stacked_traces/combine.py
index 5849150..cbdfe04 100644
--- a/pyecsca/sca/stacked_traces/combine.py
+++ b/pyecsca/sca/stacked_traces/combine.py
@@ -8,8 +8,8 @@ from math import sqrt
 from public import public
 from typing import Callable, Union, Tuple, cast
 
-from pyecsca.sca.trace.trace import CombinedTrace
-from pyecsca.sca.stacked_traces import StackedTraces
+from ...sca.trace.trace import CombinedTrace
+from ...sca.stacked_traces import StackedTraces
 
 TPB = Union[int, Tuple[int, ...]]
 CudaCTX = Tuple[
diff --git a/pyecsca/sca/stacked_traces/stacked_traces.py b/pyecsca/sca/stacked_traces/stacked_traces.py
index 09169bd..f2c67fc 100644
--- a/pyecsca/sca/stacked_traces/stacked_traces.py
+++ b/pyecsca/sca/stacked_traces/stacked_traces.py
@@ -4,7 +4,7 @@ import numpy as np
 from public import public
 from typing import Any, Mapping, Sequence
 
-from pyecsca.sca.trace_set.base import TraceSet
+from ...sca.trace_set.base import TraceSet
 
 
 @public