From 690c465838f9492a465e5f9a2692eea8dce33db1 Mon Sep 17 00:00:00 2001 From: Uncle Fatso Date: Sun, 12 Oct 2025 16:25:36 +0300 Subject: [PATCH] initial commit; double, addition, division are already optimized Signed-off-by: Uncle Fatso --- README.md | 62 +++- foundry.toml | 9 + raw_vectors.json | 227 +++++++++++++++ script/Counter.s.sol | 19 -- src/Counter.sol | 14 - src/MathTester.sol | 47 +++ src/libraries/ECMath.sol | 446 +++++++++++++++++++++++++++++ src/libraries/ECMathProjective.sol | 127 ++++++++ test/Counter.t.sol | 24 -- test/MathTester.t.sol | 72 +++++ 10 files changed, 986 insertions(+), 61 deletions(-) create mode 100644 raw_vectors.json delete mode 100644 script/Counter.s.sol delete mode 100644 src/Counter.sol create mode 100644 src/MathTester.sol create mode 100644 src/libraries/ECMath.sol create mode 100644 src/libraries/ECMathProjective.sol delete mode 100644 test/Counter.t.sol create mode 100644 test/MathTester.t.sol diff --git a/README.md b/README.md index 4db6a44..0330b34 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,63 @@ # EXODUS - EXchange Of Digital Uniformed Signatures -## Description +## Overview -For more information [visit wiki](https://git.ghostchain.io/ghostchain/ghost-exodus-draft/wiki/Description). +This module optimizes a hot spot in the verification formula used for missing‑signer recovery. The target expression is: -## Implementation +`k*P + l*Q + d*M` -An implementation of the EXODUS verifier will be in this repository. +which requires three scalar multiplications and two point additions. Naive elliptic‑curve routines make this very gas‑expensive; this work reduces cost while maintaining correctness. + +For full background and protocol details see the [project wiki](https://git.ghostchain.io/ghostchain/ghost-exodus-draft/wiki/Description). + +## Goals + +* Reduce gas for the targeted combination of scalar multiplications and additions. +* Keep implementation compact and auditable for on‑chain use. +* Maintain correctness and safety for cryptographic operations. + +## Design choices + +* Use Projective coordinates (not Jacobian) to cut down on the number of `mulmod`/`addmod` operations where possible while retaining simple formulas for point addition and doubling. +* Perform the final conversion to affine coordinates with an optimized `Extended Euclidean Algorithm` implemented in inline assembly to reduce gas compared with high‑level inversion routines. +* Benchmark against the Jacobian implementation from the [witnet elliptic‑curve‑solidity project](https://github.com/witnet/elliptic-curve-solidity) as a reference. + +## Rationale + +* Projective coordinates: fewer modular multiplications in the common path, making point operations cheaper on average. +* Assembly `Extended Euclidean Algorithm` for inversion: this algorithm in optimized inline assembly typically yields lower gas for single inversions compared with repeated `mulmod` exponentiation or other higher‑level approaches. +* Comparing to a well‑maintained Jacobian implementation gives a meaningful baseline for gas and correctness. + +## Gas Usage + +`Jacobian` is the original implementation used as a reference implementation, while `Projective` is optimized one. + +``` +╭----------------------------------------+-----------------+-------+--------+-------+---------╮ +| src/MathTester.sol:MathTester Contract | | | | | | ++=============================================================================================+ +| Deployment Cost | Deployment Size | | | | | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| 1065323 | 4715 | | | | | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| | | | | | | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| Function Name | Min | Avg | Median | Max | # Calls | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| addJacobian | 1768 | 1768 | 1768 | 1768 | 44 | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| addProjective | 992 | 992 | 992 | 992 | 44 | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| doubleJacobian | 777 | 777 | 777 | 777 | 45 | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| doubleProjective | 532 | 532 | 532 | 532 | 45 | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| toAffineJacobian | 30564 | 36206 | 36300 | 40841 | 89 | +|----------------------------------------+-----------------+-------+--------+-------+---------| +| toAffineProjective | 11838 | 13792 | 13796 | 15155 | 89 | +╰----------------------------------------+-----------------+-------+--------+-------+---------╯ +``` + +## Contributing + +All contributions are welcome — whether it's code, documentation, tests, performance benchmarks, or review. Please submit commits, issues, or pull requests; any help to improve correctness, security, or gas efficiency is greatly appreciated. diff --git a/foundry.toml b/foundry.toml index 25b918f..bb1bc14 100644 --- a/foundry.toml +++ b/foundry.toml @@ -1,6 +1,15 @@ [profile.default] +solc_version = "0.8.30" src = "src" out = "out" libs = ["lib"] +via_ir = true +optimizer = true +optimizer_runs = 4294967295 + +fs_permissions = [ + { access = "read", path = "./raw_vectors.json" }, +] + # See more config options https://github.com/foundry-rs/foundry/blob/master/crates/config/README.md#all-options diff --git a/raw_vectors.json b/raw_vectors.json new file mode 100644 index 0000000..09ce9bb --- /dev/null +++ b/raw_vectors.json @@ -0,0 +1,227 @@ +[ + { + "k": 1, + "x": "0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798", + "y": "0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8" + }, + { + "k": 2, + "x": "0xc6047f9441ed7d6d3045406e95c07cd85c778e4b8cef3ca7abac09b95c709ee5", + "y": "0x1ae168fea63dc339a3c58419466ceaeef7f632653266d0e1236431a950cfe52a" + }, + { + "k": 3, + "x": "0xf9308a019258c31049344f85f89d5229b531c845836f99b08601f113bce036f9", + "y": "0x388f7b0f632de8140fe337e62a37f3566500a99934c2231b6cb9fd7584b8e672" + }, + { + "k": 4, + "x": "0xe493dbf1c10d80f3581e4904930b1404cc6c13900ee0758474fa94abe8c4cd13", + "y": "0x51ed993ea0d455b75642e2098ea51448d967ae33bfbdfe40cfe97bdc47739922" + }, + { + "k": 5, + "x": "0x2f8bde4d1a07209355b4a7250a5c5128e88b84bddc619ab7cba8d569b240efe4", + "y": "0xd8ac222636e5e3d6d4dba9dda6c9c426f788271bab0d6840dca87d3aa6ac62d6" + }, + { + "k": 6, + "x": "0xfff97bd5755eeea420453a14355235d382f6472f8568a18b2f057a1460297556", + "y": "0xae12777aacfbb620f3be96017f45c560de80f0f6518fe4a03c870c36b075f297" + }, + { + "k": 7, + "x": "0x5cbdf0646e5db4eaa398f365f2ea7a0e3d419b7e0330e39ce92bddedcac4f9bc", + "y": "0x6aebca40ba255960a3178d6d861a54dba813d0b813fde7b5a5082628087264da" + }, + { + "k": 8, + "x": "0x2f01e5e15cca351daff3843fb70f3c2f0a1bdd05e5af888a67784ef3e10a2a01", + "y": "0x5c4da8a741539949293d082a132d13b4c2e213d6ba5b7617b5da2cb76cbde904" + }, + { + "k": 9, + "x": "0xacd484e2f0c7f65309ad178a9f559abde09796974c57e714c35f110dfc27ccbe", + "y": "0xcc338921b0a7d9fd64380971763b61e9add888a4375f8e0f05cc262ac64f9c37" + }, + { + "k": 10, + "x": "0xa0434d9e47f3c86235477c7b1ae6ae5d3442d49b1943c2b752a68e2a47e247c7", + "y": "0x893aba425419bc27a3b6c7e693a24c696f794c2ed877a1593cbee53b037368d7" + }, + { + "k": 11, + "x": "0x774ae7f858a9411e5ef4246b70c65aac5649980be5c17891bbec17895da008cb", + "y": "0xd984a032eb6b5e190243dd56d7b7b365372db1e2dff9d6a8301d74c9c953c61b" + }, + { + "k": 12, + "x": "0xd01115d548e7561b15c38f004d734633687cf4419620095bc5b0f47070afe85a", + "y": "0xa9f34ffdc815e0d7a8b64537e17bd81579238c5dd9a86d526b051b13f4062327" + }, + { + "k": 13, + "x": "0xf28773c2d975288bc7d1d205c3748651b075fbc6610e58cddeeddf8f19405aa8", + "y": "0x0ab0902e8d880a89758212eb65cdaf473a1a06da521fa91f29b5cb52db03ed81" + }, + { + "k": 14, + "x": "0x499fdf9e895e719cfd64e67f07d38e3226aa7b63678949e6e49b241a60e823e4", + "y": "0xcac2f6c4b54e855190f044e4a7b3d464464279c27a3f95bcc65f40d403a13f5b" + }, + { + "k": 15, + "x": "0xd7924d4f7d43ea965a465ae3095ff41131e5946f3c85f79e44adbcf8e27e080e", + "y": "0x581e2872a86c72a683842ec228cc6defea40af2bd896d3a5c504dc9ff6a26b58" + }, + { + "k": 16, + "x": "0xe60fce93b59e9ec53011aabc21c23e97b2a31369b87a5ae9c44ee89e2a6dec0a", + "y": "0xf7e3507399e595929db99f34f57937101296891e44d23f0be1f32cce69616821" + }, + { + "k": 17, + "x": "0xdefdea4cdb677750a420fee807eacf21eb9898ae79b9768766e4faa04a2d4a34", + "y": "0x4211ab0694635168e997b0ead2a93daeced1f4a04a95c0f6cfb199f69e56eb77" + }, + { + "k": 18, + "x": "0x5601570cb47f238d2b0286db4a990fa0f3ba28d1a319f5e7cf55c2a2444da7cc", + "y": "0xc136c1dc0cbeb930e9e298043589351d81d8e0bc736ae2a1f5192e5e8b061d58" + }, + { + "k": 19, + "x": "0x2b4ea0a797a443d293ef5cff444f4979f06acfebd7e86d277475656138385b6c", + "y": "0x85e89bc037945d93b343083b5a1c86131a01f60c50269763b570c854e5c09b7a" + }, + { + "k": 20, + "x": "0x4ce119c96e2fa357200b559b2f7dd5a5f02d5290aff74b03f3e471b273211c97", + "y": "0x12ba26dcb10ec1625da61fa10a844c676162948271d96967450288ee9233dc3a" + }, + { + "k": 112233445566778899, + "x": "0xa90cc3d3f3e146daadfc74ca1372207cb4b725ae708cef713a98edd73d99ef29", + "y": "0x5a79d6b289610c68bc3b47f3d72f9788a26a06868b4d8e433e1e2ad76fb7dc76" + }, + { + "k": 112233445566778899112233445566778899, + "x": "0xe5a2636bcfd412ebf36ec45b19bfb68a1bc5f8632e678132b885f7df99c5e9b3", + "y": "0x736c1ce161ae27b405cafd2a7520370153c2c861ac51d6c1d5985d9606b45f39" + }, + { + "k": 28948022309329048855892746252171976963209391069768726095651290785379540373584, + "x": "0xa6b594b38fb3e77c6edf78161fade2041f4e09fd8497db776e546c41567feb3c", + "y": "0x71444009192228730cd8237a490feba2afe3d27d7cc1136bc97e439d13330d55" + }, + { + "k": 57896044618658097711785492504343953926418782139537452191302581570759080747168, + "x": "0x00000000000000000000003b78ce563f89a0ed9414f5aa28ad0d96d6795f9c63", + "y": "0x3f3979bf72ae8202983dc989aec7f2ff2ed91bdd69ce02fc0700ca100e59ddf3" + }, + { + "k": 86844066927987146567678238756515930889628173209306178286953872356138621120752, + "x": "0xe24ce4beee294aa6350faa67512b99d388693ae4e7f53d19882a6ea169fc1ce1", + "y": "0x8b71e83545fc2b5872589f99d948c03108d36797c4de363ebd3ff6a9e1a95b10" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494317, + "x": "0x4ce119c96e2fa357200b559b2f7dd5a5f02d5290aff74b03f3e471b273211c97", + "y": "0xed45d9234ef13e9da259e05ef57bb3989e9d6b7d8e269698bafd77106dcc1ff5" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494318, + "x": "0x2b4ea0a797a443d293ef5cff444f4979f06acfebd7e86d277475656138385b6c", + "y": "0x7a17643fc86ba26c4cbcf7c4a5e379ece5fe09f3afd9689c4a8f37aa1a3f60b5" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494319, + "x": "0x5601570cb47f238d2b0286db4a990fa0f3ba28d1a319f5e7cf55c2a2444da7cc", + "y": "0x3ec93e23f34146cf161d67fbca76cae27e271f438c951d5e0ae6d1a074f9ded7" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494320, + "x": "0xdefdea4cdb677750a420fee807eacf21eb9898ae79b9768766e4faa04a2d4a34", + "y": "0xbdee54f96b9cae9716684f152d56c251312e0b5fb56a3f09304e660861a910b8" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494321, + "x": "0xe60fce93b59e9ec53011aabc21c23e97b2a31369b87a5ae9c44ee89e2a6dec0a", + "y": "0x081caf8c661a6a6d624660cb0a86c8efed6976e1bb2dc0f41e0cd330969e940e" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494322, + "x": "0xd7924d4f7d43ea965a465ae3095ff41131e5946f3c85f79e44adbcf8e27e080e", + "y": "0xa7e1d78d57938d597c7bd13dd733921015bf50d427692c5a3afb235f095d90d7" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494323, + "x": "0x499fdf9e895e719cfd64e67f07d38e3226aa7b63678949e6e49b241a60e823e4", + "y": "0x353d093b4ab17aae6f0fbb1b584c2b9bb9bd863d85c06a4339a0bf2afc5ebcd4" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494324, + "x": "0xf28773c2d975288bc7d1d205c3748651b075fbc6610e58cddeeddf8f19405aa8", + "y": "0xf54f6fd17277f5768a7ded149a3250b8c5e5f925ade056e0d64a34ac24fc0eae" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494325, + "x": "0xd01115d548e7561b15c38f004d734633687cf4419620095bc5b0f47070afe85a", + "y": "0x560cb00237ea1f285749bac81e8427ea86dc73a2265792ad94fae4eb0bf9d908" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494326, + "x": "0x774ae7f858a9411e5ef4246b70c65aac5649980be5c17891bbec17895da008cb", + "y": "0x267b5fcd1494a1e6fdbc22a928484c9ac8d24e1d20062957cfe28b3536ac3614" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494327, + "x": "0xa0434d9e47f3c86235477c7b1ae6ae5d3442d49b1943c2b752a68e2a47e247c7", + "y": "0x76c545bdabe643d85c4938196c5db3969086b3d127885ea6c3411ac3fc8c9358" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494328, + "x": "0xacd484e2f0c7f65309ad178a9f559abde09796974c57e714c35f110dfc27ccbe", + "y": "0x33cc76de4f5826029bc7f68e89c49e165227775bc8a071f0fa33d9d439b05ff8" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494329, + "x": "0x2f01e5e15cca351daff3843fb70f3c2f0a1bdd05e5af888a67784ef3e10a2a01", + "y": "0xa3b25758beac66b6d6c2f7d5ecd2ec4b3d1dec2945a489e84a25d3479342132b" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494330, + "x": "0x5cbdf0646e5db4eaa398f365f2ea7a0e3d419b7e0330e39ce92bddedcac4f9bc", + "y": "0x951435bf45daa69f5ce8729279e5ab2457ec2f47ec02184a5af7d9d6f78d9755" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494331, + "x": "0xfff97bd5755eeea420453a14355235d382f6472f8568a18b2f057a1460297556", + "y": "0x51ed8885530449df0c4169fe80ba3a9f217f0f09ae701b5fc378f3c84f8a0998" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494332, + "x": "0x2f8bde4d1a07209355b4a7250a5c5128e88b84bddc619ab7cba8d569b240efe4", + "y": "0x2753ddd9c91a1c292b24562259363bd90877d8e454f297bf235782c459539959" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494333, + "x": "0xe493dbf1c10d80f3581e4904930b1404cc6c13900ee0758474fa94abe8c4cd13", + "y": "0xae1266c15f2baa48a9bd1df6715aebb7269851cc404201bf30168422b88c630d" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494334, + "x": "0xf9308a019258c31049344f85f89d5229b531c845836f99b08601f113bce036f9", + "y": "0xc77084f09cd217ebf01cc819d5c80ca99aff5666cb3ddce4934602897b4715bd" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494335, + "x": "0xc6047f9441ed7d6d3045406e95c07cd85c778e4b8cef3ca7abac09b95c709ee5", + "y": "0xe51e970159c23cc65c3a7be6b99315110809cd9acd992f1edc9bce55af301705" + }, + { + "k": 115792089237316195423570985008687907852837564279074904382605163141518161494336, + "x": "0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798", + "y": "0xb7c52588d95c3b9aa25b0403f1eef75702e84bb7597aabe663b82f6f04ef2777" + } + ] diff --git a/script/Counter.s.sol b/script/Counter.s.sol deleted file mode 100644 index cdc1fe9..0000000 --- a/script/Counter.s.sol +++ /dev/null @@ -1,19 +0,0 @@ -// SPDX-License-Identifier: UNLICENSED -pragma solidity ^0.8.13; - -import {Script, console} from "forge-std/Script.sol"; -import {Counter} from "../src/Counter.sol"; - -contract CounterScript is Script { - Counter public counter; - - function setUp() public {} - - function run() public { - vm.startBroadcast(); - - counter = new Counter(); - - vm.stopBroadcast(); - } -} diff --git a/src/Counter.sol b/src/Counter.sol deleted file mode 100644 index aded799..0000000 --- a/src/Counter.sol +++ /dev/null @@ -1,14 +0,0 @@ -// SPDX-License-Identifier: UNLICENSED -pragma solidity ^0.8.13; - -contract Counter { - uint256 public number; - - function setNumber(uint256 newNumber) public { - number = newNumber; - } - - function increment() public { - number++; - } -} diff --git a/src/MathTester.sol b/src/MathTester.sol new file mode 100644 index 0000000..99e93cd --- /dev/null +++ b/src/MathTester.sol @@ -0,0 +1,47 @@ +// SPDX-License-Identifier: UNLICENSED +pragma solidity ^0.8.0; + +import {EllipticCurve} from "./libraries/ECMath.sol"; +import {EllipticCurveProjective} from "./libraries/ECMathProjective.sol"; + +contract MathTester { + // Constants are taken from https://en.bitcoin.it/wiki/Secp256k1 + uint256 public constant A = 0; + uint256 public constant B = 7; + uint256 public constant P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F; + uint256 public constant N = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141; + + function addJacobian( + uint256 x1, + uint256 y1, + uint256 x2, + uint256 y2 + ) public pure returns (uint256, uint256, uint256) { + return EllipticCurve.jacAdd(x1, y1, 1, x2, y2, 1, P); + } + + function addProjective( + uint256 x1, + uint256 y1, + uint256 x2, + uint256 y2 + ) public pure returns (uint256, uint256, uint256) { + return EllipticCurveProjective.projectiveAdd(x1, y1, 1, x2, y2, 1); + } + + function doubleJacobian(uint256 x1, uint256 y1) public pure returns (uint256, uint256, uint256) { + return EllipticCurve.jacDouble(x1, y1, 1, A, P); + } + + function doubleProjective(uint256 x1, uint256 y1) public pure returns (uint256, uint256, uint256) { + return EllipticCurveProjective.projectiveDouble(x1, y1, 1); + } + + function toAffineJacobian(uint256 x, uint256 y, uint256 z) public pure returns (uint256, uint256) { + return EllipticCurve.toAffine(x, y, z, P); + } + + function toAffineProjective(uint256 x, uint256 y, uint256 z) public pure returns (uint256, uint256) { + return EllipticCurveProjective.toAffine(x, y, z); + } +} diff --git a/src/libraries/ECMath.sol b/src/libraries/ECMath.sol new file mode 100644 index 0000000..8e7b41c --- /dev/null +++ b/src/libraries/ECMath.sol @@ -0,0 +1,446 @@ +// SPDX-License-Identifier: MIT + +pragma solidity ^0.8.0; + +/** + ** @title Elliptic Curve Library + ** @dev Library providing arithmetic operations over elliptic curves. + ** This library does not check whether the inserted points belong to the curve + ** `isOnCurve` function should be used by the library user to check the aforementioned statement. + ** @author Witnet Foundation + */ +library EllipticCurve { + // Pre-computed constant for 2 ** 255 + uint256 private constant U255_MAX_PLUS_1 = + 57896044618658097711785492504343953926634992332820282019728792003956564819968; + + /// @dev Modular euclidean inverse of a number (mod p). + /// @param _x The number + /// @param _pp The modulus + /// @return q such that x*q = 1 (mod _pp) + function invMod(uint256 _x, uint256 _pp) internal pure returns (uint256) { + require(_x != 0 && _x != _pp && _pp != 0, "Invalid number"); + uint256 q = 0; + uint256 newT = 1; + uint256 r = _pp; + uint256 t; + while (_x != 0) { + t = r / _x; + (q, newT) = (newT, addmod(q, (_pp - mulmod(t, newT, _pp)), _pp)); + (r, _x) = (_x, r - t * _x); + } + + return q; + } + + /// @dev Modular exponentiation, b^e % _pp. + /// Source: https://github.com/androlo/standard-contracts/blob/master/contracts/src/crypto/ECCMath.sol + /// @param _base base + /// @param _exp exponent + /// @param _pp modulus + /// @return r such that r = b**e (mod _pp) + function expMod( + uint256 _base, + uint256 _exp, + uint256 _pp + ) + internal pure + returns (uint256) + { + require(_pp != 0, "EllipticCurve: modulus is zero"); + + if (_base == 0) return 0; + if (_exp == 0) return 1; + + uint256 r = 1; + uint256 bit = U255_MAX_PLUS_1; + assembly { + for { + + } gt(bit, 0) { + + } { + r := mulmod( + mulmod(r, r, _pp), + exp(_base, iszero(iszero(and(_exp, bit)))), + _pp + ) + r := mulmod( + mulmod(r, r, _pp), + exp(_base, iszero(iszero(and(_exp, div(bit, 2))))), + _pp + ) + r := mulmod( + mulmod(r, r, _pp), + exp(_base, iszero(iszero(and(_exp, div(bit, 4))))), + _pp + ) + r := mulmod( + mulmod(r, r, _pp), + exp(_base, iszero(iszero(and(_exp, div(bit, 8))))), + _pp + ) + bit := div(bit, 16) + } + } + + return r; + } + + /// @dev Converts a point (x, y, z) expressed in Jacobian coordinates to affine coordinates (x', y', 1). + /// @param _x coordinate x + /// @param _y coordinate y + /// @param _z coordinate z + /// @param _pp the modulus + /// @return (x', y') affine coordinates + function toAffine( + uint256 _x, + uint256 _y, + uint256 _z, + uint256 _pp + ) + internal pure + returns (uint256, uint256) + { + uint256 zInv = invMod(_z, _pp); + uint256 zInv2 = mulmod(zInv, zInv, _pp); + uint256 x2 = mulmod(_x, zInv2, _pp); + uint256 y2 = mulmod(_y, mulmod(zInv, zInv2, _pp), _pp); + + return (x2, y2); + } + + /// @dev Derives the y coordinate from a compressed-format point x [[SEC-1]](https://www.secg.org/SEC1-Ver-1.0.pdf). + /// @param _prefix parity byte (0x02 even, 0x03 odd) + /// @param _x coordinate x + /// @param _aa constant of curve + /// @param _bb constant of curve + /// @param _pp the modulus + /// @return y coordinate y + function deriveY( + uint8 _prefix, + uint256 _x, + uint256 _aa, + uint256 _bb, + uint256 _pp + ) + internal pure + returns (uint256) + { + require( + _prefix == 0x02 || _prefix == 0x03, + "EllipticCurve:innvalid compressed EC point prefix" + ); + + // x^3 + ax + b + uint256 y2 = addmod( + mulmod(_x, mulmod(_x, _x, _pp), _pp), + addmod(mulmod(_x, _aa, _pp), _bb, _pp), + _pp + ); + y2 = expMod(y2, (_pp + 1) / 4, _pp); + // uint256 cmp = yBit ^ y_ & 1; + uint256 y = (y2 + _prefix) % 2 == 0 ? y2 : _pp - y2; + + return y; + } + + /// @dev Check whether point (x,y) is on curve defined by a, b, and _pp. + /// @param _x coordinate x of P1 + /// @param _y coordinate y of P1 + /// @param _aa constant of curve + /// @param _bb constant of curve + /// @param _pp the modulus + /// @return true if x,y in the curve, false else + function isOnCurve( + uint _x, + uint _y, + uint _aa, + uint _bb, + uint _pp + ) + internal pure + returns (bool) + { + if (0 == _x || _x >= _pp || 0 == _y || _y >= _pp) { + return false; + } + // y^2 + uint lhs = mulmod(_y, _y, _pp); + // x^3 + uint rhs = mulmod(mulmod(_x, _x, _pp), _x, _pp); + if (_aa != 0) { + // x^3 + a*x + rhs = addmod(rhs, mulmod(_x, _aa, _pp), _pp); + } + if (_bb != 0) { + // x^3 + a*x + b + rhs = addmod(rhs, _bb, _pp); + } + + return lhs == rhs; + } + + /// @dev Calculate inverse (x, -y) of point (x, y). + /// @param _x coordinate x of P1 + /// @param _y coordinate y of P1 + /// @param _pp the modulus + /// @return (x, -y) + function ecInv( + uint256 _x, + uint256 _y, + uint256 _pp + ) + internal pure + returns (uint256, uint256) + { + return (_x, (_pp - _y) % _pp); + } + + /// @dev Add two points (x1, y1) and (x2, y2) in affine coordinates. + /// @param _x1 coordinate x of P1 + /// @param _y1 coordinate y of P1 + /// @param _x2 coordinate x of P2 + /// @param _y2 coordinate y of P2 + /// @param _aa constant of the curve + /// @param _pp the modulus + /// @return (qx, qy) = P1+P2 in affine coordinates + function ecAdd( + uint256 _x1, + uint256 _y1, + uint256 _x2, + uint256 _y2, + uint256 _aa, + uint256 _pp + ) + internal pure + returns (uint256, uint256) + { + uint x = 0; + uint y = 0; + uint z = 0; + + // Double if x1==x2 else add + if (_x1 == _x2) { + // y1 = -y2 mod p + if (addmod(_y1, _y2, _pp) == 0) { + return (0, 0); + } else { + // P1 = P2 + (x, y, z) = jacDouble(_x1, _y1, 1, _aa, _pp); + } + } else { + (x, y, z) = jacAdd(_x1, _y1, 1, _x2, _y2, 1, _pp); + } + // Get back to affine + return toAffine(x, y, z, _pp); + } + + /// @dev Substract two points (x1, y1) and (x2, y2) in affine coordinates. + /// @param _x1 coordinate x of P1 + /// @param _y1 coordinate y of P1 + /// @param _x2 coordinate x of P2 + /// @param _y2 coordinate y of P2 + /// @param _aa constant of the curve + /// @param _pp the modulus + /// @return (qx, qy) = P1-P2 in affine coordinates + function ecSub( + uint256 _x1, + uint256 _y1, + uint256 _x2, + uint256 _y2, + uint256 _aa, + uint256 _pp + ) + internal pure + returns (uint256, uint256) + { + // invert square + (uint256 x, uint256 y) = ecInv(_x2, _y2, _pp); + // P1-square + return ecAdd(_x1, _y1, x, y, _aa, _pp); + } + + /// @dev Multiply point (x1, y1, z1) times d in affine coordinates. + /// @param _k scalar to multiply + /// @param _x coordinate x of P1 + /// @param _y coordinate y of P1 + /// @param _aa constant of the curve + /// @param _pp the modulus + /// @return (qx, qy) = d*P in affine coordinates + function ecMul( + uint256 _k, + uint256 _x, + uint256 _y, + uint256 _aa, + uint256 _pp + ) + internal pure + returns (uint256, uint256) + { + // Jacobian multiplication + (uint256 x1, uint256 y1, uint256 z1) = jacMul(_k, _x, _y, 1, _aa, _pp); + // Get back to affine + return toAffine(x1, y1, z1, _pp); + } + + /// @dev Adds two points (x1, y1, z1) and (x2 y2, z2). + /// @param _x1 coordinate x of P1 + /// @param _y1 coordinate y of P1 + /// @param _z1 coordinate z of P1 + /// @param _x2 coordinate x of square + /// @param _y2 coordinate y of square + /// @param _z2 coordinate z of square + /// @param _pp the modulus + /// @return (qx, qy, qz) P1+square in Jacobian + function jacAdd( + uint256 _x1, + uint256 _y1, + uint256 _z1, + uint256 _x2, + uint256 _y2, + uint256 _z2, + uint256 _pp + ) + internal pure + returns (uint256, uint256, uint256) + { + if (_x1 == 0 && _y1 == 0) return (_x2, _y2, _z2); + if (_x2 == 0 && _y2 == 0) return (_x1, _y1, _z1); + + // We follow the equations described in https://pdfs.semanticscholar.org/5c64/29952e08025a9649c2b0ba32518e9a7fb5c2.pdf Section 5 + uint[4] memory zs; // z1^2, z1^3, z2^2, z2^3 + zs[0] = mulmod(_z1, _z1, _pp); + zs[1] = mulmod(_z1, zs[0], _pp); + zs[2] = mulmod(_z2, _z2, _pp); + zs[3] = mulmod(_z2, zs[2], _pp); + + // u1, s1, u2, s2 + zs = [ + mulmod(_x1, zs[2], _pp), + mulmod(_y1, zs[3], _pp), + mulmod(_x2, zs[0], _pp), + mulmod(_y2, zs[1], _pp) + ]; + + // In case of zs[0] == zs[2] && zs[1] == zs[3], double function should be used + require( + zs[0] != zs[2] || zs[1] != zs[3], + "Use jacDouble function instead" + ); + + uint[4] memory hr; + //h + hr[0] = addmod(zs[2], _pp - zs[0], _pp); + //r + hr[1] = addmod(zs[3], _pp - zs[1], _pp); + //h^2 + hr[2] = mulmod(hr[0], hr[0], _pp); + // h^3 + hr[3] = mulmod(hr[2], hr[0], _pp); + // qx = -h^3 -2u1h^2+r^2 + uint256 qx = addmod(mulmod(hr[1], hr[1], _pp), _pp - hr[3], _pp); + qx = addmod(qx, _pp - mulmod(2, mulmod(zs[0], hr[2], _pp), _pp), _pp); + // qy = -s1*z1*h^3+r(u1*h^2 -x^3) + uint256 qy = mulmod( + hr[1], + addmod(mulmod(zs[0], hr[2], _pp), _pp - qx, _pp), + _pp + ); + qy = addmod(qy, _pp - mulmod(zs[1], hr[3], _pp), _pp); + // qz = h*z1*z2 + uint256 qz = mulmod(hr[0], mulmod(_z1, _z2, _pp), _pp); + return (qx, qy, qz); + } + + /// @dev Doubles a points (x, y, z). + /// @param _x coordinate x of P1 + /// @param _y coordinate y of P1 + /// @param _z coordinate z of P1 + /// @param _aa the a scalar in the curve equation + /// @param _pp the modulus + /// @return (qx, qy, qz) 2P in Jacobian + function jacDouble( + uint256 _x, + uint256 _y, + uint256 _z, + uint256 _aa, + uint256 _pp + ) + internal pure + returns (uint256, uint256, uint256) + { + if (_z == 0) return (_x, _y, _z); + + // We follow the equations described in https://pdfs.semanticscholar.org/5c64/29952e08025a9649c2b0ba32518e9a7fb5c2.pdf Section 5 + // Note: there is a bug in the paper regarding the m parameter, M=3*(x1^2)+a*(z1^4) + // x, y, z at this point represent the squares of _x, _y, _z + uint256 x = mulmod(_x, _x, _pp); //x1^2 + uint256 y = mulmod(_y, _y, _pp); //y1^2 + uint256 z = mulmod(_z, _z, _pp); //z1^2 + + // s + uint s = mulmod(4, mulmod(_x, y, _pp), _pp); + // m + uint m = addmod( + mulmod(3, x, _pp), + mulmod(_aa, mulmod(z, z, _pp), _pp), + _pp + ); + + // x, y, z at this point will be reassigned and rather represent qx, qy, qz from the paper + // This allows to reduce the gas cost and stack footprint of the algorithm + // qx + x = addmod(mulmod(m, m, _pp), _pp - addmod(s, s, _pp), _pp); + // qy = -8*y1^4 + M(S-T) + y = addmod( + mulmod(m, addmod(s, _pp - x, _pp), _pp), + _pp - mulmod(8, mulmod(y, y, _pp), _pp), + _pp + ); + // qz = 2*y1*z1 + z = mulmod(2, mulmod(_y, _z, _pp), _pp); + + return (x, y, z); + } + + /// @dev Multiply point (x, y, z) times d. + /// @param _d scalar to multiply + /// @param _x coordinate x of P1 + /// @param _y coordinate y of P1 + /// @param _z coordinate z of P1 + /// @param _aa constant of curve + /// @param _pp the modulus + /// @return (qx, qy, qz) d*P1 in Jacobian + function jacMul( + uint256 _d, + uint256 _x, + uint256 _y, + uint256 _z, + uint256 _aa, + uint256 _pp + ) + internal pure + returns (uint256, uint256, uint256) + { + // Early return in case that `_d == 0` + if (_d == 0) { + return (_x, _y, _z); + } + + uint256 remaining = _d; + uint256 qx = 0; + uint256 qy = 0; + uint256 qz = 1; + + // Double and add algorithm + while (remaining != 0) { + if ((remaining & 1) != 0) { + (qx, qy, qz) = jacAdd(qx, qy, qz, _x, _y, _z, _pp); + } + remaining = remaining / 2; + (_x, _y, _z) = jacDouble(_x, _y, _z, _aa, _pp); + } + return (qx, qy, qz); + } +} diff --git a/src/libraries/ECMathProjective.sol b/src/libraries/ECMathProjective.sol new file mode 100644 index 0000000..81ed896 --- /dev/null +++ b/src/libraries/ECMathProjective.sol @@ -0,0 +1,127 @@ +// SPDX-License-Identifier: MIT + +pragma solidity ^0.8.0; + +// Vectors for secp256k1 are difficult to find. These are the vectors from: +// https://web.archive.org/web/20190724010836/https://chuckbatson.wordpress.com/2014/11/26/secp256k1-test-vectors + +library EllipticCurveProjective { + // Constants are taken from https://en.bitcoin.it/wiki/Secp256k1 + uint256 private constant A = 0; + uint256 private constant B = 7; + uint256 private constant P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F; + uint256 private constant N = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141; + + uint256 private constant B3 = 21; + + function projectiveAdd( + uint256 X1, + uint256 Y1, + uint256 Z1, + uint256 X2, + uint256 Y2, + uint256 Z2 + ) internal pure returns (uint256 X3, uint256 Y3, uint256 Z3) { + // We implement the complete addition formula from Renes-Costello-Batina 2015 + // (https://eprint.iacr.org/2015/1060 Algorithm 7). + + // X3 = (X1Y2 + X2Y1)(Y1Y2 − 3bZ1Z2) − 3b(Y1Z2 + Y2Z1)(X1Z2 + X2Z1), + // 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) + uint256 t3 = addmod(X1, Y1, P); // 4. t3 ← X1 + Y1 => (X1 + Y1) + uint256 t4 = addmod(X2, Y2, P); // 5. t4 ← X2 + Y2 => (X2 + Y2) + t3 = mulmod(t3, t4, P); // 6. t3 ← t3 · t4 => ((X1 + Y1) · (X2 + Y2)) + t4 = addmod(t0, t1, P); // 7. t4 ← t0 + t1 => (X1·X2 + Y1·Y2) + t3 = addmod(t3, P - t4, P); // 8. t3 ← t3 - t4 => ((X1 + Y1)·(X2 + Y2) - X1·X2 - Y1·Y2) + t4 = addmod(Y1, Z1, P); // 9. t4 ← Y1 + Z1 => (Y1 + Z1) + X3 = addmod(Y2, Z2, P); // 10. X3 ← Y2 + Z2 => (Y2 + Z2) + t4 = mulmod(t4, X3, P); // 11. t4 ← t4 · X3 => ((Y1 + Z1) · (Y2 + Z2)) + X3 = addmod(t1, t2, P); // 12. X3 ← t1 + t2 => (Y1·Y2 + Z1·Z2) + t4 = addmod(t4, P - X3, P); // 13. t4 ← t4 - X3 => ((Y1 + Z1)·(Y2 + Z2) - Y1·Y2 - Z1·Z2) + X3 = addmod(X1, Z1, P); // 14. X3 ← X1 + Z1 => (X1 + Z1) + Y3 = addmod(X2, Z2, P); // 15. Y3 ← X2 + Z2 => (X2 + Z2) + X3 = mulmod(X3, Y3, P); // 16. X3 ← X3 · Y3 => ((X1 + Z1) · (X2 + Z2)) + Y3 = addmod(t0, t2, P); // 17. Y3 ← t0 + t2 => (X1·X2 + Z1·Z2) + Y3 = addmod(X3, P - Y3, P); // 18. Y3 ← X3 - Y3 => ((X1 + Z1)·(X2 + Z2) - X1·X2 - Z1·Z2) + X3 = addmod(t0, t0, P); // 19. X3 ← t0 + t0 => (2·X1·X2) + t0 = addmod(X3, t0, P); // 20. t0 ← X3 + t0 => (3·X1·X2) + t2 = mulmod(B3, t2, P); // 21. t2 ← B3 · t2 => 3b · Z1·Z2 + Z3 = addmod(t1, t2, P); // 22. Z3 ← t1 + t2 => Y1·Y2 + 3·b·Z1·Z2 + t1 = addmod(t1, P - t2, P); // 23. t1 ← t1 - t2 => Y1·Y2 - 3·b·Z1·Z2 + Y3 = mulmod(B3, Y3, P); // 24. Y3 ← B3 · Y3 => 3b · ((X1+Z1)(X2+Z2) - X1·X2 - Z1·Z2) + X3 = mulmod(t4, Y3, P); // 25. X3 ← t4 · Y3 => 3b·((Y1+Z1)(Y2+Z2)-Y1·Y2-Z1·Z2) · ((X1+Z1)(X2+Z2)-X1·X2-Z1·Z2) + t2 = mulmod(t3, t1, P); // 26. t2 ← t3 · t1 => ((X1+Y1)(X2+Y2)-X1·X2-Y1·Y2) · (Y1·Y2 - 3·b·Z1·Z2) + X3 = addmod(t2, P - X3, P); // 27. X3 ← t2 - X3 => (X1Y2+X2Y1)(Y1·Y2-3·b·Z1·Z2) - 3b(Y1·Z2+Y2·Z1)(X1·Z2+X2·Z1) + Y3 = mulmod(Y3, t0, P); // 28. Y3 ← Y3 · t0 => 3·b·((X1+Z1)(X2+Z2) - X1·X2 - Z1·Z2) · 3·X1·X2 = 9·b·X1·X2·(X1·Z2+X2·Z1) + t1 = mulmod(t1, Z3, P); // 29. t1 ← t1 · Z3 => (Y1·Y2-3·b·Z1·Z2) · (Y1·Y2+3·b·Z1·Z2) + Y3 = addmod(t1, Y3, P); // 30. Y3 ← t1 + Y3 => (Y1Y2+3·b·Z1·Z2)(Y1·Y2-3·b·Z1·Z2) + 9b·X1·X2·(X1·Z2+X2·Z1) + t0 = mulmod(t0, t3, P); // 31. t0 ← t0 · t3 => (3·X1·X2) · ((X1+Y1)(X2+Y2)-X1·X2-Y1·Y2) + Z3 = mulmod(Z3, t4, P); // 32. Z3 ← Z3 · t4 => (Y1·Y2+3·b·Z1·Z2) · ((Y1+Z1)(Y2+Z2)-Y1·Y2-Z1·Z2) + Z3 = addmod(Z3, t0, P); // 33. Z3 ← Z3 + t0 => (Y1·Z2+Y2·Z1)·(Y1·Y2+3·b·Z1·Z2) + 3·X1·X2·(X1·Y2+X2·Y1) + } + + function projectiveDouble( + uint256 X, + uint256 Y, + uint256 Z + ) internal pure returns (uint256 X3, uint256 Y3, uint256 Z3) { + // We implement the complete addition formula from Renes-Costello-Batina 2015 + // (https://eprint.iacr.org/2015/1060 Algorithm 9). + + // X3 = 2XY (Y^2 − 9bZ^2), (ok) + // Y3 = (Y^2 − 9bZ^2)(Y^2 + 3bZ^2) + 24bY^2Z^2, + // Z3 = 8Y^3Z + + uint256 t0 = mulmod(Y, Y, P); // 1. t0 ← Y · Y => (Y²) + Z3 = mulmod(8, t0, P); // 2. Z3 ← t0 + t0 => (8Y²) + uint256 t1 = mulmod(Y, Z, P); // 5. t1 ← Y · Z => (YZ) + uint256 t2 = mulmod(Z, Z, P); // 6. t2 ← Z · Z => (Z²) + t2 = mulmod(21, t2, P); // 7. t2 ← b3 · t2 => (3bZ²) + X3 = mulmod(t2, Z3, P); // 8. X3 ← t2 · Z3 => (3bZ²8Y²) + Y3 = addmod(t0, t2, P); // 9. Y3 ← t0 + t2 => (Y² + 3bZ²) + Z3 = mulmod(t1, Z3, P); // 10. Z3 ← t1 · Z3 => (YZ · 8Y² = 8Y³Z) + t1 = addmod(t0, P - mulmod(3, t2, P), P); // 11. t1 ← t0 - (3 · t2) => (Y² - 9bZ²) + Y3 = addmod(X3, mulmod(t1, Y3, P), P); // 12. Y3 ← t1 · (t1 · Y3) => ((Y² - 9bZ²) · (Y² + 3bZ²)) + X3 = mulmod(t1, mulmod(X, Y, P), P); // 13. X3 ← t1 · (X1 · Y1) => ((Y² - 9bZ²) · XY) + X3 = addmod(X3, X3, P); // 14. X3 ← X3 + X3 => ((Y² - 9bZ²) · 2XY) + } + + function toAffine(uint256 x, uint256 y, uint256 z) internal pure returns (uint256 _x, uint256 _y) { + uint256 t0; + assembly { + // Extended Euclidean algorithm (iterative, assembly) + // Typically lower average gas than the modexp precompile for single inversions, + // but execution time (and gas) depends on input values — not constant-time. + + let t1 := 1 + let r0 := P + let r1 := z + + for {} r1 {} { + let q := div(r0, r1) + + // t0, t1 = t1, t0 - q * t1 + let t1_new := sub(t0, mul(q, t1)) + t0 := t1 + t1 := t1_new + + // r0, r1 = r1, r0 - q * r1 + let r1_new := sub(r0, mul(q, r1)) + r0 := r1 + r1 := r1_new + } + + if slt(t0, 0) { + t0 := add(t0, P) + } + } + _x = mulmod(x, t0, P); + _y = mulmod(y, t0, P); + } +} diff --git a/test/Counter.t.sol b/test/Counter.t.sol deleted file mode 100644 index 54b724f..0000000 --- a/test/Counter.t.sol +++ /dev/null @@ -1,24 +0,0 @@ -// SPDX-License-Identifier: UNLICENSED -pragma solidity ^0.8.13; - -import {Test, console} from "forge-std/Test.sol"; -import {Counter} from "../src/Counter.sol"; - -contract CounterTest is Test { - Counter public counter; - - function setUp() public { - counter = new Counter(); - counter.setNumber(0); - } - - function test_Increment() public { - counter.increment(); - assertEq(counter.number(), 1); - } - - function testFuzz_SetNumber(uint256 x) public { - counter.setNumber(x); - assertEq(counter.number(), x); - } -} diff --git a/test/MathTester.t.sol b/test/MathTester.t.sol new file mode 100644 index 0000000..e2b0120 --- /dev/null +++ b/test/MathTester.t.sol @@ -0,0 +1,72 @@ +// SPDX-License-Identifier: UNLICENSED +pragma solidity ^0.8.13; + +import {Test, stdJson, console} from "forge-std/Test.sol"; +import {MathTester} from "../src/MathTester.sol"; + +struct Point { + uint256 k; + bytes32 x; + bytes32 y; +} + +contract MathTesterTest is Test { + MathTester public math; + Point[] public points; + + function setUp() public { + math = new MathTester(); + string memory path = "raw_vectors.json"; + string memory json = vm.readFile(path); + bytes memory data = vm.parseJson(json); + points = abi.decode(data, (Point[])); + } + + function test_projectiveAddition() public view { + uint256 len = points.length - 1; + for (uint256 i; i < len;) { + (uint256 x_p, uint256 y_p, uint256 z_p) = math.addProjective( + uint256(points[i].x), + uint256(points[i].y), + uint256(points[i+1].x), + uint256(points[i+1].y) + ); + (x_p, y_p) = math.toAffineProjective(x_p, y_p, z_p); + + (uint256 x_j, uint256 y_j, uint256 z_j) = math.addJacobian( + uint256(points[i].x), + uint256(points[i].y), + uint256(points[i+1].x), + uint256(points[i+1].y) + ); + (x_j, y_j) = math.toAffineJacobian(x_j, y_j, z_j); + + assertEq(x_p, x_j); + assertEq(y_p, y_j); + + unchecked { ++i; } + } + } + + function test_double() public view { + uint256 len = points.length; + for (uint256 i; i < len;) { + (uint256 x_p, uint256 y_p, uint256 z_p) = math.doubleProjective( + uint256(points[i].x), + uint256(points[i].y) + ); + (x_p, y_p) = math.toAffineProjective(x_p, y_p, z_p); + + (uint256 x_j, uint256 y_j, uint256 z_j) = math.doubleJacobian( + uint256(points[i].x), + uint256(points[i].y) + ); + (x_j, y_j) = math.toAffineJacobian(x_j, y_j, z_j); + + assertEq(x_p, x_j); + assertEq(y_p, y_j); + + unchecked { ++i; } + } + } +}