5 files changed, 280 insertions, 6 deletions

diff --git a/pyecsca/ec/mod.py b/pyecsca/ec/mod.py
index 301bd70..3dcb609 100644
--- a/pyecsca/ec/mod.py
+++ b/pyecsca/ec/mod.py
@@ -148,7 +148,7 @@ class Mod:
     n: Any
     __slots__ = ("x", "n")
 
-    def __new__(cls, *args, **kwargs):
+    def __new__(cls, *args, **kwargs) -> "Mod":
         if cls != Mod:
             return cls.__new__(cls, *args, **kwargs)
         if not _mod_classes:
diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index 46d4337..e7e3d22 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -4,12 +4,16 @@ 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
 
-from sympy import symbols, FF, poly
+from sympy import symbols, FF, Poly
+import networkx as nx
 
 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
+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
 
@@ -24,3 +28,153 @@ 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(n: int):
+    done = set()
+    vals = {}
+    todo = set()
+    todo.add(n)
+    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(n: int):
+    g = nx.DiGraph()
+    vals = values(n)
+    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(n: int):
+    g, vals = dep_graph(n)
+    current = set()
+    ls = []
+    for vert in nx.topological_sort(g):
+        current.update(vals[vert])
+        if vert in current:
+            current.remove(vert)
+        ls.append((vert, set(current)))
+    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, n: int) -> 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(n)
+
+    mem: Dict[int, Poly] = {}
+    for i, keep in ls:
+        val = None
+        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):
+            del mem[dl]
+        mem[i] = val
+    return mem[n]
+
+
+def divpoly(curve, n, two_torsion_multiplicity=2):
+    f = 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
diff --git a/pyproject.toml b/pyproject.toml
index c8a704d..8660496 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,7 +39,8 @@
     "datashader",
     "xarray",
     "astunparse",
-    "numba==0.57.1"
+    "numba==0.57.1",
+    "networkx"
   ]
 
 [project.urls]
diff --git a/test/sca/test_rpa.py b/test/sca/test_rpa.py
index e563452..2ab784a 100644
--- a/test/sca/test_rpa.py
+++ b/test/sca/test_rpa.py
@@ -5,7 +5,6 @@ from unittest import TestCase
 from parameterized import parameterized
 
 from pyecsca.ec.context import local
-from pyecsca.ec.coordinates import AffineCoordinateModel
 from pyecsca.ec.model import ShortWeierstrassModel
 from pyecsca.ec.curve import EllipticCurve
 from pyecsca.ec.mod import Mod
diff --git a/test/sca/test_zvp.py b/test/sca/test_zvp.py
index a6ce015..0a23f86 100644
--- a/test/sca/test_zvp.py
+++ b/test/sca/test_zvp.py
@@ -1,11 +1,14 @@
 from unittest import TestCase
+from sympy import FF
 
 from pyecsca.ec.model import ShortWeierstrassModel
-from pyecsca.sca.re.zvp import unroll_formula
+from pyecsca.ec.params import get_params
+from pyecsca.sca.re.zvp import unroll_formula, divpoly0, a_invariants, b_invariants, divpoly
 
 
 class ZVPTests(TestCase):
     def setUp(self):
+        self.secp128r1 = get_params("secg", "secp128r1", "projective")
         self.model = ShortWeierstrassModel()
         self.coords = self.model.coordinates["projective"]
         self.add = self.coords.formulas["add-2007-bl"]
@@ -19,3 +22,120 @@ class ZVPTests(TestCase):
         self.assertIsNotNone(results)
         results = unroll_formula(self.neg, 11)
         self.assertIsNotNone(results)
+
+    def test_ainvs(self):
+        ainvs = a_invariants(self.secp128r1.curve)
+        self.assertSequenceEqual(ainvs, (0,
+                                         0,
+                                         0,
+                                         340282366762482138434845932244680310780,
+                                         308990863222245658030922601041482374867))
+
+    def test_binvs(self):
+        binvs = b_invariants(self.secp128r1.curve)
+        self.assertSequenceEqual(binvs, (0,
+                                         340282366762482138434845932244680310777,
+                                         215116352601536216819152607431888567119,
+                                         340282366762482138434845932244680310774))
+
+    def test_divpoly0(self):
+        # Data from sagemath
+        coeffs = [11,
+                  0,
+                  340282366762482138434845932244680302401,
+                  211962053797180672439257756222135086642,
+                  340282366762482138434845932244678441564,
+                  115415922367823003571854983213102698477,
+                  152803211743444076787231275062278784385,
+                  68540219804769369063918923691867278088,
+                  43207172520353703997069627419519708522,
+                  83208285732019037267730920881743782729,
+                  93286967763556583502947234289842152563,
+                  324950611928652823046744874201355360259,
+                  244242343224213805514200367379671854852,
+                  307096814154284337284845014037169929735,
+                  180946781765592277412990188457219828893,
+                  301253861469456022084288029442105687698,
+                  58053323975526190296189278379252064657,
+                  224437885189054146208302696540070489578,
+                  281987318191429654256483850017931541622,
+                  21449216018131966691124843738286677726,
+                  10958264881628724646042625283328121348,
+                  104868338562600481545003572552335444641,
+                  127205813185570107009206143413997395181,
+                  116865717360861207318274706645935808417,
+                  281460458922812844939222119784601506753,
+                  336607098463310980140968249747513775735,
+                  304486486784143285234063826161805094682,
+                  194935097339732797131694429642153881938,
+                  193523171473792085604518744912658246509,
+                  204844449336357293979832621297234119270,
+                  244481753281744913785581086721299830802,
+                  46816299473081369405217767361380254657,
+                  303070923752707405164354702252828590781,
+                  222516549119176621389776816552836322766,
+                  292006660232236762950883960515487362063,
+                  53617127992846936725441702182362940200,
+                  242498306026562585655027965022211017540,
+                  25039963304689451659955607939868533124,
+                  328580435950647191774558154445103295305,
+                  24226614081978788956695324769468902511,
+                  147945052666123617872720080832548744564,
+                  287190187011075399698210761813202261601,
+                  117131681517270554750959286838283723521,
+                  35018410385280384289320020556813474742,
+                  83939964512240352730304831725346032711,
+                  147219996946006689656600631222993527180,
+                  280430477096741745234510250577626566690,
+                  32753113267385981127807026368593329576,
+                  105134319561523011785486683031223863934,
+                  206456116679151691099661865534540095270,
+                  116180470443213022739312068090342951131,
+                  245850120846480965440408943459023315919,
+                  45805943896736805301879725516256422457,
+                  226777421435695229777151315574975350291,
+                  283680841707610526659029980964566557627,
+                  53168487339451866167506032177471934158,
+                  69212302225932892622760219621519562036,
+                  183916411340675637978873336955593385541,
+                  119478537598919956688656337369481692789,
+                  234767298887335988751880131162396819780,
+                  218412162101425422347176804186940045781]
+
+        K = FF(self.secp128r1.curve.prime)
+        poly = divpoly0(self.secp128r1.curve, 11)
+        computed = list(map(K, poly.all_coeffs()))
+        self.assertListEqual(coeffs, computed)
+
+    def test_divpoly(self):
+        # Data from sagemath
+        K = FF(self.secp128r1.curve.prime)
+        coeffs_0 = {
+            (0, ): K(16020440675387382717114730680672549016),
+            (1, ): K(269851015321770885610377847857290470365),
+            (2, ): K(340282366762482138434845932244680310693),
+            (3, ): K(109469325440469337582450480850803806492),
+            (4, ): K(340282366762482138434845932244680310753),
+            (6, ): K(2)
+        }
+        self.assertDictEqual(divpoly(self.secp128r1.curve, 4, 0).as_dict(), coeffs_0)
+        coeffs_1 = {
+            (6, 1): K(4),
+            (4, 1): K(340282366762482138434845932244680310723),
+            (3, 1): K(218938650880938675164900961701607612984),
+            (2, 1): K(340282366762482138434845932244680310603),
+            (1, 1): K(199419663881059632785909763469900629947),
+            (0, 1): K(32040881350774765434229461361345098032)
+        }
+        self.assertDictEqual(divpoly(self.secp128r1.curve, 4, 1).as_dict(), coeffs_1)
+        coeffs_2 = {
+            (9, ): K(8),
+            (7, ): K(340282366762482138434845932244680310639),
+            (6, ): K(187545273439985507098415273777631738640),
+            (4, ): K(117928913205007755574446043156465405646),
+            (3, ): K(244159722710157842132157548160645018307),
+            (2, ): K(200234655086793134086408617236124137371),
+            (1, ): K(51914434605509249526780779992574428819),
+            (0, ): K(60581150995923875019702403440670701629)
+        }
+        self.assertDictEqual(divpoly(self.secp128r1.curve, 4, 2).as_dict(), coeffs_2)