Elliptic Curve Operations Within a Solidity ContractElliptic curve cryptography (ECC) operations are actually quite simple, and where we basically have a point add and a point multiply operation. With this we have a base point (\(G\)), and can then operate on this. So, let’s investigate ECC operations on the Ethereum blockchain. The code we will use is: pragma solidity ^0.8.13; import "EllipticCurve.sol"; contract Secp256k1 { uint256 public constant GX = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798; uint256 public constant GY = 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8; uint256 public constant AA = 0; uint256 public constant BB = 7; uint256 public constant PP = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F; function invMod(uint256 val, uint256 p) pure public returns (uint256) { return EllipticCurve.invMod(val,p); } function expMod(uint256 val, uint256 e, uint256 p) pure public returns (uint256) { return EllipticCurve.expMod(val,e,p); } function getY(uint8 prefix, uint256 x) pure public returns (uint256) { return EllipticCurve.deriveY(prefix,x,AA,BB,PP); } function onCurve(uint256 x, uint256 y) pure public returns (bool) { return EllipticCurve.isOnCurve(x,y,AA,BB,PP); } function inverse(uint256 x, uint256 y) pure public returns (uint256, uint256) { return EllipticCurve.ecInv(x,y,PP); } function subtract(uint256 x1, uint256 y1,uint256 x2, uint256 y2 ) pure public returns (uint256, uint256) { return EllipticCurve.ecSub(x1,y1,x2,y2,AA,PP); } function add(uint256 x1, uint256 y1,uint256 x2, uint256 y2 ) pure public returns (uint256, uint256) { return EllipticCurve.ecAdd(x1,y1,x2,y2,AA,PP); } function derivePubKey(uint256 privKey) pure public returns (uint256, uint256) { return EllipticCurve.ecMul(privKey,GX,GY,AA,PP); } } This code can be added to Remix ([here]) and then compile the code: This uses an external contract (EllipticCurve.sol) to implement the core functions. Next we start Ganache as a local blockchain: And then we deploy our smart contract to this blockchain: Now let’s look at a few of the functions we have implemented. Inverse ModWith an inverse of \(a \pmod p\), we calculate a value (\(b\)) that will give: \((a \times b) \pmod p = 1\) In our smart contract we have: function invMod(uint256 val, uint256 p) pure public returns (uint256) { return EllipticCurve.invMod(val,p); } A sample run gives: Thus the inverse of \(19 \pmod{997}\) is 105. We can prove that this is correct with: \((19 \times 105) = 1995 = 1 \pmod {997})\) Modular ExponentialWith a modular exponential, we compute: \(a = b^c \pmod p)\) In our smart contract, we have: function expMod(uint256 val, uint256 e, uint256 p) pure public returns (uint256) { return EllipticCurve.expMod(val,e,p); } If we try a sample run of \(b=65\), \(e=101\), and \(p=997\), we get: This works because: \(a = 65^{101} \pmod {997} = 12713028509888556499558777405 ... 089471876621246337890625 \pmod {997} = 752\) IsOnCurveNow we will try the secp256k1 curve, and where the parameters are: uint256 public constant GX = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798; uint256 public constant GY = 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8; uint256 public constant AA = 0; uint256 public constant BB = 7; uint256 public constant PP =0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F; In the case we have a base point of (\(G_x\), \(G_y\)), a prime number (\(p\)) of 0xFFFFFFFFFFFFFFFFFFFFFF ... FFFFFFFFFFFFFFFEFFFFFC2F (\(p=2^{256}−2^{32}−977\)), and: \(y^2=x^3+7 \pmod p\) First we can determine if a point is on the secp256k1 curve. Let’s take the point from [here]: (89565891926547004231252920425935692360644145829622209833684329913297188986597, 12158399299693830322967808612713398636155367887041628176798871954788371653930) This point represents \(2G\) (or \(G+G\)). Let’s try it in our smart contract: We can see if returns true, while if we add an incorrect point: Recovering the y co-ordinate pointWithin our elliptic curve, we have (\(x,y\)) points, but we do not actually have to store the \(y\) co-ordinate point, and only have to record it’s prefix. Now, let’s use the x co-ordinate value for \(2G\), and see if we can recover the y co-ordinate point: The result is 1215839929969383032 … 887041628176798871954788, and which is correct for \(2G\). Inverse pointWith elliptic curve we can have two points for every valid \(x\) co-ordinate value. Basically the value is the negative value for the y co-ordinate point (and converted to an integer. In the following we determine the inverse of the \(2G\) point we have previously used: The inverse point is thus: 89565891926547004231252920425935692360644145829622209833684329913297188986597, 103633689937622365100603176395974509217114616778598935862658712053120463017733 Multiplying a pointWith elliptic curve methods, we create a random value (n) for our private key, and then derive our public key value with: \(Pub= n.G\) and which is basically \(G\) added to itself \(n\) times. The code is: function derivePubKey(uint256 privKey) pure public returns (uint256,uint256) { return EllipticCurve.ecMul(privKey,GX,GY,AA,PP); } Let’s take a value of n=2, and compute 2.G: We can see we derive the \(2G\) point. Point addingOne of the basic operations in elliptic curve cryptography is the adding of two points: \(P = P_1 + P_2\) Let’s add \(G\) and \(2G\): In this case we get \(3.G\), and we can check here: CodeThe code is here: Addtional codeThe additional code from [here]: // 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 constant private 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, "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, "Invalid 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); } } |