Shamir-Straus trick for pair and quartet

Signed-off-by: Uncle Fatso <uncle.fatso@ghostchain.io>

README.md

@@ -38,29 +38,39 @@ For full background and protocol details see the [project wiki](https://git.ghos
 +==================================================================================================+
 | Deployment Cost                        | Deployment Size |         |         |         |         |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| 2011993                                | 9094            |         |         |         |         |
+| 2624055                                | 11923           |         |         |         |         |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
 |                                        |                 |         |         |         |         |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
 | Function Name                          | Min             | Avg     | Median  | Max     | # Calls |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| addJacobian                            | 1800            | 1800    | 1800    | 1800    | 88      |
+| addJacobian                            | 1896            | 1896    | 1896    | 1896    | 88      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| addProjective                          | 1015            | 1015    | 1015    | 1015    | 44      |
+| addProjective                          | 1110            | 1110    | 1110    | 1110    | 44      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| addProjectiveMixed                     | 967             | 967     | 967     | 967     | 44      |
+| addProjectiveMixed                     | 1084            | 1084    | 1084    | 1084    | 44      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| doubleJacobian                         | 794             | 794     | 794     | 794     | 45      |
+| doubleJacobian                         | 889             | 889     | 889     | 889     | 45      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| doubleProjective                       | 546             | 546     | 546     | 546     | 45      |
+| doubleProjective                       | 641             | 641     | 641     | 641     | 45      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| mulEcTriplet                           | 205548          | 1140310 | 1546125 | 1958623 | 43      |
+| isOnCurve                              | 280             | 280     | 280     | 280     | 262     |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| mulProjectiveTriplet                   | 6436            | 193298  | 334976  | 385511  | 43      |
+| mulEcPair                              | 103258          | 731053  | 1058263 | 1269945 | 44      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| toAffineJacobian                       | 40275           | 47740   | 47859   | 53863   | 133     |
+| mulEcQuartet                           | 308151          | 1561360 | 1982690 | 2672830 | 42      |
 |----------------------------------------+-----------------+---------+---------+---------+---------|
-| toAffineProjective                     | 11056           | 13787   | 13904   | 16105   | 176     |
+| mulEcTriplet                           | 205666          | 1140527 | 1546387 | 1958933 | 43      |
+|----------------------------------------+-----------------+---------+---------+---------+---------|
+| mulProjectivePair                      | 4265            | 178529  | 311013  | 338160  | 44      |
+|----------------------------------------+-----------------+---------+---------+---------+---------|
+| mulProjectiveQuartet                   | 12001           | 206222  | 345559  | 396103  | 42      |
+|----------------------------------------+-----------------+---------+---------+---------+---------|
+| mulProjectiveTriplet                   | 6622            | 199854  | 346604  | 397492  | 43      |
+|----------------------------------------+-----------------+---------+---------+---------+---------|
+| toAffineJacobian                       | 40300           | 47765   | 47884   | 53888   | 133     |
+|----------------------------------------+-----------------+---------+---------+---------+---------|
+| toAffineProjective                     | 11129           | 13876   | 13864   | 16178   | 262     |
 ╰----------------------------------------+-----------------+---------+---------+---------+---------╯
 ```
 

src/MathTester.sol

@@ -81,6 +81,63 @@ contract MathTester {
         return EllipticCurveProjective.mulAddProjectiveTriplet(x1, y1, k1, x2, y2, k2, x3, y3, k3);
     }
 
+    function mulEcPair(
+        uint256 x1,
+        uint256 y1,
+        uint256 k1,
+        uint256 x2,
+        uint256 y2,
+        uint256 k2
+    ) public pure returns(uint256, uint256) {
+        (x1, y1) = EllipticCurve.ecMul(k1, x1, y1, A, P);
+        (x2, y2) = EllipticCurve.ecMul(k2, x2, y2, A, P);
+        (x1, y1) = EllipticCurve.ecAdd(x1, y1, x2, y2, A, P);
+        return (x1, y1);
+    }
+
+    function mulProjectivePair(
+        uint256 x1,
+        uint256 y1,
+        uint256 k1,
+        uint256 x2,
+        uint256 y2,
+        uint256 k2
+    ) public pure returns(uint256, uint256, uint256) {
+        return EllipticCurveProjective.mulAddProjectivePair(x1, y1, k1, x2, y2, k2);
+    }
+
+    function mulEcQuartet(
+        uint256 x1, uint256 y1, uint256 k1,
+        uint256 x2, uint256 y2, uint256 k2,
+        uint256 x3, uint256 y3, uint256 k3,
+        uint256 x4, uint256 y4, uint256 k4
+    ) public pure returns(uint256, uint256) {
+        (x1, y1) = EllipticCurve.ecMul(k1, x1, y1, A, P);
+        (x2, y2) = EllipticCurve.ecMul(k2, x2, y2, A, P);
+        (x3, y3) = EllipticCurve.ecMul(k3, x3, y3, A, P);
+        (x4, y4) = EllipticCurve.ecMul(k4, x4, y4, A, P);
+
+        (x1, y1) = EllipticCurve.ecAdd(x1, y1, x2, y2, A, P);
+        (x3, y3) = EllipticCurve.ecAdd(x3, y3, x4, y4, A, P);
+        (x1, y1) = EllipticCurve.ecAdd(x1, y1, x3, y3, A, P);
+
+        return (x1, y1);
+    }
+
+    function mulProjectiveQuartet(
+        uint256 x1, uint256 y1, uint256 k1,
+        uint256 x2, uint256 y2, uint256 k2,
+        uint256 x3, uint256 y3, uint256 k3,
+        uint256 x4, uint256 y4, uint256 k4
+    ) public pure returns (uint256, uint256, uint256) {
+        return EllipticCurveProjective.mulAddProjectiveQuartet(
+            x1, y1, k1,
+            x2, y2, k2,
+            x3, y3, k3,
+            x4, y4, k4
+        );
+    }
+
     function toAffineJacobian(uint256 x, uint256 y, uint256 z) public pure returns (uint256, uint256) {
         return EllipticCurve.toAffine(x, y, z, P);
     }

src/libraries/ECMathProjective.sol

@@ -10,23 +10,24 @@ pragma solidity ^0.8.0;
 
 library EllipticCurveProjective {
     // Constants are taken from https://en.bitcoin.it/wiki/Secp256k1
-    uint256 private constant A = 0;
+    uint256 public constant B = 7;
-    uint256 private constant B = 7;
+    uint256 public constant P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F;
-    uint256 private constant P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F;
+    uint256 public constant N = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141;
-    uint256 private constant N = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141;
 
     uint256 private constant B3 = 21;
 
     function findMaxBitLength(
         uint256 k1,
         uint256 k2,
-        uint256 k3
+        uint256 k3,
+        uint256 k4
     ) internal pure returns (uint256 bits) {
         assembly {
             // Find maximum of the three scalars
             let max := k1
             if gt(k2, max) { max := k2 }
             if gt(k3, max) { max := k3 }
+            if gt(k4, max) { max := k4 }
 
             // if (v >> 128 != 0) { v >>= 128; bits += 128; }
             if gt(max, 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF) {
@@ -71,6 +72,61 @@ library EllipticCurveProjective {
         }
     }
 
+    function mulAddProjectivePair(
+        uint256 x1, uint256 y1, uint256 k1,
+        uint256 x2, uint256 y2, uint256 k2
+    ) internal pure returns (uint256 x3, uint256 y3, uint256 z3) {
+        // We implement the Straus-Shamir trick described in
+        // Trading Inversions for Multiplications in Elliptic Curve Cryptography.
+        // (https://eprint.iacr.org/2003/257.pdf Page 7).
+
+        uint256 bits = findMaxBitLength(k1, k2, 0, 0);
+
+        uint256[4] memory precomputedXs;
+        uint256[4] memory precomputedYs;
+        uint256[4] memory precomputedZs;
+
+        precomputedXs[1] = x2; precomputedYs[1] = y2; precomputedZs[1] = 1; // 01: P2
+        precomputedXs[2] = x1; precomputedYs[2] = y1; precomputedZs[2] = 1; // 10: P1
+
+        (x3, y3, z3) = projectiveAddMixed(x1, y1, 1, x2, y2); // 11: P1+P2
+        precomputedXs[3] = x3;
+        precomputedYs[3] = y3;
+        precomputedZs[3] = z3;
+
+        x3 = 0;
+        y3 = 1;
+        z3 = 0;
+
+        for (; bits > 0;) {
+            unchecked { --bits; }
+
+            (x3, y3, z3) = projectiveDouble(x3, y3, z3);
+
+            uint8 mask;
+            assembly {
+                mask := or(
+                    shl(1, and(shr(bits, k1), 1)),
+                    and(shr(bits, k2), 1)
+                )
+            }
+
+            if (mask == 0) {
+                continue;
+            }
+
+            if (mask == 3) {
+                (x3, y3, z3) = projectiveAdd(
+                    x3, y3, z3, precomputedXs[mask], precomputedYs[mask], precomputedZs[mask]
+                );
+            } else {
+                (x3, y3, z3) = projectiveAddMixed(
+                    x3, y3, z3, precomputedXs[mask], precomputedYs[mask]
+                );
+            }
+        }
+    }
+
     function mulAddProjectiveTriplet(
         uint256 x1, uint256 y1, uint256 k1,
         uint256 x2, uint256 y2, uint256 k2,
@@ -80,7 +136,7 @@ library EllipticCurveProjective {
         // Trading Inversions for Multiplications in Elliptic Curve Cryptography.
         // (https://eprint.iacr.org/2003/257.pdf Page 7).
 
-        uint256 bits = findMaxBitLength(k1, k2, k3);
+        uint256 bits = findMaxBitLength(k1, k2, k3, 0);
 
         uint256[8] memory precomputedXs;
         uint256[8] memory precomputedYs;
@@ -143,6 +199,118 @@ library EllipticCurveProjective {
         }
     }
 
+    function mulAddProjectiveQuartet(
+        uint256 x1, uint256 y1, uint256 k1,
+        uint256 x2, uint256 y2, uint256 k2,
+        uint256 x3, uint256 y3, uint256 k3,
+        uint256 x4, uint256 y4, uint256 k4
+    ) internal pure returns (uint256 x5, uint256 y5, uint256 z5) {
+        // We implement the Straus-Shamir trick described in
+        // Trading Inversions for Multiplications in Elliptic Curve Cryptography.
+        // (https://eprint.iacr.org/2003/257.pdf Page 7).
+
+        uint256 bits = findMaxBitLength(k1, k2, k3, k4);
+
+        uint256[16] memory precomputedXs;
+        uint256[16] memory precomputedYs;
+        uint256[16] memory precomputedZs;
+
+        precomputedXs[1] = x4; precomputedYs[1] = y4; precomputedZs[1] = 1; // 0001: P4
+        precomputedXs[2] = x3; precomputedYs[2] = y3; precomputedZs[2] = 1; // 0010: P3
+        precomputedXs[4] = x2; precomputedYs[4] = y2; precomputedZs[4] = 1; // 0100: P2
+        precomputedXs[8] = x1; precomputedYs[8] = y1; precomputedZs[8] = 1; // 1000: P1
+
+        (precomputedXs[3], precomputedYs[3], precomputedZs[3]) = projectiveAddMixed(x4, y4, 1, x3, y3); // 0011: P4+P3
+        (precomputedXs[6], precomputedYs[6], precomputedZs[6]) = projectiveAddMixed(x3, y3, 1, x2, y2); // 0110: P3+P2
+        (precomputedXs[12], precomputedYs[12], precomputedZs[12]) = projectiveAddMixed(x2, y2, 1, x1, y1); // 1100: P2+P1
+        (precomputedXs[9], precomputedYs[9], precomputedZs[9]) = projectiveAddMixed(x4, y4, 1, x1, y1); // 1001: P4+P1
+        (precomputedXs[10], precomputedYs[10], precomputedZs[10]) = projectiveAddMixed(x1, y1, 1, x3, y3); // 1010: P1+P3
+        (precomputedXs[5], precomputedYs[5], precomputedZs[5]) = projectiveAddMixed(x4, y4, 1, x2, y2); // 0101: P4+P2
+
+        (precomputedXs[7], precomputedYs[7], precomputedZs[7]) = projectiveAddMixed(
+            precomputedXs[3],
+            precomputedYs[3],
+            precomputedZs[3],
+            x2,
+            y2
+        ); // 0111: P4+P3+P2
+
+        (precomputedXs[11], precomputedYs[11], precomputedZs[11]) = projectiveAddMixed(
+            precomputedXs[3],
+            precomputedYs[3],
+            precomputedZs[3],
+            x1,
+            y1
+        ); // 1011: P4+P3+P1
+
+        (precomputedXs[13], precomputedYs[13], precomputedZs[13]) = projectiveAddMixed(
+            precomputedXs[12],
+            precomputedYs[12],
+            precomputedZs[12],
+            x4,
+            y4
+        ); // 1101: P4+P2+P1
+
+        (precomputedXs[14], precomputedYs[14], precomputedZs[14]) = projectiveAddMixed(
+            precomputedXs[12],
+            precomputedYs[12],
+            precomputedZs[12],
+            x3,
+            y3
+        ); // 1110: P3+P2+P1
+
+        (precomputedXs[15], precomputedYs[15], precomputedZs[15]) = projectiveAddMixed(
+            precomputedXs[14],
+            precomputedYs[14],
+            precomputedZs[14],
+            x4,
+            y4
+        ); // 1111: P4+P3+P2+P1
+
+        x5 = 0;
+        y5 = 1;
+        z5 = 0;
+
+        for (; bits > 0;) {
+            unchecked { --bits; }
+
+            (x5, y5, z5) = projectiveDouble(x5, y5, z5);
+
+            uint8 mask;
+            uint16 bitmask;
+
+            assembly {
+                mask := or(
+                    shl(3, and(shr(bits, k1), 1)),
+                    or(
+                        shl(2, and(shr(bits, k2), 1)),
+                        or(
+                            shl(1, and(shr(bits, k3), 1)),
+                            and(shr(bits, k4), 1)
+                        )
+                    )
+                )
+
+                // bitmask for 1,2,4, 8 is 278 = 0b0000000100110110
+                bitmask := and(shl(278, mask), 1)
+            }
+
+            if (mask == 0) {
+                continue;
+            }
+
+            if (bitmask == 1) {
+                (x5, y5, z5) = projectiveAddMixed(
+                    x5, y5, z5, precomputedXs[mask], precomputedYs[mask]
+                );
+            } else {
+                (x5, y5, z5) = projectiveAdd(
+                    x5, y5, z5, precomputedXs[mask], precomputedYs[mask], precomputedZs[mask]
+                );
+            }
+        }
+    }
+
     function projectiveAddMixed(
         uint256 X1,
         uint256 Y1,
@@ -197,8 +365,6 @@ library EllipticCurveProjective {
         // Y3 = (Y1Y2 + 3bZ1Z2)(Y1Y2 − 3bZ1Z2) + 9bX1X2(X1Z2 + X2Z1),
         // Z3 = (Y1Z2 + Y2Z1)(Y1Y2 + 3bZ1Z2) + 3X1X2(X1Y2 + X2Y1),
 
-        // TODO: Can we optimize it even more?
-
         uint256 t0 = mulmod(X1, X2, P);           // 1.  t0 ← X1 · X2  => (X1·X2)
         uint256 t1 = mulmod(Y1, Y2, P);           // 2.  t1 ← Y1 · Y2  => (Y1·Y2)
         uint256 t2 = mulmod(Z1, Z2, P);           // 3.  t2 ← Z1 · Z2  => (Z1·Z2)
@@ -292,4 +458,14 @@ library EllipticCurveProjective {
         _x = mulmod(x, t0, P);
         _y = mulmod(y, t0, P);
     }
+
+    function isOnCurve(uint256 x, uint256 y) internal pure returns (bool) {
+        // TODO: check if something is zero
+
+        uint lhs = mulmod(y, y, P);               // y^2
+        uint rhs = mulmod(mulmod(x, x, P), x, P); // x^3
+        rhs = addmod(rhs, B, P);                  // x^3 + 7
+
+        return lhs == rhs;
+    }
 }

test/MathTester.t.sol

@@ -22,6 +22,78 @@ contract MathTesterTest is Test {
         points = abi.decode(data, (Point[]));
     }
 
+    function test_quartet() public view {
+        uint256 len = points.length - 3;
+        for (uint256 i; i < len;) {
+            (uint256 x_p, uint256 y_p, uint256 z_p) = math.mulProjectiveQuartet(
+                uint256(points[i].x),
+                uint256(points[i].y),
+                uint256(points[i].k),
+                uint256(points[i+1].x),
+                uint256(points[i+1].y),
+                uint256(points[i+1].k),
+                uint256(points[i+2].x),
+                uint256(points[i+2].y),
+                uint256(points[i+2].k),
+                uint256(points[i+3].x),
+                uint256(points[i+3].y),
+                uint256(points[i+3].k)
+            );
+            (x_p, y_p) = math.toAffineProjective(x_p, y_p, z_p);
+
+            (uint256 x_j, uint256 y_j) = math.mulEcQuartet(
+                uint256(points[i].x),
+                uint256(points[i].y),
+                uint256(points[i].k),
+                uint256(points[i+1].x),
+                uint256(points[i+1].y),
+                uint256(points[i+1].k),
+                uint256(points[i+2].x),
+                uint256(points[i+2].y),
+                uint256(points[i+2].k),
+                uint256(points[i+3].x),
+                uint256(points[i+3].y),
+                uint256(points[i+3].k)
+            );
+
+            assertEq(x_p, x_j);
+            assertEq(y_p, y_j);
+            assertEq(math.isOnCurve(x_p, y_p), true);
+
+            unchecked { ++i; }
+        }
+    }
+
+    function test_pair() public view {
+        uint256 len = points.length - 1;
+        for (uint256 i; i < len;) {
+            (uint256 x_p, uint256 y_p, uint256 z_p) = math.mulProjectivePair(
+                uint256(points[i].x),
+                uint256(points[i].y),
+                uint256(points[i].k),
+                uint256(points[i+1].x),
+                uint256(points[i+1].y),
+                uint256(points[i+1].k)
+            );
+            (x_p, y_p) = math.toAffineProjective(x_p, y_p, z_p);
+
+            (uint256 x_j, uint256 y_j) = math.mulEcPair(
+                uint256(points[i].x),
+                uint256(points[i].y),
+                uint256(points[i].k),
+                uint256(points[i+1].x),
+                uint256(points[i+1].y),
+                uint256(points[i+1].k)
+            );
+
+            assertEq(x_p, x_j);
+            assertEq(y_p, y_j);
+            assertEq(math.isOnCurve(x_p, y_p), true);
+
+            unchecked { ++i; }
+        }
+    }
+
     function test_triplet() public view {
         uint256 len = points.length - 2;
         for (uint256 i; i < len;) {
@@ -49,10 +121,10 @@ contract MathTesterTest is Test {
                 uint256(points[i+2].y),
                 uint256(points[i+2].k)
             );
-            // (x_j, y_j) = math.toAffineJacobian(x_j, y_j, z_j);
 
             assertEq(x_p, x_j);
             assertEq(y_p, y_j);
+            assertEq(math.isOnCurve(x_p, y_p), true);
 
             unchecked { ++i; }
         }
@@ -79,6 +151,7 @@ contract MathTesterTest is Test {
 
             assertEq(x_p, x_j);
             assertEq(y_p, y_j);
+            assertEq(math.isOnCurve(x_p, y_p), true);
 
             unchecked { ++i; }
         }
@@ -105,6 +178,7 @@ contract MathTesterTest is Test {
 
             assertEq(x_p, x_j);
             assertEq(y_p, y_j);
+            assertEq(math.isOnCurve(x_p, y_p), true);
 
             unchecked { ++i; }
         }
@@ -127,6 +201,7 @@ contract MathTesterTest is Test {
 
             assertEq(x_p, x_j);
             assertEq(y_p, y_j);
+            assertEq(math.isOnCurve(x_p, y_p), true);
 
             unchecked { ++i; }
         }