8 files changed, 91 insertions, 31 deletions

diff --git a/docs/libraries.rst b/docs/libraries.rst
index 6953a32..9631d1a 100644
--- a/docs/libraries.rst
+++ b/docs/libraries.rst
@@ -52,14 +52,17 @@ KeyGen:
  - `Comb <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/crypto/generators/ECKeyPairGenerator.java#L94>`__ via ``ECKeyPairGenerator.generateKeyPair -> ECKeyPairGenerator.createBasePointMultiplier``.
  - `Jacobian-Modified <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECCurve.java#L676>`__ via ``ECCurve.FP_DEFAULT_COORDS``.
    SECP curves use Jacobian, SECT curves use Lambda-Projective.
- - Formulas unknown.
+ - Formulas unknown: `add-bc-r1rv76-jac <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/add-bc-r1rv76-jac.op3>`__,
+`dbl-bc-r1rv76-jac <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/dbl-bc-r1rv76-jac.op3>`__,
+`add-bc-r1rv76-mod <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/add-bc-r1rv76-mod.op3>`__,
+`dbl-bc-r1rv76-mod <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/dbl-bc-r1rv76-mod.op3>`__
 
 Derive:
  - Short-Weierstrass
  - `GLV if possible, else Window NAF <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECCurve.java#L154>`__ via ``ECDHBasicAgreement.calculateAgreement -> ECPoint.multiply -> ECCurve.getMultiplier -> ECCurve.createDefaultMultiplier``.
  - `Jacobian-Modified <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECCurve.java#L676>`__ via ``ECCurve.FP_DEFAULT_COORDS``.
    SECP curves use Jacobian, SECT curves use Lambda-Projective.
- - Formulas unknown.
+ - Formulas same as KeyGen.
 
 ECDSA
 ^^^^^
@@ -69,7 +72,7 @@ KeyGen:
  - `Comb <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/crypto/generators/ECKeyPairGenerator.java#L94>`__ via ``ECKeyPairGenerator.generateKeyPair -> ECKeyPairGenerator.createBasePointMultiplier``.
  - `Jacobian-Modified <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECCurve.java#L676>`__ via ``ECCurve.FP_DEFAULT_COORDS``.
    SECP curves use Jacobian, SECT curves use Lambda-Projective.
- - Formulas unknown.
+ - Formulas same as KeyGen.
 
 Sign:
  - Short-Weierstrass
@@ -77,14 +80,14 @@ Sign:
    ``ECDSASigner.generateSignature -> ECDSASigner.createBasePointMultiplier``.
  - `Jacobian-Modified <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECCurve.java#L676>`__ via ``ECCurve.FP_DEFAULT_COORDS``.
    SECP curves use Jacobian, SECT curves use Lambda-Projective.
- - Formulas unknown.
+ - Formulas same as KeyGen.
 
 Verify:
  - Short-Weierstrass
  - `Multi-scalar GLV if possible, else multi-scalar Window NAF with Shamir's trick <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECAlgorithms.java#L72>`__ via ``ECDSASigner.verifySignature -> ECAlgorithms.sumOfTwoMultiples``.
  - `Jacobian-Modified <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECCurve.java#L676>`__ via ``ECCurve.FP_DEFAULT_COORDS``.
    SECP curves use Jacobian, SECT curves use Lambda-Projective.
- - Formulas unknown.
+ - Formulas same as KeyGen.
 
 X25519
 ^^^^^^
@@ -101,8 +104,9 @@ Derive:
  - `Ladder <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L93>`__ via
    ``X25519.calculateAgreement -> X25519.scalarMult``.
  - `xz <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L68>`__.
- - `dbl-1987-m-3 <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L73>`__ and
-   some `ladd-1987 <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L111>`__ formula.
+ - Unknown formulas: `ladd-bc-r1rv76-x25519 <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/ladd-bc-r1rv76-x25519.op3>`__,
+`dbl-bc-r1rv76-x25519 <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/dbl-bc-r1rv76-x25519.op3>`__.
+Code: `dbl <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L73>`__ and `ladd <https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L111>`__
 
 Ed25519
 ^^^^^^^
@@ -149,7 +153,8 @@ P-224
  - `Comb <https://github.com/google/boringssl/blob/bfa8369795b7533a222a72b7a1bc928941cd66bf/crypto/fipsmodule/ec/p224-64.c#L995>`__ via ``mul_base -> ec_GFp_nistp224_point_mul_base``.
    `Fixed Window <https://github.com/google/boringssl/blob/bfa8369795b7533a222a72b7a1bc928941cd66bf/crypto/fipsmodule/ec/p224-64.c#L947C13-L947C38>`__ via ``mul -> ec_GFp_nistp224_point_mul``.
  - `Jacobian <https://github.com/google/boringssl/blob/bfa8369795b7533a222a72b7a1bc928941cd66bf/crypto/fipsmodule/ec/p224-64.c#L580>`__,
- - Formulas unknown.
+ - Formulas unknown: `add-boringssl-p224 <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/add-boringssl-p224.op3>`__,
+`dbl-boringssl-p224 <https://github.com/J08nY/pyecsca/blob/master/test/data/formulas/dbl-boringssl-p224.op3>`__.
 
 P-256
 ^^^^^
@@ -1406,7 +1411,7 @@ Verify:
  - Coordinates and formulas same as in KeyGen.
 
 x25519
-------
+^^^^^^
 
 KeyGen:
  - Montgomery
diff --git a/pyecsca/ec/point.py b/pyecsca/ec/point.py
index b39cef7..280d746 100644
--- a/pyecsca/ec/point.py
+++ b/pyecsca/ec/point.py
@@ -5,7 +5,7 @@ from typing import Mapping, TYPE_CHECKING
 from public import public
 
 from .context import ResultAction
-from .coordinates import AffineCoordinateModel, CoordinateModel, EFDCoordinateModel
+from .coordinates import AffineCoordinateModel, CoordinateModel
 from .mod import Mod, Undefined
 from .op import CodeOp
 
@@ -44,7 +44,9 @@ class Point:
 
     def __init__(self, model: CoordinateModel, **coords: Mod):
         if not set(model.variables) == set(coords.keys()):
-            raise ValueError(f"Wrong coordinate values for coordinate model, expected {model.variables} got {coords.keys()}.")
+            raise ValueError(
+                f"Wrong coordinate values for coordinate model, expected {model.variables} got {coords.keys()}."
+            )
         self.coordinate_model = model
         self.coords = coords
         field = None
@@ -82,7 +84,9 @@ class Point:
                     pass
             result_variables = set(map(lambda x: x.result, ops))
             if not result_variables.issuperset(affine_model.variables):
-                raise NotImplementedError(f"Coordinate model does have affine mapping function ({result_variables})")
+                raise NotImplementedError(
+                    f"Coordinate model does have affine mapping function ({result_variables})"
+                )
             result = {}
             locls = {**self.coords}
             for op in ops:
@@ -98,7 +102,10 @@ class Point:
             return action.exit(Point(affine_model, **result))
 
     def to_model(
-        self, coordinate_model: CoordinateModel, curve: "EllipticCurve"
+        self,
+        coordinate_model: CoordinateModel,
+        curve: "EllipticCurve",
+        randomized: bool = False,
     ) -> "Point":
         """Convert an affine point into a given coordinate model, if possible."""
         if not isinstance(self.coordinate_model, AffineCoordinateModel):
@@ -122,10 +129,18 @@ class Point:
                     continue
             result = {}
             for var in coordinate_model.variables:
-                if var in locls:  # Try this first.
-                    result[var] = Mod(locls[var], curve.prime) if not isinstance(locls[var], Mod) else locls[var]
+                if var in locls:
+                    result[var] = (
+                        Mod(locls[var], curve.prime)
+                        if not isinstance(locls[var], Mod)
+                        else locls[var]
+                    )
                 else:
                     raise NotImplementedError
+            if randomized:
+                lmbd = Mod.random(curve.prime)
+                for var, value in result.items():
+                    result[var] = value * lmbd**coordinate_model.homogweights[var]
             return action.exit(Point(coordinate_model, **result))
 
     def equals_affine(self, other: "Point") -> bool:
@@ -183,7 +198,13 @@ class Point:
         return self.coords == other.coords
 
     def __hash__(self):
-        return hash((self.coordinate_model, tuple(self.coords.keys()), tuple(self.coords.values())))
+        return hash(
+            (
+                self.coordinate_model,
+                tuple(self.coords.keys()),
+                tuple(self.coords.values()),
+            )
+        )
 
     def __str__(self):
         args = ", ".join([f"{key}={val}" for key, val in self.coords.items()])
@@ -205,7 +226,8 @@ class InfinityPoint(Point):
         return InfinityPoint(AffineCoordinateModel(self.coordinate_model.curve_model))
 
     def to_model(
-        self, coordinate_model: CoordinateModel, curve: "EllipticCurve"
+        self, coordinate_model: CoordinateModel, curve: "EllipticCurve",
+            randomized: bool = False
     ) -> "InfinityPoint":
         return InfinityPoint(coordinate_model)
 
diff --git a/pyecsca/sca/re/structural.py b/pyecsca/sca/re/structural.py
index 7b2e3ca..c79e604 100644
--- a/pyecsca/sca/re/structural.py
+++ b/pyecsca/sca/re/structural.py
@@ -1,4 +1,5 @@
 """"""
+import warnings
 from typing import Dict
 from public import public
 
@@ -12,7 +13,7 @@ from operator import itemgetter, attrgetter
 @public
 def formula_similarity(one: Formula, other: Formula) -> Dict[str, float]:
     if one.coordinate_model != other.coordinate_model:
-        raise ValueError("Mismatched coordinate model.")
+        warnings.warn("Mismatched coordinate model.")
 
     one_unroll = unroll_formula(one)
     other_unroll = unroll_formula(other)
@@ -41,7 +42,7 @@ def formula_similarity_fuzz(
     one: Formula, other: Formula, curve: EllipticCurve, samples: int = 1000
 ) -> Dict[str, float]:
     if one.coordinate_model != other.coordinate_model:
-        raise ValueError("Mismatched coordinate model.")
+        warnings.warn("Mismatched coordinate model.")
 
     output_matches = 0.0
     iv_matches = 0.0
diff --git a/pyecsca/sca/re/zvp.py b/pyecsca/sca/re/zvp.py
index 7e3d790..4e12549 100644
--- a/pyecsca/sca/re/zvp.py
+++ b/pyecsca/sca/re/zvp.py
@@ -3,7 +3,7 @@ Provides functionality inspired by the Zero-value point attack [ZVP]_.
 
 Implements ZVP point construction from [FFD]_.
 """
-from typing import List, Set, Tuple, Dict
+from typing import List, Set, Tuple, Dict, Union
 from public import public
 from astunparse import unparse
 
@@ -17,7 +17,7 @@ from ...ec.point import Point
 
 
 @public
-def unroll_formula(formula: Formula) -> List[Tuple[str, Poly]]:
+def unroll_formula_expr(formula: Formula) -> List[Tuple[str, Expr]]:
     """
     Unroll a given formula symbolically to obtain symbolic expressions for its intermediate values.
 
@@ -33,6 +33,12 @@ def unroll_formula(formula: Formula) -> List[Tuple[str, Poly]]:
         for var in formula.coordinate_model.variables
         for i in range(1, formula.num_inputs + 1)
     }
+    for coord_assumption in formula.coordinate_model.assumptions:
+        assumption_string = unparse(coord_assumption).strip()
+        lhs, rhs = assumption_string.split(" = ")
+        if lhs in params:
+            expr = sympify(rhs, evaluate=False)
+            params[lhs] = expr
     for assumption_string in formula.assumptions_str:
         lhs, rhs = assumption_string.split(" == ")
         if lhs in formula.parameters:
@@ -43,24 +49,40 @@ def unroll_formula(formula: Formula) -> List[Tuple[str, Poly]]:
             params[lhs] = expr
 
     locls = {**params, **inputs}
-    values: List[Tuple[str, Poly]] = []
+    values: List[Tuple[str, Expr]] = []
     for op in formula.code:
         result: Expr = op(**locls)  # type: ignore
         locls[op.result] = result
+        values.append((op.result, result))
+    return values
+
+
+@public
+def unroll_formula(formula: Formula) -> List[Tuple[str, Poly]]:
+    """
+    Unroll a given formula symbolically to obtain symbolic expressions (as Polynomials) for its intermediate values.
+
+    :param formula: Formula to unroll.
+    :return: List of symbolic intermediate values, with associated variable names.
+    """
+    values = unroll_formula_expr(formula)
+    polys = []
+    for name, result in values:
         if result.free_symbols:
             gens = list(result.free_symbols)
             gens.sort(key=str)
             poly = Poly(result, *gens)
-            values.append((op.result, poly))
+            polys.append((name, poly))
         else:
             # TODO: We cannot create a Poly here, because the result does not have free symbols (i.e. it is a constant)
             pass
-
-    return values
+    return polys
 
 
 @public
-def map_to_affine(formula: Formula, polys: List[Tuple[str, Poly]]) -> List[Tuple[str, Poly]]:
+def map_to_affine(
+    formula: Formula, polys: List[Tuple[str, Poly]]
+) -> List[Tuple[str, Poly]]:
     """
     Map unrolled polynomials of a formula to affine form, using some assumptions along the way (e.g. `Z = 1`).
 
diff --git a/test/data/formulas/add-openssl-z256.op3 b/test/data/formulas/add-openssl-z256.op3
index bc0ccfc..3819546 100644
--- a/test/data/formulas/add-openssl-z256.op3
+++ b/test/data/formulas/add-openssl-z256.op3
@@ -11,7 +11,7 @@ H = U2 - U1
 Rsqr = R^2
 Z3 = H * Z1
 Hsqr = H^2
-Z3 = H * Z2
+Z3 = Z3 * Z2
 Hcub = Hsqr * H
 U2 = U1 * Hsqr
 Hsqr = 2 * U2
diff --git a/test/data/formulas/add-sunec-v21-ed25519 b/test/data/formulas/add-sunec-v21-ed25519
index 078aa39..5814455 100644
--- a/test/data/formulas/add-sunec-v21-ed25519
+++ b/test/data/formulas/add-sunec-v21-ed25519
@@ -1,2 +1,2 @@
 source Java JDK 21 https://github.com/openjdk/jdk/blob/jdk-21-ga/src/jdk.crypto.ec/share/classes/sun/security/ec/ed/Ed25519Operations.java#L147
-coords extended
+coords extended-1
diff --git a/test/data/formulas/dbl-sunec-v21-ed25519 b/test/data/formulas/dbl-sunec-v21-ed25519
index 6814542..f20095f 100644
--- a/test/data/formulas/dbl-sunec-v21-ed25519
+++ b/test/data/formulas/dbl-sunec-v21-ed25519
@@ -1,2 +1,2 @@
 source Java JDK 21 https://github.com/openjdk/jdk/blob/jdk-21-ga/src/jdk.crypto.ec/share/classes/sun/security/ec/ed/Ed25519Operations.java#L184
-coords extended
+coords extended-1
diff --git a/test/sca/test_structural.py b/test/sca/test_structural.py
index be1f45f..0587f7a 100644
--- a/test/sca/test_structural.py
+++ b/test/sca/test_structural.py
@@ -11,8 +11,18 @@ from pyecsca.ec.formula import (
 )
 from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel, TwistedEdwardsModel
 from pyecsca.ec.params import get_params
-from pyecsca.sca.re.structural import formula_similarity, formula_similarity_fuzz
-import itertools
+from pyecsca.sca.re.structural import formula_similarity
+
+
+def test_formula_similarity(secp128r1):
+    add_bl = secp128r1.curve.coordinate_model.formulas["add-2007-bl"]
+    add_rcb = secp128r1.curve.coordinate_model.formulas["add-2015-rcb"]
+    out = formula_similarity(add_bl, add_rcb)
+    assert out["output"] == 0
+    assert out["ivs"] < 0.5
+    out_same = formula_similarity(add_bl, add_bl)
+    assert out_same["output"] == 1
+    assert out_same["ivs"] == 1
 
 
 @pytest.mark.parametrize(