Add filtering of RPA polynomials

2 files changed, 54 insertions, 22 deletions

diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index 206fc5e..e1d4cbc 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -6,11 +6,11 @@ Provides functionality inspired by the Zero-value point attack.
 
 Implements ZVP point construction from [FFD]_.
 """
-from typing import List, Set
+from typing import List, Set, Tuple
 from public import public
 from astunparse import unparse
 
-from sympy import FF, Poly, Monomial, Symbol, Expr, sympify, symbols
+from sympy import FF, Poly, Monomial, Symbol, Expr, sympify, symbols, div
 
 from ...ec.curve import EllipticCurve
 from ...ec.divpoly import mult_by_n
@@ -21,6 +21,10 @@ from ...ec.point import Point
 
 @public
 def unroll_formula(formula: Formula, affine: bool = False) -> List[Poly]:
+    return [poly for _,poly in unroll_formula_with_names(formula, affine)]
+
+
+def unroll_formula_with_names(formula: Formula, affine: bool = False) -> List[Tuple[str, Poly]]:
     """
     Unroll a given formula symbolically to obtain symbolic expressions for its intermediate values.
 
@@ -68,10 +72,10 @@ def unroll_formula(formula: Formula, affine: bool = False) -> List[Poly]:
     for op in formula.code:
         result = op(**locls)
         locls[op.result] = result
-        values.append(result)
+        values.append((op.result,result))
 
     results = []
-    for value in values:
+    for result_var, value in values:
         if affine:
             expr = value
             for lhs, rhs in subs_map.items():
@@ -80,12 +84,12 @@ def unroll_formula(formula: Formula, affine: bool = False) -> List[Poly]:
                 gens = list(expr.free_symbols)
                 gens.sort(key=str)
                 poly = Poly(expr, *gens)
-                results.append(poly)
+                results.append((result_var, poly))
             else:
                 # Skip if no variables remain (constant poly)
                 continue
         else:
-            results.append(Poly(value))
+            results.append((result_var, Poly(value)))
     return results
 
 
@@ -101,10 +105,10 @@ def compute_factor_set(formula: Formula, affine: bool = False) -> Set[Poly]:
     :param affine: Whether to transform the unrolled polynomials (and thus the resulting factors) into affine form.
     :return: The set of factors present in the formula.
     """
-    unrolled = unroll_formula(formula, affine=affine)
+    unrolled = unroll_formula_with_names(formula, affine=affine)
     factors = set()
     # Go over all of the unrolled intermediates
-    for poly in unrolled:
+    for name, 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
@@ -120,9 +124,37 @@ def compute_factor_set(formula: Formula, affine: bool = False) -> Set[Poly]:
             # Make sure the poly has canonical gens order
             canonical = reduced.reorder()
             factors.add(canonical)
+
+    factors = filter_out_rpa_polynomials(factors, formula, unrolled)
     return factors
 
 
+def filter_out_rpa_polynomials(factor_set: Set[Poly], formula: Formula, unrolled: List[Tuple[str, Poly]]) -> Set[Poly]:
+    """
+    Remove polynomials from factorset that imply RPA
+    """
+    
+    # Find polynomials that define the output variables
+    # We save the latest occurence of output variable in the list of ops
+    output_polynomials = {}
+    for name, polynomial in unrolled:
+        if name in formula.outputs:
+            output_polynomials[name] = polynomial
+
+    filtered_factorset = set()
+    for poly in factor_set:
+        # ignore poly = constant*input_variable
+        if poly.is_monomial and poly.is_linear:
+            continue
+        divisible = False
+        for _, output_poly in output_polynomials.items():
+            if div(output_poly, poly)[1] == 0:
+                divisible = True
+        if not divisible:
+            filtered_factorset.add(poly)
+    return filtered_factorset
+                
+    
 def curve_equation(x: Symbol, curve: EllipticCurve, symbolic: bool = True) -> Expr:
     """
     Get the "ysquared" curve polynomial in :paramref:`~.x` for the :paramref:`~.curve`,
diff --git a/test/sca/test_zvp.py b/test/sca/test_zvp.py
index e0e4cfe..7548963 100644
--- a/test/sca/test_zvp.py
+++ b/test/sca/test_zvp.py
@@ -30,30 +30,30 @@ def test_factor_set(formula):
 
     expected_factors = {
         "add-2007-bl": {
-            "y2",
-            "y1",
-            "y1 + y2",
-            "x2",
-            "x1",
+            #"y2", RPA
+            #"y1", RPA
+            #"y1 + y2", RPA
+            #"x2", RPA
+            #"x1", RPA
             "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",
             "a + x1^2 + x1*x2 + x2^2",
-            "a^2 + 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 + 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", RPA
             "2*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^3 + 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 + 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" RPA
         },
         "add-2015-rcb": {
-            "y2",
+            #"y2", RPA
             "y2 + 1",
-            "y1",
+            #"y1", RPA
             "y1 + 1",
             "y1 + y2",
-            "x2",
+            #"x2", RPA
             "x2 + 1",
             "x2 + y2",
-            "x1",
+            #"x1", RPA
             "x1 + 1",
             "x1 + y1",
             "x1 + x2",
@@ -67,9 +67,9 @@ def test_factor_set(formula):
             "x1*x2 + y1*y2",
             "3*x1*x2 + a",
             "a^2 - x1*x2*a - 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"
+            #"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", RPA
+            #"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", RPA
+            #"-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" RPA
         }
     }
     if formula.name in expected_factors: