Add more formulas from BouncyCastle.

14 files changed, 293 insertions, 19 deletions

diff --git a/pyecsca/sca/re/structural.py b/pyecsca/sca/re/structural.py
index 4fa8508..1a8a200 100644
--- a/pyecsca/sca/re/structural.py
+++ b/pyecsca/sca/re/structural.py
@@ -43,16 +43,21 @@ def formula_similarity_fuzz(
     output_matches = 0.0
     iv_matches = 0.0
     for _ in range(samples):
-        P = curve.affine_random().to_model(one.coordinate_model, curve)
-        Q = curve.affine_random().to_model(other.coordinate_model, curve)
+        Paff = curve.affine_random()
+        Qaff = curve.affine_random()
+        Raff = curve.affine_add(Paff, Qaff)
+        P = Paff.to_model(one.coordinate_model, curve)
+        Q = Qaff.to_model(one.coordinate_model, curve)
+        R = Raff.to_model(one.coordinate_model, curve)
+        inputs = (P, Q, R)[:one.num_inputs]
         with local(DefaultContext()) as ctx:
-            res_one = one(curve.prime, P, Q, **curve.parameters)
+            res_one = one(curve.prime, *inputs, **curve.parameters)
         action_one = ctx.actions.get_by_index([0])
         ivs_one = set(
             map(attrgetter("value"), sum(action_one[0].intermediates.values(), []))
         )
         with local(DefaultContext()) as ctx:
-            res_other = other(curve.prime, P, Q, **curve.parameters)
+            res_other = other(curve.prime, *inputs, **curve.parameters)
         action_other = ctx.actions.get_by_index([0])
         ivs_other = set(
             map(attrgetter("value"), sum(action_other[0].intermediates.values(), []))
diff --git a/test/data/formulas/add-bc-r1rv76 b/test/data/formulas/add-bc-r1rv76-jac
index c33b915..ad7844b 100644
--- a/test/data/formulas/add-bc-r1rv76
+++ b/test/data/formulas/add-bc-r1rv76-jac
@@ -1 +1,2 @@
 source BouncyCastle r1rv76 https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECPoint.java#L749
+coords jacobian
diff --git a/test/data/formulas/add-bc-r1rv76.op3 b/test/data/formulas/add-bc-r1rv76-jac.op3
index 5e7f521..5e7f521 100644
--- a/test/data/formulas/add-bc-r1rv76.op3
+++ b/test/data/formulas/add-bc-r1rv76-jac.op3
diff --git a/test/data/formulas/add-bc-r1rv76-mod b/test/data/formulas/add-bc-r1rv76-mod
new file mode 100644
index 0000000..87158b5
--- /dev/null
+++ b/test/data/formulas/add-bc-r1rv76-mod
@@ -0,0 +1,2 @@
+source BouncyCastle r1rv76 https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECPoint.java#L749
+coords modified
diff --git a/test/data/formulas/add-bc-r1rv76-mod.op3 b/test/data/formulas/add-bc-r1rv76-mod.op3
new file mode 100644
index 0000000..7adff36
--- /dev/null
+++ b/test/data/formulas/add-bc-r1rv76-mod.op3
@@ -0,0 +1,26 @@
+Z1Squared = Z1^2
+U2 = Z1Squared * X2
+Z1Cubed = Z1Squared * Z1
+S2 = Z1Cubed * Y2
+Z2Squared = Z2^2
+U1 = Z2Squared * X1
+Z2Cubed = Z2Squared * Z2
+S1 = Z2Cubed * Y1
+H = U1 - U2
+R = S1 - S2
+HSquared = H^2
+G = HSquared * H
+V = HSquared * U1
+t0 = 2 * V
+t1 = R^2
+t2 = t1 + G
+X3 = t2 - t0
+t3 = V - X3
+t4 = G * S1
+t5 = t3 * R
+Y3 = t5 - t4
+Z3 = H * Z1
+Z3 = Z3 * Z2
+Z3Squared = Z3^2
+W = Z3Squared^2
+T3 = W * a
diff --git a/test/data/formulas/dbl-bc-r1rv76-jac b/test/data/formulas/dbl-bc-r1rv76-jac
new file mode 100644
index 0000000..8fca349
--- /dev/null
+++ b/test/data/formulas/dbl-bc-r1rv76-jac
@@ -0,0 +1,2 @@
+source BouncyCastle r1rv76 https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECPoint.java#L877
+coords jacobian
diff --git a/test/data/formulas/dbl-bc-r1rv76-jac.op3 b/test/data/formulas/dbl-bc-r1rv76-jac.op3
new file mode 100644
index 0000000..c4ebef2
--- /dev/null
+++ b/test/data/formulas/dbl-bc-r1rv76-jac.op3
@@ -0,0 +1,19 @@
+Y1Squared = Y1^2
+T = Y1Squared^2
+X1Squared = X1^2
+M = 3 * X1Squared
+Z1Squared = Z1^2
+Z1Pow4 = Z1Squared^2
+t0 = Z1Pow4 * a
+M = M + t0
+t1 = X1 * Y1Squared
+S = 4 * t1
+t2 = 2 * S
+t3 = M^2
+X3 = t3 - t2
+t4 = 8 * T
+t5 = S - X3
+t6 = t5 * M
+Y3 = t6 - t4
+Z3 = 2 * Y1
+Z3 = Z3 * Z1
diff --git a/test/data/formulas/dbl-bc-r1rv76-mod b/test/data/formulas/dbl-bc-r1rv76-mod
new file mode 100644
index 0000000..f9b9843
--- /dev/null
+++ b/test/data/formulas/dbl-bc-r1rv76-mod
@@ -0,0 +1,2 @@
+source BouncyCastle r1rv76 https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/ECPoint.java#L1321
+coords modified
diff --git a/test/data/formulas/dbl-bc-r1rv76-mod.op3 b/test/data/formulas/dbl-bc-r1rv76-mod.op3
new file mode 100644
index 0000000..201408e
--- /dev/null
+++ b/test/data/formulas/dbl-bc-r1rv76-mod.op3
@@ -0,0 +1,18 @@
+X1Squared = X1^2
+t0 = 3 * X1Squared
+M = t0 + T1
+_2Y1 = 2 * Y1
+_2Y1Squared = _2Y1 * Y1
+t1 = X1 * _2Y1Squared
+S = 2 * t1
+t2 = M^2
+t3 = 2 * S
+X3 = t2 - t3
+_4T = _2Y1Squared^2
+_8T = 2 * _4T
+t4 = S - X3
+t5 = M * t4
+Y3 = t5 - _8T
+t6 = _8T * T1
+T3 = 2 * t6
+Z3 = _2Y1 * Z1
diff --git a/test/data/formulas/dbl-bc-r1rv76-x25519 b/test/data/formulas/dbl-bc-r1rv76-x25519
new file mode 100644
index 0000000..d5d73fb
--- /dev/null
+++ b/test/data/formulas/dbl-bc-r1rv76-x25519
@@ -0,0 +1,4 @@
+source BouncyCastle r1rv76 https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L73
+parameter a24
+assume a24 = (a+2)/4
+coords xz
diff --git a/test/data/formulas/dbl-bc-r1rv76-x25519.op3 b/test/data/formulas/dbl-bc-r1rv76-x25519.op3
new file mode 100644
index 0000000..86f2e50
--- /dev/null
+++ b/test/data/formulas/dbl-bc-r1rv76-x25519.op3
@@ -0,0 +1,9 @@
+pa = X1 + Z1
+pb = X1 - Z1
+pa = pa^2
+pb = pb^2
+X3 = pa * pb
+pa = pa - pb
+Z3 = pa * a24
+Z3 = Z3 + pb
+Z3 = Z3 * pa
diff --git a/test/data/formulas/ladd-bc-r1rv76-x25519 b/test/data/formulas/ladd-bc-r1rv76-x25519
new file mode 100644
index 0000000..01cc06d
--- /dev/null
+++ b/test/data/formulas/ladd-bc-r1rv76-x25519
@@ -0,0 +1,5 @@
+source BouncyCastle r1rv76 https://github.com/bcgit/bc-java/blob/r1rv76/core/src/main/java/org/bouncycastle/math/ec/rfc7748/X25519.java#L111
+parameter a24
+assume a24 = (a+2)/4
+assume Z1 = 1
+coords xz
diff --git a/test/data/formulas/ladd-bc-r1rv76-x25519.op3 b/test/data/formulas/ladd-bc-r1rv76-x25519.op3
new file mode 100644
index 0000000..4142bff
--- /dev/null
+++ b/test/data/formulas/ladd-bc-r1rv76-x25519.op3
@@ -0,0 +1,18 @@
+t1 = X3 + Z3
+X3 = X3 - Z3
+Z3 = X2 + Z2
+X2 = X2 - Z2
+t1 = t1 * X2
+X3 = X3 * Z3
+Z3 = Z3^3
+X2 = X2^2
+t2 = Z3 - X2
+Z2 = t2 * a24
+Z2 = Z2 + X2
+Z4 = Z2 * t2
+X4 = X2 * Z3
+_X3 = t1 + X3
+Z3 = t1 - X3
+X5 = _X3^2
+Z3 = Z3^2
+Z5 = X1 * Z3
diff --git a/test/sca/test_structural.py b/test/sca/test_structural.py
index db75b12..1fe0288 100644
--- a/test/sca/test_structural.py
+++ b/test/sca/test_structural.py
@@ -1,9 +1,18 @@
+import pytest
 from importlib_resources import files, as_file
 import test.data.formulas
-from pyecsca.ec.formula import AdditionEFDFormula
-from pyecsca.ec.model import ShortWeierstrassModel
+from pyecsca.ec.formula import (
+    AdditionEFDFormula,
+    LadderEFDFormula,
+    DoublingEFDFormula,
+    AdditionFormula,
+    DoublingFormula,
+    LadderFormula,
+)
+from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel
 from pyecsca.ec.params import get_params
 from pyecsca.sca.re.structural import formula_similarity, formula_similarity_fuzz
+import itertools
 
 
 def test_formula_match():
@@ -11,26 +20,180 @@ def test_formula_match():
     coords = model.coordinates["jacobian"]
     secp128r1 = get_params("secg", "secp128r1", "jacobian")
     with as_file(
-        files(test.data.formulas).joinpath("add-bc-r1rv76")
+        files(test.data.formulas).joinpath("dbl-bc-r1rv76-jac")
     ) as meta_path, as_file(
-        files(test.data.formulas).joinpath("add-bc-r1rv76.op3")
+        files(test.data.formulas).joinpath("dbl-bc-r1rv76-jac.op3")
     ) as op3_path:
-        bc_formula = AdditionEFDFormula(meta_path, op3_path, "add-bc-r1rv76", coords)
+        bc_formula = DoublingEFDFormula(
+            meta_path, op3_path, "dbl-bc-r1rv76-jac", coords
+        )
     print()
     for other_name, other_formula in coords.formulas.items():
-        if not other_name.startswith("add"):
+        if not other_name.startswith("dbl"):
             continue
-        print(other_name, "fuzz", formula_similarity_fuzz(other_formula, bc_formula, secp128r1.curve, 1000), "symbolic", formula_similarity(other_formula, bc_formula))
+        print(
+            other_name,
+            "fuzz",
+            formula_similarity_fuzz(other_formula, bc_formula, secp128r1.curve, 1000),
+            "symbolic",
+            formula_similarity(other_formula, bc_formula),
+        )
 
 
-def test_efd_formula_match():
-    model = ShortWeierstrassModel()
-    coords = model.coordinates["modified"]
+def test_formula_match1():
+    model = MontgomeryModel()
+    coords = model.coordinates["xz"]
+    curve25519 = get_params("other", "Curve25519", "xz")
+    with as_file(
+        files(test.data.formulas).joinpath("dbl-bc-r1rv76-x25519")
+    ) as meta_path, as_file(
+        files(test.data.formulas).joinpath("dbl-bc-r1rv76-x25519.op3")
+    ) as op3_path:
+        bc_formula = DoublingEFDFormula(
+            meta_path, op3_path, "dbl-bc-r1rv76-x25519", coords
+        )
     print()
     for other_name, other_formula in coords.formulas.items():
-        if not other_name.startswith("add"):
+        if not other_name.startswith("dbl"):
             continue
-        for one_name, one_formula in coords.formulas.items():
-            if not one_name.startswith("add"):
-                continue
-            print(one_name, other_name, formula_similarity(one_formula, other_formula))
+        print(
+            other_name,
+            "fuzz",
+            formula_similarity_fuzz(other_formula, bc_formula, curve25519.curve, 1000),
+            "symbolic",
+            formula_similarity(other_formula, bc_formula),
+        )
+
+
+def test_efd_formula_match():
+    model = ShortWeierstrassModel()
+    coords = model.coordinates["projective"]
+    secp128r1 = get_params("secg", "secp128r1", "projective")
+    print()
+    adds = list(filter(lambda tup: tup[0].startswith("add"), coords.formulas.items()))
+    for one, other in itertools.combinations_with_replacement(adds, 2):
+        one_name, one_formula = one
+        other_name, other_formula = other
+        print(
+            one_name,
+            other_name,
+            "fuzz",
+            formula_similarity_fuzz(one_formula, other_formula, secp128r1.curve, 1000),
+            "symbolic",
+            formula_similarity(one_formula, other_formula),
+        )
+
+
+@pytest.mark.parametrize(
+    "name,model,coords,param_spec,formula_type",
+    [
+        [
+            "add-bc-r1rv76-jac",
+            ShortWeierstrassModel,
+            "jacobian",
+            ("secg", "secp128r1"),
+            AdditionEFDFormula,
+        ],
+        [
+            "add-bc-r1rv76-mod",
+            ShortWeierstrassModel,
+            "modified",
+            ("secg", "secp128r1"),
+            AdditionEFDFormula,
+        ],
+        [
+            "dbl-bc-r1rv76-jac",
+            ShortWeierstrassModel,
+            "jacobian",
+            ("secg", "secp128r1"),
+            DoublingEFDFormula,
+        ],
+        [
+            "dbl-bc-r1rv76-mod",
+            ShortWeierstrassModel,
+            "modified",
+            ("secg", "secp128r1"),
+            DoublingEFDFormula,
+        ],
+        [
+            "dbl-bc-r1rv76-x25519",
+            MontgomeryModel,
+            "xz",
+            ("other", "Curve25519"),
+            DoublingEFDFormula,
+        ],
+        [
+            "ladd-bc-r1rv76-x25519",
+            MontgomeryModel,
+            "xz",
+            ("other", "Curve25519"),
+            LadderEFDFormula,
+        ],
+    ],
+)
+def test_formula_correctness(name, model, coords, param_spec, formula_type):
+    model = model()
+    coordinate_model = model.coordinates[coords]
+    with as_file(files(test.data.formulas).joinpath(name)) as meta_path, as_file(
+        files(test.data.formulas).joinpath(name + ".op3")
+    ) as op3_path:
+        formula = formula_type(meta_path, op3_path, name, coordinate_model)
+    params = get_params(*param_spec, coords)
+    scale = coordinate_model.formulas.get("z", None)
+    if scale is None:
+        scale = coordinate_model.formulas.get("scale", None)
+    for _ in range(10):
+        Paff = params.curve.affine_random()
+        P2aff = params.curve.affine_double(Paff)
+        Qaff = params.curve.affine_random()
+        Q2aff = params.curve.affine_double(Qaff)
+        Raff = params.curve.affine_add(Paff, Qaff)
+        R2aff = params.curve.affine_double(Raff)
+        QRaff = params.curve.affine_add(Qaff, Raff)
+        P = Paff.to_model(coordinate_model, params.curve)
+        P2 = P2aff.to_model(coordinate_model, params.curve)
+        Q = Qaff.to_model(coordinate_model, params.curve)
+        Q2 = Q2aff.to_model(coordinate_model, params.curve)  # noqa
+        R = Raff.to_model(coordinate_model, params.curve)
+        R2 = R2aff.to_model(coordinate_model, params.curve)  # noqa
+        QR = QRaff.to_model(coordinate_model, params.curve)
+        inputs = (P, Q, R)[: formula.num_inputs]
+        res = formula(params.curve.prime, *inputs, **params.curve.parameters)
+        if issubclass(formula_type, AdditionFormula):
+            try:
+                assert res[0].to_affine() == Raff
+            except NotImplementedError:
+                assert (
+                    scale(params.curve.prime, res[0], **params.curve.parameters)[0] == R
+                )
+        elif issubclass(formula_type, DoublingFormula):
+            try:
+                assert res[0].to_affine() == P2aff
+            except NotImplementedError:
+                assert (
+                    scale(params.curve.prime, res[0], **params.curve.parameters)[0]
+                    == P2
+                )
+        elif issubclass(formula_type, LadderFormula):
+            try:
+                # assert res[0].to_affine() == Q2aff
+                assert res[1].to_affine() == QRaff
+            except NotImplementedError:
+                # print(scale(params.curve.prime, res[0], **params.curve.parameters)[0])
+                # print(scale(params.curve.prime, res[1], **params.curve.parameters)[0])
+                # print(P)
+                # print(Q)
+                # print(R)
+                # print(P2)
+                # print(Q2)
+                # print(R2)
+                # print(QR)
+                # print("------------------------------------")
+                assert (
+                    scale(params.curve.prime, res[1], **params.curve.parameters)[0]
+                    == QR
+                )
+                # assert (
+                #    scale(params.curve.prime, res[0], **params.curve.parameters)[0]
+                #    == Q2
+                # )