/* * Copyright (C) 2024 Michael Brown . * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation; either version 2 of the * License, or any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA * 02110-1301, USA. * * You can also choose to distribute this program under the terms of * the Unmodified Binary Distribution Licence (as given in the file * COPYING.UBDL), provided that you have satisfied its requirements. */ FILE_LICENCE ( GPL2_OR_LATER_OR_UBDL ); /** @file * * Weierstrass elliptic curves * * The implementation is based upon Algorithm 1 from "Complete * addition formulas for prime order elliptic curves" (Joost Renes, * Craig Costello, and Lejla Batina), available from * * https://www.microsoft.com/en-us/research/wp-content/uploads/2016/06/complete-2.pdf * * The steps within the algorithm have been reordered and temporary * variables shuffled to reduce stack usage, and calculations are * carried out modulo small multiples of the field prime in order to * elide reductions after intermediate addition and subtraction * operations. * * The algorithm is encoded using a bytecode representation, since * this substantially reduces the code size compared to direct * implementation of the big integer operations. */ #include #include /** Big integer register names */ enum weierstrass_register { /* * Read-only registers */ /* Curve constant "a" (for multiply), zero (for add/subtract) */ WEIERSTRASS_a = 0, /* Curve constant "3b" */ WEIERSTRASS_3b, /* Augend (x,y,z) co-ordinates */ WEIERSTRASS_x1, WEIERSTRASS_y1, WEIERSTRASS_z1, /* Addend (x,y,z) co-ordinates */ WEIERSTRASS_x2, WEIERSTRASS_y2, WEIERSTRASS_z2, /* * Read-write registers */ /* Temporary working registers */ WEIERSTRASS_Wt, WEIERSTRASS_Wxy, WEIERSTRASS_Wyz, WEIERSTRASS_Wzx, /* Low half of multiplication product */ WEIERSTRASS_Wp, /* Result (x,y,z) co-ordinates */ WEIERSTRASS_x3, WEIERSTRASS_y3, WEIERSTRASS_z3, /* Number of registers */ WEIERSTRASS_NUM_REGISTERS = 16 }; /** Zero register (for add/subtract operations */ #define WEIERSTRASS_zero WEIERSTRASS_a /** Construct big integer register index */ #define WEIERSTRASS_REGISTER( name ) _C2 ( WEIERSTRASS_, name ) /** Bytecode operation codes */ enum weierstrass_opcode { /** Subtract big integers (and add nothing)*/ WEIERSTRASS_OP_SUB_0N = 0, /** Subtract big integers (and add 2N) */ WEIERSTRASS_OP_SUB_2N = WEIERSTRASS_2N, /** Subtract big integers (and add 4N) */ WEIERSTRASS_OP_SUB_4N = WEIERSTRASS_4N, /** Add big integers */ WEIERSTRASS_OP_ADD, /** Multiply big integers (and perform Montgomery reduction) */ WEIERSTRASS_OP_MUL, }; /** * Define a bytecode operation * * @v opcode Operation code * @v dest Destination big integer register name * @v left Left source big integer register name * @v right Right source big integer register name */ #define WEIERSTRASS_OP( opcode, dest, left, right ) \ ( ( (opcode) << 12 ) | \ ( WEIERSTRASS_REGISTER ( dest ) << 8 ) | \ ( WEIERSTRASS_REGISTER ( left ) << 4 ) | \ ( WEIERSTRASS_REGISTER ( right ) << 0 ) ) /** Extract bytecode operation code */ #define WEIERSTRASS_OPCODE( op ) ( ( (op) >> 12 ) & 0xf ) /** Extract destination big integer register */ #define WEIERSTRASS_DEST( op ) ( ( (op) >> 8 ) & 0xf ) /** Extract left source big integer register */ #define WEIERSTRASS_LEFT( op ) ( ( (op) >> 4 ) & 0xf ) /** Extract right source big integer register */ #define WEIERSTRASS_RIGHT( op ) ( ( (op) >> 0 ) & 0xf ) /** Define a three-argument addition operation */ #define WEIERSTRASS_ADD3( dest, augend, addend ) \ WEIERSTRASS_OP ( WEIERSTRASS_OP_ADD, dest, augend, addend ) /** Define a two-argument addition operation */ #define WEIERSTRASS_ADD2( augend, addend ) \ WEIERSTRASS_ADD3 ( augend, augend, addend ) /** Define a move operation */ #define WEIERSTRASS_MOV( dest, source ) \ WEIERSTRASS_ADD3( dest, source, zero ) /** Define a three-argument subtraction operation */ #define WEIERSTRASS_SUB3( dest, minuend, subtrahend, multiple ) \ WEIERSTRASS_OP ( _C2 ( WEIERSTRASS_OP_SUB_, multiple ), \ dest, minuend, subtrahend ) /** Define a two-argument subtraction operation */ #define WEIERSTRASS_SUB2( minuend, subtrahend, multiple ) \ WEIERSTRASS_SUB3 ( minuend, minuend, subtrahend, multiple ) /** Define a stop operation */ #define WEIERSTRASS_STOP WEIERSTRASS_SUB2 ( zero, zero, 0N ) /** Define a three-argument multiplication operation */ #define WEIERSTRASS_MUL3( dest, multiplicand, multiplier ) \ WEIERSTRASS_OP ( WEIERSTRASS_OP_MUL, dest, multiplicand, multiplier ) /** Define a two-argument multiplication operation */ #define WEIERSTRASS_MUL2( multiplicand, multiplier ) \ WEIERSTRASS_MUL3 ( multiplicand, multiplicand, multiplier ) /** * Initialise curve * * @v curve Weierstrass curve */ static void weierstrass_init ( struct weierstrass_curve *curve ) { unsigned int size = curve->size; bigint_t ( size ) __attribute__ (( may_alias )) *prime = ( ( void * ) curve->prime[0] ); bigint_t ( size ) __attribute__ (( may_alias )) *fermat = ( ( void * ) curve->fermat ); bigint_t ( size ) __attribute__ (( may_alias )) *square = ( ( void * ) curve->square ); bigint_t ( size ) __attribute__ (( may_alias )) *one = ( ( void * ) curve->one ); bigint_t ( size ) __attribute__ (( may_alias )) *a = ( ( void * ) curve->a ); bigint_t ( size ) __attribute__ (( may_alias )) *b3 = ( ( void * ) curve->b3 ); bigint_t ( size ) __attribute__ (( may_alias )) *mont = ( ( void * ) curve->mont[0] ); bigint_t ( size ) __attribute__ (( may_alias )) *temp = ( ( void * ) curve->prime[1] ); bigint_t ( size * 2 ) __attribute__ (( may_alias )) *product = ( ( void * ) temp ); bigint_t ( size ) __attribute__ (( may_alias )) *two = ( ( void * ) temp ); static const uint8_t one_raw[] = { 1 }; static const uint8_t two_raw[] = { 2 }; size_t len = curve->len; unsigned int i; /* Initialise field prime */ bigint_init ( prime, curve->prime_raw, len ); DBGC ( curve, "WEIERSTRASS %s N = %s\n", curve->name, bigint_ntoa ( prime ) ); /* Calculate Montgomery constant R^2 mod N */ bigint_reduce ( prime, square ); DBGC ( curve, "WEIERSTRASS %s R^2 = %s mod N\n", curve->name, bigint_ntoa ( square ) ); /* Calculate constant "3b" */ bigint_init ( b3, curve->b_raw, len ); DBGC ( curve, "WEIERSTRASS %s b = %s\n", curve->name, bigint_ntoa ( b3 ) ); bigint_copy ( b3, a ); bigint_add ( b3, b3 ); bigint_add ( a, b3 ); /* Initialise "a" */ bigint_init ( a, curve->a_raw, len ); DBGC ( curve, "WEIERSTRASS %s a = %s\n", curve->name, bigint_ntoa ( a ) ); /* Initialise "1" */ bigint_init ( one, one_raw, sizeof ( one_raw ) ); /* Convert relevant constants to Montgomery form * * We rely on the fact that the prime multiples have not yet * been calculated, and so can be used as a temporary buffer. */ for ( i = 0 ; i < WEIERSTRASS_NUM_MONT ; i++ ) { static const char *names[] = { " ", " a", "3b" }; bigint_multiply ( &mont[i], square, product ); bigint_montgomery ( prime, product, &mont[i] ); DBGC ( curve, "WEIERSTRASS %s %sR = %s mod N\n", curve->name, names[i], bigint_ntoa ( &mont[i] ) ); } /* Calculate constant "N-2" * * We rely on the fact that the prime multiples have not yet * been calculated, and so can be used as a temporary buffer. */ bigint_copy ( prime, fermat ); bigint_init ( two, two_raw, sizeof ( two_raw ) ); bigint_subtract ( two, fermat ); DBGC ( curve, "WEIERSTRASS %s N-2 = %s\n", curve->name, bigint_ntoa ( fermat ) ); /* Calculate multiples of field prime */ for ( i = 1 ; i < WEIERSTRASS_NUM_MULTIPLES ; i++ ) { bigint_copy ( &prime[ i - 1 ], &prime[i] ); bigint_add ( &prime[i], &prime[i] ); DBGC ( curve, "WEIERSTRASS %s %dN = %s\n", curve->name, ( 1 << i ), bigint_ntoa ( &prime[i] ) ); } } /** * Execute bytecode instruction * * @v curve Weierstrass curve * @v regs Registers * @v size Big integer size * @v op Operation */ static void weierstrass_exec ( const struct weierstrass_curve *curve, void **regs, unsigned int size, unsigned int op ) { const bigint_t ( size ) __attribute__ (( may_alias )) *prime = ( ( const void * ) curve->prime[0] ); bigint_t ( size * 2 ) __attribute__ (( may_alias )) *product = regs[WEIERSTRASS_Wp]; bigint_t ( size ) __attribute__ (( may_alias )) *dest; const bigint_t ( size ) __attribute__ (( may_alias )) *left; const bigint_t ( size ) __attribute__ (( may_alias )) *right; const bigint_t ( size ) __attribute__ (( may_alias )) *addend; const bigint_t ( size ) __attribute__ (( may_alias )) *subtrahend; unsigned int op_code; unsigned int op_dest; unsigned int op_left; unsigned int op_right; /* Decode instruction */ op_code = WEIERSTRASS_OPCODE ( op ); op_dest = WEIERSTRASS_DEST ( op ); op_left = WEIERSTRASS_LEFT ( op ); op_right = WEIERSTRASS_RIGHT ( op ); dest = regs[op_dest]; left = regs[op_left]; right = regs[op_right]; /* Check destination is a writable register */ assert ( op_dest >= WEIERSTRASS_Wt ); /* Handle multiplications */ if ( op_code == WEIERSTRASS_OP_MUL ) { assert ( op_left != WEIERSTRASS_Wp ); assert ( op_right != WEIERSTRASS_Wp ); bigint_multiply ( left, right, product ); bigint_montgomery_relaxed ( prime, product, dest ); DBGCP ( curve, "WEIERSTRASS %s R%d := R%d x R%d = %s\n", curve->name, op_dest, op_left, op_right, bigint_ntoa ( dest ) ); return; } /* Copy left source, if required */ if ( op_dest != op_left ) bigint_copy ( left, dest ); /* Do nothing more if addend/subtrahend is zero */ if ( ! op_right ) { DBGCP ( curve, "WEIERSTRASS %s R%d := R%d = %s\n", curve->name, op_dest, op_left, bigint_ntoa ( dest ) ); return; } /* Determine addend and subtrahend */ addend = NULL; subtrahend = NULL; if ( op_code == WEIERSTRASS_OP_ADD ) { DBGCP ( curve, "WEIERSTRASS %s R%d := R%d + R%d = ", curve->name, op_dest, op_left, op_right ); addend = ( ( const void * ) right ); } else { subtrahend = ( ( const void * ) right ); if ( op_code > WEIERSTRASS_OP_SUB_0N ) { DBGCP ( curve, "WEIERSTRASS %s R%d := R%d - R%d + " "%dN = ", curve->name, op_dest, op_left, op_right, ( 1 << op_code ) ); addend = ( ( const void * ) curve->prime[op_code] ); } else { DBGCP ( curve, "WEIERSTRASS %s R%d := R%d - R%d = ", curve->name, op_dest, op_left, op_right ); } } /* Perform addition and subtraction */ if ( addend ) bigint_add ( addend, dest ); if ( subtrahend ) bigint_subtract ( subtrahend, dest ); DBGCP ( curve, "%s\n", bigint_ntoa ( dest ) ); } /** * Add points on curve * * @v curve Weierstrass curve * @v augend0 Element 0 of point (x1,y1,z1) to be added * @v addend0 Element 0 of point (x2,y2,z2) to be added * @v result0 Element 0 of point (x3,y3,z3) to hold result * * Points are represented in projective coordinates, with all values * in Montgomery form and in the range [0,4N) where N is the field * prime. * * The augend may have the same value as the addend (i.e. this routine * may be used to perform point doubling as well as point addition), * and either or both may be the point at infinity. * * The result may overlap either input, since the inputs are fully * consumed before the result is written. */ static void weierstrass_add_raw ( const struct weierstrass_curve *curve, const bigint_element_t *augend0, const bigint_element_t *addend0, bigint_element_t *result0 ) { unsigned int size = curve->size; const bigint_t ( size ) __attribute__ (( may_alias )) *prime = ( ( const void * ) curve->prime[0] ); const bigint_t ( size ) __attribute__ (( may_alias )) *a = ( ( const void * ) curve->a ); const bigint_t ( size ) __attribute__ (( may_alias )) *b3 = ( ( const void * ) curve->b3 ); const weierstrass_t ( size ) __attribute__ (( may_alias )) *augend = ( ( const void * ) augend0 ); const weierstrass_t ( size ) __attribute__ (( may_alias )) *addend = ( ( const void * ) addend0 ); weierstrass_t ( size ) __attribute__ (( may_alias )) *result = ( ( void * ) result0 ); struct { bigint_t ( size ) Wt; bigint_t ( size ) Wxy; bigint_t ( size ) Wyz; bigint_t ( size ) Wzx; bigint_t ( size * 2 ) Wp; } temp; void *regs[WEIERSTRASS_NUM_REGISTERS]; unsigned int schedule; const uint16_t *op; unsigned int i; /* On entry, we assume that x1, x2, y1, y2, z1, z2 are all in * the range [0,4N). Additions will extend the range. * Subtractions will extend the range (and require an addition * of a suitable multiple of the modulus to ensure that the * result is a positive value). Relaxed Montgomery * multiplications will reduce the range to [0,2N). The * outputs x3, y3, z3 will be in the range [0,4N) and * therefore usable as subsequent inputs. */ static const uint16_t ops[] = { /* [Wxy] Qxy = (x1+y1)*(x2+y2) (mod 2N) */ WEIERSTRASS_ADD3 ( Wt, x1, y1 ), WEIERSTRASS_ADD3 ( Wxy, x2, y2 ), WEIERSTRASS_MUL2 ( Wxy, Wt ), /* [Wyz] Qyz = (y1+z1)*(y2+z2) (mod 2N) */ WEIERSTRASS_ADD3 ( Wt, y1, z1 ), WEIERSTRASS_ADD3 ( Wyz, y2, z2 ), WEIERSTRASS_MUL2 ( Wyz, Wt ), /* [Wzx] Qzx = (z1+x1)*(z2+x2) (mod 2N) */ WEIERSTRASS_ADD3 ( Wt, z1, x1 ), WEIERSTRASS_ADD3 ( Wzx, z2, x2 ), WEIERSTRASS_MUL2 ( Wzx, Wt ), /* [x3] Px = x1*x2 (mod 2N) */ WEIERSTRASS_MUL3 ( x3, x1, x2 ), /* [y3] Py = y1*y2 (mod 2N) */ WEIERSTRASS_MUL3 ( y3, y1, y2 ), /* [z3] Pz = z1*z2 (mod 2N) */ WEIERSTRASS_MUL3 ( z3, z1, z2 ), /* [Wxy] Rxy = Qxy - Px - Py (mod 6N) * = (x1+y1)*(x2+y2) - x1*x2 - y1*y2 (mod 6N) * = x1*y2 + x2*y1 (mod 6N) */ WEIERSTRASS_SUB2 ( Wxy, x3, 0N ), WEIERSTRASS_SUB2 ( Wxy, y3, 4N ), /* [Wyz] Ryz = Qyz - Py - Pz (mod 6N) * = (y1+z1)*(y2+z2) - y1*y2 - z1*z2 (mod 6N) * = y1*z2 + y2*z1 (mod 6N) */ WEIERSTRASS_SUB2 ( Wyz, y3, 0N ), WEIERSTRASS_SUB2 ( Wyz, z3, 4N ), /* [Wzx] Rzx = Qzx - Pz - Px (mod 6N) * = (z1+x1)*(z2+x2) - z1*z2 - x1*x2 (mod 6N) * = x1*z2 + x2*z1 (mod 6N) */ WEIERSTRASS_SUB2 ( Wzx, z3, 0N ), WEIERSTRASS_SUB2 ( Wzx, x3, 4N ), /* [Wt] aRzx = a * Rzx (mod 2N) * = a * (x1*z2 + x2*z1) (mod 2N) */ WEIERSTRASS_MUL3 ( Wt, a, Wzx ), /* [Wp] 3bPz = 3b * Pz (mod 2N) * = 3b*z1*z2 (mod 2N) */ WEIERSTRASS_MUL3 ( Wp, 3b, z3 ), /* [Wp] Sy = aRzx + 3bPz (mod 4N) * = a*(x1*z2 + x2*z1) + 3b*z1*z2 (mod 4N) */ WEIERSTRASS_ADD2 ( Wp, Wt ), /* [Wt] Syz = Py + Sy (mod 6N) * = y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2 (mod 6N) */ WEIERSTRASS_ADD3 ( Wt, y3, Wp ), /* [y3] Sxy = Py - Sy (mod 6N) * = y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2 (mod 6N) */ WEIERSTRASS_SUB2 ( y3, Wp, 4N ), /* [z3] aPz = a * Pz (mod 2N) * = a * z1*z2 (mod 2N) */ WEIERSTRASS_MUL2 ( z3, a ), /* [Wzx] 3bRzx = 3b * Rzx (mod 2N) * = 3b * (x1*z2 + x2*z1) (mod 2N) */ WEIERSTRASS_MUL2 ( Wzx, 3b ), /* [x3] aPzx' = Px - aPz (mod 4N) * = x1*x2 - a*z1*z2 (mod 4N) */ WEIERSTRASS_SUB2 ( x3, z3, 2N ), /* [Wp] Szx = a * aPzx' (mod 2N) * = a * (x1*x2 - a*z1*z2) (mod 2N) * = a*x1*x2 - (a^2)*z1*z2 (mod 2N) */ WEIERSTRASS_MUL3 ( Wp, a, x3 ), /* [x3] Px = aPzx' + aPz (mod 6N) * = x1*x2 - a*z1*z2 + a*z1*z2 (mod 6N) * = x1*x2 (mod 6N) */ WEIERSTRASS_ADD2 ( x3, z3 ), /* [Wzx] Tzx = 3bRzx + Szx (mod 4N) * = a*x1*x2 + 3b*(x1*z2 + x2*z1) - * (a^2)*z1*z2 (mod 4N) */ WEIERSTRASS_ADD2 ( Wzx, Wp ), /* [z3] aPzx = Px + aPz (mod 8N) * = x1*x2 + a*z1*z2 (mod 8N) */ WEIERSTRASS_ADD2 ( z3, x3 ), /* [x3] 2Px = Px + Px (mod 12N) * = 2*x1*x2 (mod 12N) */ WEIERSTRASS_ADD2 ( x3, x3 ), /* [x3] Tyz = 2Px + aPzx (mod 20N) * = 2*x1*x2 + x1*x2 + a*z1*z2 (mod 20N) * = 3*x1*x2 + a*z1*z2 (mod 20N) */ WEIERSTRASS_ADD2 ( x3, z3 ), /* [z3] Syz = Syz (mod 6N) * = y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2 (mod 6N) */ WEIERSTRASS_MOV ( z3, Wt ), /* [Wt] Tyz = Tyz (mod 20N) * = 3*x1*x2 + a*z1*z2 (mod 20N) */ WEIERSTRASS_MOV ( Wt, x3 ), /* [x3] Ux = Rxy * Sxy (mod 2N) * = (x1*y2 + x2*y1) * * (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2) (mod 2N) */ WEIERSTRASS_MUL3 ( x3, Wxy, y3 ), /* [y3] Uy = Syz * Sxy (mod 2N) * = (y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2) * * (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2) (mod 2N) */ WEIERSTRASS_MUL2 ( y3, z3 ), /* [z3] Uz = Ryz * Syz (mod 2N) * = (y1*z2 + y2*z1) * * (y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2) (mod 2N) */ WEIERSTRASS_MUL2 ( z3, Wyz ), /* [Wp] Vx = Ryz * Tzx (mod 2N) * = (y1*z2 + y2*z1) * * (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2) * (mod 2N) */ WEIERSTRASS_MUL3 ( Wp, Wyz, Wzx ), /* [x3] x3 = Ux - Vx (mod 4N) * = ((x1*y2 + x2*y1) * * (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2)) - * ((y1*z2 + y2*z1) * * (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2)) * (mod 4N) */ WEIERSTRASS_SUB2 ( x3, Wp, 2N ), /* [Wp] Vy = Tyz * Tzx (mod 2N) * = (3*x1*x2 + a*z1*z2) * * (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2) * (mod 2N) */ WEIERSTRASS_MUL3 ( Wp, Wt, Wzx ), /* [y3] y3 = Vy + Uy (mod 4N) * = ((3*x1*x2 + a*z1*z2) * * (a*x1*x2 + 3b*(x1*z2 + x2*z1) - (a^2)*z1*z2)) + * ((y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2) * * (y1*y2 - a*(x1*z2 + x2*z1) - 3b*z1*z2)) * (mod 4N) */ WEIERSTRASS_ADD2 ( y3, Wp ), /* [Wp] Vz = Rxy * Tyz (mod 2N) * = (x1*y2 + x2*y1) * (3*x1*x2 + a*z1*z2) (mod 2N) */ WEIERSTRASS_MUL3 ( Wp, Wxy, Wt ), /* [z3] z3 = Uz + Vz (mod 4N) * = ((y1*z2 + y2*z1) * * (y1*y2 + a*(x1*z2 + x2*z1) + 3b*z1*z2)) + * ((x1*y2 + x2*y1) * (3*x1*x2 + a*z1*z2)) * (mod 4N) */ WEIERSTRASS_ADD2 ( z3, Wp ), /* Stop */ WEIERSTRASS_STOP }; /* Initialise register list */ regs[WEIERSTRASS_a] = ( ( void * ) a ); regs[WEIERSTRASS_3b] = ( ( void * ) b3 ); regs[WEIERSTRASS_x1] = ( ( void * ) &augend->x ); regs[WEIERSTRASS_x2] = ( ( void * ) &addend->x ); regs[WEIERSTRASS_x3] = ( ( void * ) &result->x ); regs[WEIERSTRASS_Wt] = &temp; schedule = ( ( ( 1 << WEIERSTRASS_NUM_REGISTERS ) - 1 ) - ( 1 << WEIERSTRASS_a ) - ( 1 << WEIERSTRASS_3b ) - ( 1 << WEIERSTRASS_x1 ) - ( 1 << WEIERSTRASS_x2 ) - ( 1 << WEIERSTRASS_x3 ) - ( 1 << WEIERSTRASS_Wt ) ); for ( i = 0 ; schedule ; i++, schedule >>= 1 ) { if ( schedule & 1 ) regs[i] = ( regs[ i - 1 ] + sizeof ( *prime ) ); } DBGC2 ( curve, "WEIERSTRASS %s augend (%s,", curve->name, bigint_ntoa ( &augend->x ) ); DBGC2 ( curve, "%s,", bigint_ntoa ( &augend->y ) ); DBGC2 ( curve, "%s)\n", bigint_ntoa ( &augend->z ) ); DBGC2 ( curve, "WEIERSTRASS %s addend (%s,", curve->name, bigint_ntoa ( &addend->x ) ); DBGC2 ( curve, "%s,", bigint_ntoa ( &addend->y ) ); DBGC2 ( curve, "%s)\n", bigint_ntoa ( &addend->z ) ); /* Sanity checks */ assert ( regs[WEIERSTRASS_a] == a ); assert ( regs[WEIERSTRASS_3b] == b3 ); assert ( regs[WEIERSTRASS_x1] == &augend->x ); assert ( regs[WEIERSTRASS_y1] == &augend->y ); assert ( regs[WEIERSTRASS_z1] == &augend->z ); assert ( regs[WEIERSTRASS_x2] == &addend->x ); assert ( regs[WEIERSTRASS_y2] == &addend->y ); assert ( regs[WEIERSTRASS_z2] == &addend->z ); assert ( regs[WEIERSTRASS_x3] == &result->x ); assert ( regs[WEIERSTRASS_y3] == &result->y ); assert ( regs[WEIERSTRASS_z3] == &result->z ); assert ( regs[WEIERSTRASS_Wt] == &temp.Wt ); assert ( regs[WEIERSTRASS_Wxy] == &temp.Wxy ); assert ( regs[WEIERSTRASS_Wyz] == &temp.Wyz ); assert ( regs[WEIERSTRASS_Wzx] == &temp.Wzx ); assert ( regs[WEIERSTRASS_Wp] == &temp.Wp ); /* Execute bytecode instruction sequence */ for ( op = ops ; *op != WEIERSTRASS_STOP ; op++ ) weierstrass_exec ( curve, regs, size, *op ); DBGC2 ( curve, "WEIERSTRASS %s result (%s,", curve->name, bigint_ntoa ( &result->x ) ); DBGC2 ( curve, "%s,", bigint_ntoa ( &result->y ) ); DBGC2 ( curve, "%s)\n", bigint_ntoa ( &result->z ) ); } /** * Add points on curve * * @v curve Weierstrass curve * @v augend Point (x1,y1,z1) to be added * @v addend Point (x2,y2,z2) to be added * @v result0 Point (x3,y3,z3) to hold result */ #define weierstrass_add( curve, augend, addend, result ) do { \ weierstrass_add_raw ( (curve), (augend)->all.element, \ (addend)->all.element, \ (result)->all.element ); \ } while ( 0 ) /** * Add points on curve as part of a Montgomery ladder * * @v operand Element 0 of first input operand (may overlap result) * @v result Element 0 of second input operand and result * @v size Number of elements in operands and result * @v ctx Operation context * @v tmp Temporary working space (not used) */ static void weierstrass_add_ladder ( const bigint_element_t *operand0, bigint_element_t *result0, unsigned int size, const void *ctx, void *tmp __unused ) { const struct weierstrass_curve *curve = ctx; const weierstrass_t ( curve->size ) __attribute__ (( may_alias )) *operand = ( ( const void * ) operand0 ); weierstrass_t ( curve->size ) __attribute__ (( may_alias )) *result = ( ( void * ) result0 ); /* Add curve points */ assert ( size == bigint_size ( &operand->all ) ); assert ( size == bigint_size ( &result->all ) ); weierstrass_add ( curve, operand, result, result ); } /** * Verify point is on curve * * @v curve Weierstrass curve * @v point0 Element 0 of point (x,y,z) to be verified * @ret rc Return status code * * As with point addition, points are represented in projective * coordinates, with all values in Montgomery form and in the range * [0,4N) where N is the field prime. */ static int weierstrass_verify_raw ( const struct weierstrass_curve *curve, const bigint_element_t *point0 ) { unsigned int size = curve->size; const bigint_t ( size ) __attribute__ (( may_alias )) *prime = ( ( const void * ) curve->prime[0] ); const bigint_t ( size ) __attribute__ (( may_alias )) *a = ( ( const void * ) curve->a ); const bigint_t ( size ) __attribute__ (( may_alias )) *b3 = ( ( const void * ) curve->b3 ); const weierstrass_t ( size ) __attribute__ (( may_alias )) *point = ( ( const void * ) point0 ); struct { bigint_t ( size ) Wt; bigint_t ( size * 2 ) Wp; } temp; void *regs[WEIERSTRASS_NUM_REGISTERS]; const uint16_t *op; /* Calculate 3*(x^3 + a*x + b - y^2) */ static const uint16_t ops[] = { /* [Wt] Tx = x^2 (mod 2N) */ WEIERSTRASS_MUL3 ( Wt, x1, x1 ), /* [Wt] Txa = Tx + a (mod 3N) * = x^2 + a (mod 3N) */ WEIERSTRASS_MOV ( Wp, a ), WEIERSTRASS_ADD2 ( Wt, Wp ), /* [Wt] Txax = Txa * x (mod 2N) * = (x^2 + a)*x (mod 2N) * = x^3 + a*x (mod 2N) */ WEIERSTRASS_MUL2 ( Wt, x1 ), /* [Wp] Ty = y^2 (mod 2N) */ WEIERSTRASS_MUL3 ( Wp, y1, y1 ), /* [Wt] Txaxy = Txax - Ty (mod 4N) * = x^3 + a*x - y^2 (mod 4N) */ WEIERSTRASS_SUB2 ( Wt, Wp, 2N ), /* [Wp] 2Txaxy = Txaxy + Txaxy (mod 8N) * = 2*(x^3 + a*x - y^2) (mod 8N) */ WEIERSTRASS_ADD3 ( Wp, Wt, Wt ), /* [Wt] 3Txaxy = 2Txaxy + Txaxy (mod 12N) * = 3*(x^3 + a*x - y^2) (mod 12N) */ WEIERSTRASS_ADD2 ( Wt, Wp ), /* [Wt] 3Txaxyb = 3Txaxy + 3b (mod 13N) * = 3*(x^3 + a*x - y^2) + 3b (mod 13N) * = 3*(x^3 + a*x + b - y^2) (mod 13N) */ WEIERSTRASS_ADD2 ( Wt, 3b ), /* Stop */ WEIERSTRASS_STOP }; /* Initialise register list */ regs[WEIERSTRASS_a] = ( ( void * ) a ); regs[WEIERSTRASS_3b] = ( ( void * ) b3 ); regs[WEIERSTRASS_x1] = ( ( void * ) &point->x ); regs[WEIERSTRASS_y1] = ( ( void * ) &point->y ); regs[WEIERSTRASS_Wt] = &temp.Wt; regs[WEIERSTRASS_Wp] = &temp.Wp; /* Execute bytecode instruction sequence */ for ( op = ops ; *op != WEIERSTRASS_STOP ; op++ ) weierstrass_exec ( curve, regs, size, *op ); /* Check that result is zero (modulo the field prime) */ bigint_grow ( &temp.Wt, &temp.Wp ); bigint_montgomery ( prime, &temp.Wp, &temp.Wt ); if ( ! bigint_is_zero ( &temp.Wt ) ) { DBGC ( curve, "WEIERSTRASS %s base point is not on curve\n", curve->name ); return -EINVAL; } return 0; } /** * Verify point is on curve * * @v curve Weierstrass curve * @v point Point (x,y,z) to be verified * @ret rc Return status code */ #define weierstrass_verify( curve, point ) ( { \ weierstrass_verify_raw ( (curve), (point)->all.element ); \ } ) /** * Multiply curve point by scalar * * @v curve Weierstrass curve * @v base Base point (or NULL to use generator) * @v scalar Scalar multiple * @v result Result point to fill in * @ret rc Return status code */ int weierstrass_multiply ( struct weierstrass_curve *curve, const void *base, const void *scalar, void *result ) { unsigned int size = curve->size; size_t len = curve->len; const bigint_t ( size ) __attribute__ (( may_alias )) *prime = ( ( const void * ) curve->prime[0] ); const bigint_t ( size ) __attribute__ (( may_alias )) *prime2 = ( ( const void * ) curve->prime[WEIERSTRASS_2N] ); const bigint_t ( size ) __attribute__ (( may_alias )) *fermat = ( ( const void * ) curve->fermat ); const bigint_t ( size ) __attribute__ (( may_alias )) *square = ( ( const void * ) curve->square ); const bigint_t ( size ) __attribute__ (( may_alias )) *one = ( ( const void * ) curve->one ); struct { union { weierstrass_t ( size ) result; bigint_t ( size * 2 ) product_in; }; union { weierstrass_t ( size ) multiple; bigint_t ( size * 2 ) product_out; }; bigint_t ( bigint_required_size ( len ) ) scalar; } temp; size_t offset; unsigned int i; int rc; /* Initialise curve, if not already done * * The least significant element of the field prime must be * odd, and so the least significant element of the * (initialised) first multiple of the field prime must be * non-zero. */ if ( ! prime2->element[0] ) weierstrass_init ( curve ); /* Use generator if applicable */ if ( ! base ) base = curve->base; /* Convert input to projective coordinates in Montgomery form */ DBGC ( curve, "WEIERSTRASS %s base (", curve->name ); for ( i = 0, offset = 0 ; i < WEIERSTRASS_AXES ; i++, offset += len ) { bigint_init ( &temp.multiple.axis[i], ( base + offset ), len ); DBGC ( curve, "%s%s", ( i ? "," : "" ), bigint_ntoa ( &temp.multiple.axis[i] ) ); bigint_multiply ( &temp.multiple.axis[i], square, &temp.product_in ); bigint_montgomery_relaxed ( prime, &temp.product_in, &temp.multiple.axis[i] ); } bigint_copy ( one, &temp.multiple.z ); DBGC ( curve, ")\n" ); /* Verify point is on curve */ if ( ( rc = weierstrass_verify ( curve, &temp.multiple ) ) != 0 ) return rc; /* Construct identity element (the point at infinity) */ memset ( &temp.result, 0, sizeof ( temp.result ) ); bigint_copy ( one, &temp.result.y ); /* Initialise scalar */ bigint_init ( &temp.scalar, scalar, len ); DBGC ( curve, "WEIERSTRASS %s scalar %s\n", curve->name, bigint_ntoa ( &temp.scalar ) ); /* Perform multiplication via Montgomery ladder */ bigint_ladder ( &temp.result.all, &temp.multiple.all, &temp.scalar, weierstrass_add_ladder, curve, NULL ); /* Invert result Z co-ordinate (via Fermat's little theorem) */ bigint_copy ( one, &temp.multiple.z ); bigint_ladder ( &temp.multiple.z, &temp.result.z, fermat, bigint_mod_exp_ladder, prime, &temp.product_out ); /* Convert result back to affine co-ordinates */ DBGC ( curve, "WEIERSTRASS %s result (", curve->name ); for ( i = 0, offset = 0 ; i < WEIERSTRASS_AXES ; i++, offset += len ) { bigint_multiply ( &temp.result.axis[i], &temp.multiple.z, &temp.product_out ); bigint_montgomery_relaxed ( prime, &temp.product_out, &temp.result.axis[i] ); bigint_grow ( &temp.result.axis[i], &temp.product_out ); bigint_montgomery ( prime, &temp.product_out, &temp.result.axis[i] ); DBGC ( curve, "%s%s", ( i ? "," : "" ), bigint_ntoa ( &temp.result.axis[i] ) ); bigint_done ( &temp.result.axis[i], ( result + offset ), len ); } DBGC ( curve, ")\n" ); return 0; }