Compute formula factor set.

3 files changed, 127 insertions, 32 deletions

diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index 536c599..93aab6f 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -11,35 +11,58 @@ from public import public
 
 from sympy import symbols, FF, Poly, Monomial, Symbol, Expr
 
-from ...ec.context import DefaultContext, local
 from ...ec.curve import EllipticCurve
 from ...ec.divpoly import mult_by_n
 from ...ec.formula import Formula
-from ...ec.mod import SymbolicMod, Mod
+from ...ec.mod import Mod
+from ...ec.params import DomainParameters
 from ...ec.point import Point
-from ...misc.cfg import TemporaryConfig
 
 
 @public
-def unroll_formula(formula: Formula, prime: int) -> List[Poly]:
+def unroll_formula(formula: Formula) -> List[Poly]:
     """
     Unroll a given formula symbolically to obtain symbolic expressions for its intermediate values.
 
     :param formula: Formula to unroll.
-    :param prime: Field to unroll over, necessary for technical reasons.
     :return: List of symbolic intermediate values.
     """
-    field = FF(prime)
-    inputs = [Point(formula.coordinate_model,
-                    **{var: SymbolicMod(symbols(var + str(i)), prime) for var in formula.coordinate_model.variables})
-              for i in
-              range(1, 1 + formula.num_inputs)]
-    params = {var: SymbolicMod(symbols(var), prime) for var in formula.coordinate_model.curve_model.parameter_names}
-    with local(DefaultContext()) as ctx, TemporaryConfig() as cfg:
-        cfg.ec.mod_implementation = "symbolic"
-        formula(prime, *inputs, **params)
-    return [Poly(op_result.value.x, *op_result.value.x.free_symbols, domain=field) for op_result in
-            ctx.actions.get_by_index([0])[0].op_results]
+    params = {var: symbols(var) for var in formula.coordinate_model.curve_model.parameter_names}
+    inputs = {f"{var}{i}": symbols(f"{var}{i}") for var in formula.coordinate_model.variables for i in
+              range(1, formula.num_inputs + 1)}
+    locals = {**params, **inputs}
+    values = []
+    for op in formula.code:
+        result = op(**locals)
+        locals[op.result] = result
+        values.append(result)
+    return [Poly(value) for value in values]
+
+
+def compute_factor_set(formula: Formula) -> Set[Poly]:
+    """
+
+    :param formula:
+    :return:
+    """
+    unrolled = unroll_formula(formula)
+    factors = set()
+    # Go over all of the unrolled intermediates
+    for poly in unrolled:
+        # Factor the intermediate, don't worry about the coeff
+        coeff, factor_list = poly.factor_list()
+        # Go over all the factors of the intermediate, forget the power
+        for factor, power in factor_list:
+            # Remove unnecessary variables from the Poly
+            reduced = factor.exclude()
+            # If there are only lowercase variables remaining, those are only curve parameters
+            # so we do not care about the polynomial
+            if all(str(gen).islower() for gen in reduced.gens):  # type: ignore[attr-defined]
+                continue
+            # Divide out the GCD of the coefficients from the poly
+            _, reduced = reduced.primitive()
+            factors.add(reduced)
+    return factors
 
 
 def curve_equation(x: Symbol, curve: EllipticCurve, symbolic: bool = True) -> Expr:
@@ -65,6 +88,7 @@ def subs_curve_equation(poly: Poly, curve: EllipticCurve) -> Poly:
     :param curve:
     :return:
     """
+    poly = Poly(poly, domain=FF(curve.prime))
     gens = poly.gens  # type: ignore[attr-defined]
     terms = []
     for term in poly.terms():
@@ -90,6 +114,7 @@ def subs_curve_params(poly: Poly, curve: EllipticCurve) -> Poly:
     :param curve:
     :return:
     """
+    poly = Poly(poly, domain=FF(curve.prime))
     for name, value in curve.parameters.items():
         symbol = symbols(name)
         if symbol in poly.gens:  # type: ignore[attr-defined]
@@ -106,6 +131,7 @@ def subs_dlog(poly: Poly, k: int, curve: EllipticCurve):
     :param curve:
     :return:
     """
+    poly = Poly(poly, domain=FF(curve.prime))
     X1, X2 = symbols("X1,X2")
     gens = poly.gens  # type: ignore[attr-defined]
     if X2 not in gens or X1 not in gens:
@@ -118,13 +144,23 @@ def subs_dlog(poly: Poly, k: int, curve: EllipticCurve):
     mx, my = mult_by_n(curve, k)
     u, v = mx[0].subs("x", X1), mx[1].subs("x", X1)
 
+    # The polynomials are quite dense, hence it makes sense
+    # to compute all of the u and v powers in advance and
+    # just use them, because they will likely all be needed.
+    # Note, this has a memory cost...
+    u_powers = [1]
+    v_powers = [1]
+    for i in range(1, max_degree + 1):
+        u_powers.append(u_powers[i - 1] * u)
+        v_powers.append(v_powers[i - 1] * v)
+
     res = 0
     for term in poly.terms():
         powers = list(term[0])
         u_power = powers[X2i]
-        u_factor = u**u_power
+        u_factor = u_powers[u_power]
         v_power = max_degree - u_power
-        v_factor = v**v_power
+        v_factor = v_powers[v_power]
         powers[X2i] = 0
         monom = Monomial(powers, gens).as_expr() * term[1]
         res += Poly(monom, *new_gens, domain=poly.domain) * u_factor * v_factor
@@ -154,6 +190,7 @@ def eliminate_y(poly: Poly, curve: EllipticCurve) -> Poly:
     :param curve:
     :return:
     """
+    poly = Poly(poly, domain=FF(curve.prime))
     Y1, Y2 = symbols("Y1,Y2")
     gens = poly.gens  # type: ignore[attr-defined]
     Y1i = gens.index(Y1) if Y1 in gens else None
@@ -185,11 +222,11 @@ def eliminate_y(poly: Poly, curve: EllipticCurve) -> Poly:
 
 
 @public
-def zvp_point(poly: Poly, curve: EllipticCurve, k: int) -> Set[Point]:
+def zvp_points(poly: Poly, curve: EllipticCurve, k: int) -> Set[Point]:
     """
     Find a set of ZVP points for a given intermediate value and dlog relationship.
 
-    :param poly: The polynomial to zero out, obtained as a result of :py:meth:`.unroll_formula`.
+    :param poly: The polynomial to zero out, obtained as a result of :py:meth:`.unroll_formula` (or its factor).
     :param curve: The curve to compute over.
     :param k: The discrete-log relationship between the two points, i.e. (X2, Y2) = [k](X1, Y1)
     :return: The set of points (X1, Y1).
@@ -197,6 +234,7 @@ def zvp_point(poly: Poly, curve: EllipticCurve, k: int) -> Set[Point]:
     # If input poly is trivial (only in params), abort early
     if not set(symbols("X1,X2,Y1,Y2")).intersection(poly.gens):  # type: ignore[attr-defined]
         return set()
+    poly = Poly(poly, domain=FF(curve.prime))
     # Start with removing all squares of Y1, Y2
     subbed = subs_curve_equation(poly, curve)
     # Remove the Zs by setting them to 1
diff --git a/test/sca/perf_zvp.py b/test/sca/perf_zvp.py
index 2713009..09491cb 100644
--- a/test/sca/perf_zvp.py
+++ b/test/sca/perf_zvp.py
@@ -3,7 +3,7 @@ import click
 
 from pyecsca.ec.mod import has_gmp
 from pyecsca.misc.cfg import TemporaryConfig
-from pyecsca.sca.re.zvp import zvp_point, unroll_formula
+from pyecsca.sca.re.zvp import zvp_points, unroll_formula
 from pyecsca.ec.params import get_params
 from test.utils import Profiler
 
@@ -38,7 +38,7 @@ def main(profiler, mod, operations, directory):
         )
         with Profiler(profiler, directory, f"zvp_p128_{operations}_{mod}"):
             for _ in range(operations):
-                zvp_point(poly, p128.curve, k)
+                zvp_points(poly, p128.curve, k)
 
 
 if __name__ == "__main__":
diff --git a/test/sca/test_zvp.py b/test/sca/test_zvp.py
index 0a625f9..30eade1 100644
--- a/test/sca/test_zvp.py
+++ b/test/sca/test_zvp.py
@@ -4,7 +4,7 @@ from pyecsca.ec.coordinates import AffineCoordinateModel
 from pyecsca.ec.mod import Mod
 from pyecsca.ec.point import Point
 from pyecsca.sca.re.zvp import unroll_formula, subs_curve_equation, remove_z, eliminate_y, subs_dlog, subs_curve_params, \
-    zvp_point
+    zvp_points, compute_factor_set
 from pyecsca.ec.context import local, DefaultContext
 from sympy import symbols, Poly, sympify, FF
 
@@ -14,15 +14,72 @@ def formula(secp128r1, request):
     return secp128r1.curve.coordinate_model.formulas[request.param]
 
 
-def test_unroll(secp128r1, formula):
-    results = unroll_formula(formula, secp128r1.curve.prime)
+def test_unroll(formula):
+    results = unroll_formula(formula)
     assert results is not None
     for res in results:
         assert isinstance(res, Poly)
 
 
+@pytest.mark.xfail(reason="Not removing Zs so far.")
+def test_factor_set(formula):
+    factor_set = compute_factor_set(formula)
+    assert factor_set is not None
+    assert isinstance(factor_set, set)
+
+    expected_factors = {
+        "add-2007-bl": {
+            "Y2",
+            "Y1",
+            "Y1 + Y2",
+            "X2",
+            "X1",
+            "X1 + X2",
+            "Y1^2 + 2*Y1*Y2 + Y2^2 + X1 + X2",
+            "Y1^2 + 2*Y1*Y2 + Y2^2 + 2*X1 + 2*X2",
+            "X1^2 + X1*X2 + X2^2",
+            "X1^2 + X1*X2 + X2^2 + a",
+            "X1^4 + 2*X1^3*X2 + 3*X1^2*X2^2 + 2*X1*X2^3 + X2^4 - X1*Y1^2 - X2*Y1^2 - 2*X1*Y1*Y2 - 2*X2*Y1*Y2 - X1*Y2^2 - X2*Y2^2 + 2*X1^2*a + 2*X1*X2*a + 2*X2^2*a + a^2",
+            "2*X1^4 + 4*X1^3*X2 + 6*X1^2*X2^2 + 4*X1*X2^3 + 2*X2^4 - 3*X1*Y1^2 - 3*X2*Y1^2 - 6*X1*Y1*Y2 - 6*X2*Y1*Y2 - 3*X1*Y2^2 - 3*X2*Y2^2 + 4*X1^2*a + 4*X1*X2*a + 4*X2^2*a + 2*a^2",
+            "2*X1^6 + 6*X1^5*X2 + 12*X1^4*X2^2 + 14*X1^3*X2^3 + 12*X1^2*X2^4 + 6*X1*X2^5 + 2*X2^6 - 3*X1^3*Y1^2 - 6*X1^2*X2*Y1^2 - 6*X1*X2^2*Y1^2 - 3*X2^3*Y1^2 - 6*X1^3*Y1*Y2 - 12*X1^2*X2*Y1*Y2 - 12*X1*X2^2*Y1*Y2 - 6*X2^3*Y1*Y2 - 3*X1^3*Y2^2 - 6*X1^2*X2*Y2^2 - 6*X1*X2^2*Y2^2 - 3*X2^3*Y2^2 + 6*X1^4*a + 12*X1^3*X2*a + 18*X1^2*X2^2*a + 12*X1*X2^3*a + 6*X2^4*a + Y1^4 + 4*Y1^3*Y2 + 6*Y1^2*Y2^2 + 4*Y1*Y2^3 + Y2^4 - 3*X1*Y1^2*a - 3*X2*Y1^2*a - 6*X1*Y1*Y2*a - 6*X2*Y1*Y2*a - 3*X1*Y2^2*a - 3*X2*Y2^2*a + 6*X1^2*a^2 + 6*X1*X2*a^2 + 6*X2^2*a^2 + 2*a^3"
+        },
+        "add-2016-rcb": {
+            "Y2",
+            "Y2 + 1",
+            "Y1",
+            "Y1 + 1",
+            "Y1 + Y2",
+            "X2",
+            "X2 + 1",
+            "X2 + Y2",
+            "X1",
+            "X1 + 1",
+            "X1 + Y1",
+            "X1 + X2",
+            "X1*a + X2*a + 3*b",
+            "-Y1*Y2 + X1*a + X2*a + 3*b",
+            "Y1*Y2 + 1",
+            "Y1*Y2 + X1*a + X2*a + 3*b",
+            "X2*Y1 + X1*Y2",
+            "-X1*X2 + a",
+            "X1*X2 + 1",
+            "X1*X2 + Y1*Y2",
+            "3*X1*X2 + a",
+            "X1*X2*a - a^2 + 3*X1*b + 3*X2*b",
+            "-X2*Y1^2*Y2 - X1*Y1*Y2^2 + 2*X1*X2*Y1*a + X2^2*Y1*a + X1^2*Y2*a + 2*X1*X2*Y2*a - Y1*a^2 - Y2*a^2 + 3*X1*Y1*b + 6*X2*Y1*b + 6*X1*Y2*b + 3*X2*Y2*b",
+            "3*X1*X2^2*Y1 + 3*X1^2*X2*Y2 + Y1^2*Y2 + Y1*Y2^2 + X1*Y1*a + 2*X2*Y1*a + 2*X1*Y2*a + X2*Y2*a + 3*Y1*b + 3*Y2*b",
+            "-3*X1^2*X2^2*a - Y1^2*Y2^2 + X1^2*a^2 + 4*X1*X2*a^2 + X2^2*a^2 - 9*X1^2*X2*b - 9*X1*X2^2*b + a^3 + 3*X1*a*b + 3*X2*a*b + 9*b^2"
+        }
+    }
+    if formula.name in expected_factors:
+        expected_set = set(map(lambda s: Poly(s), expected_factors[formula.name]))
+        print(expected_set)
+        print(len(expected_set), len(factor_set))
+        assert factor_set == expected_set
+
+
 def test_curve_elimination(secp128r1, formula):
-    unrolled = unroll_formula(formula, secp128r1.curve.prime)
+    unrolled = unroll_formula(formula)
     subbed = subs_curve_equation(unrolled[-1], secp128r1.curve)
     assert subbed is not None
     Y1, Y2 = symbols("Y1,Y2")
@@ -35,14 +92,14 @@ def test_curve_elimination(secp128r1, formula):
 
 
 def test_remove_z(secp128r1, formula):
-    unrolled = unroll_formula(formula, secp128r1.curve.prime)
+    unrolled = unroll_formula(formula)
     removed = remove_z(unrolled[-1])
     for gen in removed.gens:
         assert not str(gen).startswith("Z")
 
 
 def test_eliminate_y(secp128r1, formula):
-    unrolled = unroll_formula(formula, secp128r1.curve.prime)
+    unrolled = unroll_formula(formula)
     subbed = subs_curve_equation(unrolled[-1], secp128r1.curve)
     eliminated = eliminate_y(subbed, secp128r1.curve)
     assert eliminated is not None
@@ -54,7 +111,7 @@ def test_eliminate_y(secp128r1, formula):
 
 
 def test_full(secp128r1, formula):
-    unrolled = unroll_formula(formula, secp128r1.curve.prime)
+    unrolled = unroll_formula(formula)
     subbed = subs_curve_equation(unrolled[-1], secp128r1.curve)
     removed = remove_z(subbed)
     eliminated = eliminate_y(removed, secp128r1.curve)
@@ -72,10 +129,10 @@ def test_full(secp128r1, formula):
 
 @pytest.mark.slow
 def test_zvp(secp128r1, formula):
-    unrolled = unroll_formula(formula, secp128r1.curve.prime)
+    unrolled = unroll_formula(formula)
     # Try all intermediates, zvp_point should return empty set if ZVP points do not exist
     for poly in unrolled:
-        points = zvp_point(poly, secp128r1.curve, 5)
+        points = zvp_points(poly, secp128r1.curve, 5)
         assert isinstance(points, set)
 
         # If points are produced, try them all.
@@ -102,5 +159,5 @@ def test_points(secp128r1, poly_str, point, k):
                y=Mod(point[1], secp128r1.curve.prime))
     poly_expr = sympify(poly_str)
     poly = Poly(poly_expr, domain=FF(secp128r1.curve.prime))
-    res = zvp_point(poly, secp128r1.curve, k)
+    res = zvp_points(poly, secp128r1.curve, k)
     assert pt in res