Return a point deterministically for zer ZVP polynomial.

2 files changed, 53 insertions, 65 deletions

diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index f1c5f3d..c4e1496 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -205,7 +205,9 @@ def is_homogeneous(polynomial: Poly, weighted_variables: List[Tuple[str, int]])
 
 
 @public
-def compute_factor_set(formula: Formula, affine: bool = True, filter_nonhomo: bool = True) -> Set[Poly]:
+def compute_factor_set(
+    formula: Formula, affine: bool = True, filter_nonhomo: bool = True
+) -> Set[Poly]:
     """
     Compute a set of factors present in the :paramref:`~.compute_factor_set.formula`.
 
@@ -505,20 +507,30 @@ def zvp_points(poly: Poly, curve: EllipticCurve, k: int, n: int) -> Set[Point]:
     return points
 
 
+def _deterministic_point_x(curve: EllipticCurve) -> int:
+    x = Mod(1, curve.prime)
+    while True:
+        points = curve.affine_lift_x(x)
+        if points:
+            return int(x)
+        x += 1
+
+
 def solve_easy_dcp(xonly_polynomial: Poly, curve: EllipticCurve) -> Set[Point]:
     points = set()
     final = subs_curve_params(xonly_polynomial, curve)
-    if has_pari:
+    # Solve either via pari or if not available sympy.
+    if final.is_zero:
+        roots = {_deterministic_point_x(curve)}
+    elif final.total_degree() == 0:
+        roots = set()
+    elif has_pari:
         pari = cypari2.Pari()
         polynomial = pari(str(final.expr).replace("**", "^"))
-        if final.degree() == 0:
-            if final.is_zero:
-                # TODO: When the polynomial is zero, sympy does not give any roots, we want to return the generator.
-                pass
-            return set()
-        roots = list(map(int, pari.polrootsmod(polynomial, curve.prime)))
+        roots = set(map(int, pari.polrootsmod(polynomial, curve.prime)))
     else:
         roots = final.ground_roots().keys()
+
     for root in roots:
         points.update(curve.affine_lift_x(Mod(int(root), curve.prime)))
     return points
@@ -526,6 +538,7 @@ def solve_easy_dcp(xonly_polynomial: Poly, curve: EllipticCurve) -> Set[Point]:
 
 def solve_hard_dcp(xonly_polynomial: Poly, curve: EllipticCurve, k: int) -> Set[Point]:
     points = set()
+    # Solve either via pari or if not available sympy.
     if has_pari:
         roots = solve_hard_dcp_cypari(xonly_polynomial, curve, k)
     else:
@@ -533,26 +546,32 @@ def solve_hard_dcp(xonly_polynomial: Poly, curve: EllipticCurve, k: int) -> Set[
         dlog = subs_dlog(xonly_polynomial, k, curve)
         # Put in concrete curve parameters
         final = subs_curve_params(dlog, curve)
-        # Find the roots (X1)
-        roots = final.ground_roots().keys()
-        # Finally lift the roots to find the points (if any)
-        # TODO: When the polynomial is zero, sympy does not give any roots, we want to return the generator.
+        if final.is_zero:
+            roots = {_deterministic_point_x(curve)}
+        else:
+            # Find the roots (X1)
+            roots = final.ground_roots().keys()
+
+    # 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)))
     return points
 
 
-def solve_hard_dcp_cypari(xonly_polynomial: Poly, curve: EllipticCurve, k: int) -> Set[int]:
+def solve_hard_dcp_cypari(
+    xonly_polynomial: Poly, curve: EllipticCurve, k: int
+) -> Set[int]:
     try:
         a, b = int(curve.parameters["a"]), int(curve.parameters["b"])
         xonly_polynomial = subs_curve_params(xonly_polynomial, curve)
+        if xonly_polynomial.is_zero:
+            return {_deterministic_point_x(curve)}
 
         # k^2 * degree
         # k=25, deg=6, 128bit -> 3765, a 20MB
         # k=32, deg=6, 128bit -> 6150, a 32MB
         # k=10, deg=6, 128bit -> 606, a 4MB
-        outdegree = k**2 * xonly_polynomial.degree()
-        # TODO: When can outdegree be 0?
+        outdegree = k**2 * xonly_polynomial.total_degree()
         # Magic heuristic, plus some constant term for very small polys
         stacksize = 2 * (outdegree * (40 * curve.prime.bit_length())) + 1000000
         stacksizemax = 15 * stacksize
@@ -568,14 +587,13 @@ def solve_hard_dcp_cypari(xonly_polynomial: Poly, curve: EllipticCurve, k: int)
         subspoly = 0
         x = pari("x")
         num, den = pari.subst(mul[0], x, x1), pari.subst(mul[1], x, x1)
-        for deg in range(polydeg+1):
+        for deg in range(polydeg + 1):
             monomial = pari.polcoef(polynomial, deg, x2)
-            monomial *= polypower(pari, num, deg)
-            monomial *= polypower(pari, den, polydeg-deg)
+            monomial *= num**deg
+            monomial *= den**(polydeg - deg)
             subspoly += monomial
         if subspoly == pari.zero():
-            # TODO: Return the generator here when the API changes.
-            return set()
+            return {_deterministic_point_x(curve)}
         res = set(map(int, pari.polrootsmod(subspoly, curve.prime)))
     except cypari2.PariError as err:
         raise ValueError("PariError " + err.errtext())
diff --git a/test/sca/test_zvp.py b/test/sca/test_zvp.py
index 8800b11..6413f2d 100644
--- a/test/sca/test_zvp.py
+++ b/test/sca/test_zvp.py
@@ -1,11 +1,8 @@
 import pytest
 
 from pyecsca.ec.coordinates import AffineCoordinateModel
-from pyecsca.ec.curve import EllipticCurve
 from pyecsca.ec.mod import Mod
-from pyecsca.ec.model import ShortWeierstrassModel
 from pyecsca.ec.mult import LTRMultiplier, AccumulationOrder
-from pyecsca.ec.params import DomainParameters
 from pyecsca.ec.point import Point
 from pyecsca.sca.re.zvp import (
     unroll_formula,
@@ -18,7 +15,8 @@ from pyecsca.sca.re.zvp import (
     zvp_points,
     compute_factor_set,
     addition_chain,
-    precomp_zvp_points, solve_easy_dcp,
+    solve_easy_dcp,
+    solve_hard_dcp,
 )
 from pyecsca.ec.context import local, DefaultContext
 from sympy import symbols, Poly, sympify, FF
@@ -254,46 +252,6 @@ def test_addition_chain(secp128r1):
     assert len(res) == 25
 
 
-@pytest.fixture()
-def small_params():
-    model = ShortWeierstrassModel()
-    coords = model.coordinates["projective"]
-    p = 0xC50DE883F0E7B167
-    a = Mod(0x4833D7AA73FA6694, p)
-    b = Mod(0xA6C44A61C5323F6A, p)
-    gx = Mod(0x5FD1F7D38D4F2333, p)
-    gy = Mod(0x21F43957D7E20CEB, p)
-    n = 0xC50DE885003B80EB
-    h = 1
-    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))
-    params = DomainParameters(curve, g, n, h)
-    return params
-
-
-def test_precomp(small_params):
-    chain = addition_chain(
-        123456789,
-        small_params,
-        LTRMultiplier,
-        lambda add, dbl: LTRMultiplier(
-            add, dbl, None, False, AccumulationOrder.PeqRP, True, True
-        ),
-    )
-    res = precomp_zvp_points(
-        chain,
-        {
-            "add": small_params.curve.coordinate_model.formulas["add-2007-bl"],
-            "dbl": small_params.curve.coordinate_model.formulas["dbl-2007-bl"],
-        },
-        small_params,
-    )
-    assert res is not None
-    assert len(res) == len(chain)
-
-
 @pytest.mark.parametrize("k", [7, 25, 31])
 def test_big_boy(secp128r1, k):
     poly_expr = sympify("x1*x2 + y1*y2")
@@ -302,9 +260,21 @@ def test_big_boy(secp128r1, k):
     assert res is not None
 
 
-@pytest.mark.parametrize("numero", [1, 0])
+@pytest.mark.parametrize("numero", [5, 1, 0])
 def test_small_boy(secp128r1, numero):
     x = symbols("x1")
     poly = Poly(numero, x, domain=FF(secp128r1.curve.prime))
     res = solve_easy_dcp(poly, secp128r1.curve)
     assert res is not None
+    if numero == 0:
+        assert len(res) >= 1
+
+
+# secp128r1 has a = -3 so the first poly is zero as well
+@pytest.mark.parametrize("poly_expr", ["a+3", "0"])
+def test_zero_boy(secp128r1, poly_expr):
+    x, a = symbols("x1 a")
+    poly = Poly(sympify(poly_expr), x, a, domain=FF(secp128r1.curve.prime))
+    res = solve_hard_dcp(poly, secp128r1.curve, 5)
+    assert res is not None
+    assert len(res) >= 1