diff --git a/src/bench_internal.c b/src/bench_internal.c index a960549b94f..bbbc7c3fc79 100644 --- a/src/bench_internal.c +++ b/src/bench_internal.c @@ -193,7 +193,7 @@ void bench_group_double_var(void* arg) { bench_inv_t *data = (bench_inv_t*)arg; for (i = 0; i < 200000; i++) { - secp256k1_gej_double_var(&data->gej_x, &data->gej_x); + secp256k1_gej_double_var(&data->gej_x, &data->gej_x, NULL); } } @@ -202,7 +202,7 @@ void bench_group_add_var(void* arg) { bench_inv_t *data = (bench_inv_t*)arg; for (i = 0; i < 200000; i++) { - secp256k1_gej_add_var(&data->gej_x, &data->gej_x, &data->gej_y); + secp256k1_gej_add_var(&data->gej_x, &data->gej_x, &data->gej_y, NULL); } } @@ -220,7 +220,7 @@ void bench_group_add_affine_var(void* arg) { bench_inv_t *data = (bench_inv_t*)arg; for (i = 0; i < 200000; i++) { - secp256k1_gej_add_ge_var(&data->gej_x, &data->gej_x, &data->ge_y); + secp256k1_gej_add_ge_var(&data->gej_x, &data->gej_x, &data->ge_y, NULL); } } diff --git a/src/ecmult_gen_impl.h b/src/ecmult_gen_impl.h index 4697753ac81..a5447de88a8 100644 --- a/src/ecmult_gen_impl.h +++ b/src/ecmult_gen_impl.h @@ -40,7 +40,7 @@ static void secp256k1_ecmult_gen_context_build(secp256k1_ecmult_gen_context_t *c VERIFY_CHECK(secp256k1_ge_set_xo_var(&nums_ge, &nums_x, 0)); secp256k1_gej_set_ge(&nums_gej, &nums_ge); /* Add G to make the bits in x uniformly distributed. */ - secp256k1_gej_add_ge_var(&nums_gej, &nums_gej, &secp256k1_ge_const_g); + secp256k1_gej_add_ge_var(&nums_gej, &nums_gej, &secp256k1_ge_const_g, NULL); } /* compute prec. */ @@ -54,18 +54,18 @@ static void secp256k1_ecmult_gen_context_build(secp256k1_ecmult_gen_context_t *c /* Set precj[j*16 .. j*16+15] to (numsbase, numsbase + gbase, ..., numsbase + 15*gbase). */ precj[j*16] = numsbase; for (i = 1; i < 16; i++) { - secp256k1_gej_add_var(&precj[j*16 + i], &precj[j*16 + i - 1], &gbase); + secp256k1_gej_add_var(&precj[j*16 + i], &precj[j*16 + i - 1], &gbase, NULL); } /* Multiply gbase by 16. */ for (i = 0; i < 4; i++) { - secp256k1_gej_double_var(&gbase, &gbase); + secp256k1_gej_double_var(&gbase, &gbase, NULL); } /* Multiply numbase by 2. */ - secp256k1_gej_double_var(&numsbase, &numsbase); + secp256k1_gej_double_var(&numsbase, &numsbase, NULL); if (j == 62) { /* In the last iteration, numsbase is (1 - 2^j) * nums instead. */ secp256k1_gej_neg(&numsbase, &numsbase); - secp256k1_gej_add_var(&numsbase, &numsbase, &nums_gej); + secp256k1_gej_add_var(&numsbase, &numsbase, &nums_gej, NULL); } } secp256k1_ge_set_all_gej_var(1024, prec, precj); diff --git a/src/ecmult_impl.h b/src/ecmult_impl.h index 1b2856f83df..d6aa2ea7db0 100644 --- a/src/ecmult_impl.h +++ b/src/ecmult_impl.h @@ -24,62 +24,107 @@ #define WINDOW_G 16 #endif -/** Fill a table 'pre' with precomputed odd multiples of a. W determines the size of the table. - * pre will contains the values [1*a,3*a,5*a,...,(2^(w-1)-1)*a], so it needs place for - * 2^(w-2) entries. - * - * There are two versions of this function: - * - secp256k1_ecmult_precomp_wnaf_gej, which operates on group elements in jacobian notation, - * fast to precompute, but slower to use in later additions. - * - secp256k1_ecmult_precomp_wnaf_ge, which operates on group elements in affine notations, - * (much) slower to precompute, but a bit faster to use in later additions. - * To compute a*P + b*G, we use the jacobian version for P, and the affine version for G, as - * G is constant, so it only needs to be done once in advance. - */ -static void secp256k1_ecmult_table_precomp_gej_var(secp256k1_gej_t *pre, const secp256k1_gej_t *a, int w) { - secp256k1_gej_t d; - int i; - pre[0] = *a; - secp256k1_gej_double_var(&d, &pre[0]); - for (i = 1; i < (1 << (w-2)); i++) { - secp256k1_gej_add_var(&pre[i], &d, &pre[i-1]); - } -} - -static void secp256k1_ecmult_table_precomp_ge_storage_var(secp256k1_ge_storage_t *pre, const secp256k1_gej_t *a, int w) { - secp256k1_gej_t d; - int i; - const int table_size = 1 << (w-2); - secp256k1_gej_t *prej = (secp256k1_gej_t *)checked_malloc(sizeof(secp256k1_gej_t) * table_size); - secp256k1_ge_t *prea = (secp256k1_ge_t *)checked_malloc(sizeof(secp256k1_ge_t) * table_size); - prej[0] = *a; - secp256k1_gej_double_var(&d, a); - for (i = 1; i < table_size; i++) { - secp256k1_gej_add_var(&prej[i], &d, &prej[i-1]); - } - secp256k1_ge_set_all_gej_var(table_size, prea, prej); - for (i = 0; i < table_size; i++) { - secp256k1_ge_to_storage(&pre[i], &prea[i]); - } - free(prej); - free(prea); -} - /** The number of entries a table with precomputed multiples needs to have. */ #define ECMULT_TABLE_SIZE(w) (1 << ((w)-2)) +/** Fill a table 'prej' with precomputed odd multiples of a. Prej will contain + * the values [1*a,3*a,...,(2*n-1)*a], so it space for n values. zr[0] will + * contain prej[0].z / a.z. The other zr[i] values = prej[i].z / prej[i-1].z. + * Prej's Z values are undefined, except for the last value. + */ +static void secp256k1_ecmult_odd_multiples_table(int n, secp256k1_gej_t *prej, secp256k1_fe_t *zr, const secp256k1_gej_t *a) { + secp256k1_gej_t d; + secp256k1_ge_t a_ge, d_ge; + int i; + + VERIFY_CHECK(!a->infinity); + + secp256k1_gej_double_var(&d, a, NULL); + + /* + * Perform the additions on an isomorphism where 'd' is affine: drop the z coordinate + * of 'd', and scale the 1P starting value's x/y coordinates without changing its z. + */ + d_ge.x = d.x; + d_ge.y = d.y; + d_ge.infinity = 0; + + secp256k1_ge_set_gej_zinv(&a_ge, a, &d.z); + prej[0].x = a_ge.x; + prej[0].y = a_ge.y; + prej[0].z = a->z; + prej[0].infinity = 0; + + zr[0] = d.z; + for (i = 1; i < n; i++) { + secp256k1_gej_add_ge_var(&prej[i], &prej[i-1], &d_ge, &zr[i]); + } + + /* + * Each point in 'prej' has a z coordinate too small by a factor of 'd.z'. Only + * the final point's z coordinate is actually used though, so just update that. + */ + secp256k1_fe_mul(&prej[n-1].z, &prej[n-1].z, &d.z); +} + +/** Fill a table 'pre' with precomputed odd multiples of a. + * + * There are two versions of this function: + * - secp256k1_ecmult_odd_multiples_table_globalz_windowa which brings its + * resulting point set to a single constant Z denominator, stores the X and Y + * coordinates as ge_storage points in pre, and stores the global Z in rz. + * It only operates on tables sized for WINDOW_A wnaf multiples. + * - secp256k1_ecmult_odd_multiples_table_storage_var, which converts its + * resulting point set to actually affine points, and stores those in pre. + * It operates on tables of any size, but uses heap-allocated temporaries. + * + * To compute a*P + b*G, we compute a table for P using the first function, + * and for G using the second (which requires an inverse, but it only needs to + * happen once). + */ +static void secp256k1_ecmult_odd_multiples_table_globalz_windowa(secp256k1_ge_t *pre, secp256k1_fe_t *globalz, const secp256k1_gej_t *a) { + secp256k1_gej_t prej[ECMULT_TABLE_SIZE(WINDOW_A)]; + secp256k1_fe_t zr[ECMULT_TABLE_SIZE(WINDOW_A)]; + + /* Compute the odd multiples in Jacobian form. */ + secp256k1_ecmult_odd_multiples_table(ECMULT_TABLE_SIZE(WINDOW_A), prej, zr, a); + /* Bring them to the same Z denominator. */ + secp256k1_ge_globalz_set_table_gej(ECMULT_TABLE_SIZE(WINDOW_A), pre, globalz, prej, zr); +} + +static void secp256k1_ecmult_odd_multiples_table_storage_var(int n, secp256k1_ge_storage_t *pre, const secp256k1_gej_t *a) { + secp256k1_gej_t *prej = checked_malloc(sizeof(secp256k1_gej_t) * n); + secp256k1_ge_t *prea = checked_malloc(sizeof(secp256k1_ge_t) * n); + secp256k1_fe_t *zr = checked_malloc(sizeof(secp256k1_fe_t) * n); + int i; + + /* Compute the odd multiples in Jacobian form. */ + secp256k1_ecmult_odd_multiples_table(n, prej, zr, a); + /* Convert them in batch to affine coordinates. */ + secp256k1_ge_set_table_gej_var(n, prea, prej, zr); + /* Convert them to compact storage form. */ + for (i = 0; i < n; i++) { + secp256k1_ge_to_storage(&pre[i], &prea[i]); + } + + free(prea); + free(prej); + free(zr); +} + /** The following two macro retrieves a particular odd multiple from a table * of precomputed multiples. */ -#define ECMULT_TABLE_GET_GEJ(r,pre,n,w) do { \ +#define ECMULT_TABLE_GET_GE(r,pre,n,w) do { \ VERIFY_CHECK(((n) & 1) == 1); \ VERIFY_CHECK((n) >= -((1 << ((w)-1)) - 1)); \ VERIFY_CHECK((n) <= ((1 << ((w)-1)) - 1)); \ if ((n) > 0) { \ *(r) = (pre)[((n)-1)/2]; \ } else { \ - secp256k1_gej_neg((r), &(pre)[(-(n)-1)/2]); \ + secp256k1_ge_neg((r), &(pre)[(-(n)-1)/2]); \ } \ } while(0) + #define ECMULT_TABLE_GET_GE_STORAGE(r,pre,n,w) do { \ VERIFY_CHECK(((n) & 1) == 1); \ VERIFY_CHECK((n) >= -((1 << ((w)-1)) - 1)); \ @@ -112,7 +157,7 @@ static void secp256k1_ecmult_context_build(secp256k1_ecmult_context_t *ctx) { ctx->pre_g = (secp256k1_ge_storage_t (*)[])checked_malloc(sizeof((*ctx->pre_g)[0]) * ECMULT_TABLE_SIZE(WINDOW_G)); /* precompute the tables with odd multiples */ - secp256k1_ecmult_table_precomp_ge_storage_var(*ctx->pre_g, &gj, WINDOW_G); + secp256k1_ecmult_odd_multiples_table_storage_var(ECMULT_TABLE_SIZE(WINDOW_G), *ctx->pre_g, &gj); #ifdef USE_ENDOMORPHISM { @@ -124,9 +169,9 @@ static void secp256k1_ecmult_context_build(secp256k1_ecmult_context_t *ctx) { /* calculate 2^128*generator */ g_128j = gj; for (i = 0; i < 128; i++) { - secp256k1_gej_double_var(&g_128j, &g_128j); + secp256k1_gej_double_var(&g_128j, &g_128j, NULL); } - secp256k1_ecmult_table_precomp_ge_storage_var(*ctx->pre_g_128, &g_128j, WINDOW_G); + secp256k1_ecmult_odd_multiples_table_storage_var(ECMULT_TABLE_SIZE(WINDOW_G), *ctx->pre_g_128, &g_128j); } #endif } @@ -208,11 +253,11 @@ static int secp256k1_ecmult_wnaf(int *wnaf, const secp256k1_scalar_t *a, int w) } static void secp256k1_ecmult(const secp256k1_ecmult_context_t *ctx, secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_scalar_t *na, const secp256k1_scalar_t *ng) { - secp256k1_gej_t tmpj; - secp256k1_gej_t pre_a[ECMULT_TABLE_SIZE(WINDOW_A)]; + secp256k1_ge_t pre_a[ECMULT_TABLE_SIZE(WINDOW_A)]; secp256k1_ge_t tmpa; + secp256k1_fe_t Z; #ifdef USE_ENDOMORPHISM - secp256k1_gej_t pre_a_lam[ECMULT_TABLE_SIZE(WINDOW_A)]; + secp256k1_ge_t pre_a_lam[ECMULT_TABLE_SIZE(WINDOW_A)]; secp256k1_scalar_t na_1, na_lam; /* Splitted G factors. */ secp256k1_scalar_t ng_1, ng_128; @@ -252,12 +297,21 @@ static void secp256k1_ecmult(const secp256k1_ecmult_context_t *ctx, secp256k1_ge bits = bits_na; #endif - /* calculate odd multiples of a */ - secp256k1_ecmult_table_precomp_gej_var(pre_a, a, WINDOW_A); + /* Calculate odd multiples of a. + * All multiples are brought to the same Z 'denominator', which is stored + * in Z. Due to secp256k1' isomorphism we can do all operations pretending + * that the Z coordinate was 1, use affine addition formulae, and correct + * the Z coordinate of the result once at the end. + * The exception is the precomputed G table points, which are actually + * affine. Compared to the base used for other points, they have a Z ratio + * of 1/Z, so we can use secp256k1_gej_add_zinv_var, which uses the same + * isomorphism to efficiently add with a known Z inverse. + */ + secp256k1_ecmult_odd_multiples_table_globalz_windowa(pre_a, &Z, a); #ifdef USE_ENDOMORPHISM for (i = 0; i < ECMULT_TABLE_SIZE(WINDOW_A); i++) { - secp256k1_gej_mul_lambda(&pre_a_lam[i], &pre_a[i]); + secp256k1_ge_mul_lambda(&pre_a_lam[i], &pre_a[i]); } /* split ng into ng_1 and ng_128 (where gn = gn_1 + gn_128*2^128, and gn_1 and gn_128 are ~128 bit) */ @@ -281,37 +335,41 @@ static void secp256k1_ecmult(const secp256k1_ecmult_context_t *ctx, secp256k1_ge secp256k1_gej_set_infinity(r); - for (i = bits-1; i >= 0; i--) { + for (i = bits - 1; i >= 0; i--) { int n; - secp256k1_gej_double_var(r, r); + secp256k1_gej_double_var(r, r, NULL); #ifdef USE_ENDOMORPHISM if (i < bits_na_1 && (n = wnaf_na_1[i])) { - ECMULT_TABLE_GET_GEJ(&tmpj, pre_a, n, WINDOW_A); - secp256k1_gej_add_var(r, r, &tmpj); + ECMULT_TABLE_GET_GE(&tmpa, pre_a, n, WINDOW_A); + secp256k1_gej_add_ge_var(r, r, &tmpa, NULL); } if (i < bits_na_lam && (n = wnaf_na_lam[i])) { - ECMULT_TABLE_GET_GEJ(&tmpj, pre_a_lam, n, WINDOW_A); - secp256k1_gej_add_var(r, r, &tmpj); + ECMULT_TABLE_GET_GE(&tmpa, pre_a_lam, n, WINDOW_A); + secp256k1_gej_add_ge_var(r, r, &tmpa, NULL); } if (i < bits_ng_1 && (n = wnaf_ng_1[i])) { ECMULT_TABLE_GET_GE_STORAGE(&tmpa, *ctx->pre_g, n, WINDOW_G); - secp256k1_gej_add_ge_var(r, r, &tmpa); + secp256k1_gej_add_zinv_var(r, r, &tmpa, &Z); } if (i < bits_ng_128 && (n = wnaf_ng_128[i])) { ECMULT_TABLE_GET_GE_STORAGE(&tmpa, *ctx->pre_g_128, n, WINDOW_G); - secp256k1_gej_add_ge_var(r, r, &tmpa); + secp256k1_gej_add_zinv_var(r, r, &tmpa, &Z); } #else if (i < bits_na && (n = wnaf_na[i])) { - ECMULT_TABLE_GET_GEJ(&tmpj, pre_a, n, WINDOW_A); - secp256k1_gej_add_var(r, r, &tmpj); + ECMULT_TABLE_GET_GE(&tmpa, pre_a, n, WINDOW_A); + secp256k1_gej_add_ge_var(r, r, &tmpa, NULL); } if (i < bits_ng && (n = wnaf_ng[i])) { ECMULT_TABLE_GET_GE_STORAGE(&tmpa, *ctx->pre_g, n, WINDOW_G); - secp256k1_gej_add_ge_var(r, r, &tmpa); + secp256k1_gej_add_zinv_var(r, r, &tmpa, &Z); } #endif } + + if (!r->infinity) { + secp256k1_fe_mul(&r->z, &r->z, &Z); + } } #endif diff --git a/src/group.h b/src/group.h index 0b08b3b9910..1d9ef9d2aa7 100644 --- a/src/group.h +++ b/src/group.h @@ -62,6 +62,17 @@ static void secp256k1_ge_set_gej(secp256k1_ge_t *r, secp256k1_gej_t *a); /** Set a batch of group elements equal to the inputs given in jacobian coordinates */ static void secp256k1_ge_set_all_gej_var(size_t len, secp256k1_ge_t *r, const secp256k1_gej_t *a); +/** Set a batch of group elements equal to the inputs given in jacobian + * coordinates (with known z-ratios). zr must contain the known z-ratios such + * that mul(a[i].z, zr[i+1]) == a[i+1].z. zr[0] is ignored. */ +static void secp256k1_ge_set_table_gej_var(size_t len, secp256k1_ge_t *r, const secp256k1_gej_t *a, const secp256k1_fe_t *zr); + +/** Bring a batch inputs given in jacobian coordinates (with known z-ratios) to + * the same global z "denominator". zr must contain the known z-ratios such + * that mul(a[i].z, zr[i+1]) == a[i+1].z. zr[0] is ignored. The x and y + * coordinates of the result are stored in r, the common z coordinate is + * stored in globalz. */ +static void secp256k1_ge_globalz_set_table_gej(size_t len, secp256k1_ge_t *r, secp256k1_fe_t *globalz, const secp256k1_gej_t *a, const secp256k1_fe_t *zr); /** Set a group element (jacobian) equal to the point at infinity. */ static void secp256k1_gej_set_infinity(secp256k1_gej_t *r); @@ -81,23 +92,26 @@ static void secp256k1_gej_neg(secp256k1_gej_t *r, const secp256k1_gej_t *a); /** Check whether a group element is the point at infinity. */ static int secp256k1_gej_is_infinity(const secp256k1_gej_t *a); -/** Set r equal to the double of a. */ -static void secp256k1_gej_double_var(secp256k1_gej_t *r, const secp256k1_gej_t *a); +/** Set r equal to the double of a. If rzr is not-NULL, r->z = a->z * *rzr (where infinity means an implicit z = 0). */ +static void secp256k1_gej_double_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, secp256k1_fe_t *rzr); -/** Set r equal to the sum of a and b. */ -static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_gej_t *b); +/** Set r equal to the sum of a and b. If rzr is non-NULL, r->z = a->z * *rzr (a cannot be infinity in that case). */ +static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_gej_t *b, secp256k1_fe_t *rzr); /** Set r equal to the sum of a and b (with b given in affine coordinates, and not infinity). */ static void secp256k1_gej_add_ge(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b); /** Set r equal to the sum of a and b (with b given in affine coordinates). This is more efficient than secp256k1_gej_add_var. It is identical to secp256k1_gej_add_ge but without constant-time - guarantee, and b is allowed to be infinity. */ -static void secp256k1_gej_add_ge_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b); + guarantee, and b is allowed to be infinity. If rzr is non-NULL, r->z = a->z * *rzr (a cannot be infinity in that case). */ +static void secp256k1_gej_add_ge_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b, secp256k1_fe_t *rzr); + +/** Set r equal to the sum of a and b (with the inverse of b's Z coordinate passed as bzinv). */ +static void secp256k1_gej_add_zinv_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b, const secp256k1_fe_t *bzinv); #ifdef USE_ENDOMORPHISM /** Set r to be equal to lambda times a, where lambda is chosen in a way such that this is very fast. */ -static void secp256k1_gej_mul_lambda(secp256k1_gej_t *r, const secp256k1_gej_t *a); +static void secp256k1_ge_mul_lambda(secp256k1_ge_t *r, const secp256k1_ge_t *a); #endif /** Clear a secp256k1_gej_t to prevent leaking sensitive information. */ diff --git a/src/group_impl.h b/src/group_impl.h index 0f64576fbb5..45bafa6623d 100644 --- a/src/group_impl.h +++ b/src/group_impl.h @@ -23,6 +23,16 @@ static const secp256k1_ge_t secp256k1_ge_const_g = SECP256K1_GE_CONST( 0xFD17B448UL, 0xA6855419UL, 0x9C47D08FUL, 0xFB10D4B8UL ); +static void secp256k1_ge_set_gej_zinv(secp256k1_ge_t *r, const secp256k1_gej_t *a, const secp256k1_fe_t *zi) { + secp256k1_fe_t zi2; + secp256k1_fe_t zi3; + secp256k1_fe_sqr(&zi2, zi); + secp256k1_fe_mul(&zi3, &zi2, zi); + secp256k1_fe_mul(&r->x, &a->x, &zi2); + secp256k1_fe_mul(&r->y, &a->y, &zi3); + r->infinity = a->infinity; +} + static void secp256k1_ge_set_infinity(secp256k1_ge_t *r) { r->infinity = 1; } @@ -92,17 +102,55 @@ static void secp256k1_ge_set_all_gej_var(size_t len, secp256k1_ge_t *r, const se for (i = 0; i < len; i++) { r[i].infinity = a[i].infinity; if (!a[i].infinity) { - secp256k1_fe_t zi2, zi3; - secp256k1_fe_t *zi = &azi[count++]; - secp256k1_fe_sqr(&zi2, zi); - secp256k1_fe_mul(&zi3, &zi2, zi); - secp256k1_fe_mul(&r[i].x, &a[i].x, &zi2); - secp256k1_fe_mul(&r[i].y, &a[i].y, &zi3); + secp256k1_ge_set_gej_zinv(&r[i], &a[i], &azi[count++]); } } free(azi); } +static void secp256k1_ge_set_table_gej_var(size_t len, secp256k1_ge_t *r, const secp256k1_gej_t *a, const secp256k1_fe_t *zr) { + size_t i = len - 1; + secp256k1_fe_t zi; + + if (len < 1) + return; + + /* Compute the inverse of the last z coordinate, and use it to compute the last affine output. */ + secp256k1_fe_inv(&zi, &a[i].z); + secp256k1_ge_set_gej_zinv(&r[i], &a[i], &zi); + + /* Work out way backwards, using the z-ratios to scale the x/y values. */ + while (i > 0) { + secp256k1_fe_mul(&zi, &zi, &zr[i]); + i--; + secp256k1_ge_set_gej_zinv(&r[i], &a[i], &zi); + } +} + +static void secp256k1_ge_globalz_set_table_gej(size_t len, secp256k1_ge_t *r, secp256k1_fe_t *globalz, const secp256k1_gej_t *a, const secp256k1_fe_t *zr) { + size_t i = len - 1; + secp256k1_fe_t zs; + + if (len < 1) + return; + + /* The z of the final point gives us the "global Z" for the table. */ + r[i].x = a[i].x; + r[i].y = a[i].y; + *globalz = a[i].z; + r[i].infinity = 0; + zs = zr[i]; + + /* Work our way backwards, using the z-ratios to scale the x/y values. */ + while (i > 0) { + if (i != len - 1) { + secp256k1_fe_mul(&zs, &zs, &zr[i]); + } + i--; + secp256k1_ge_set_gej_zinv(&r[i], &a[i], &zs); + } +} + static void secp256k1_gej_set_infinity(secp256k1_gej_t *r) { r->infinity = 1; secp256k1_fe_set_int(&r->x, 0); @@ -210,7 +258,7 @@ static int secp256k1_ge_is_valid_var(const secp256k1_ge_t *a) { return secp256k1_fe_equal_var(&y2, &x3); } -static void secp256k1_gej_double_var(secp256k1_gej_t *r, const secp256k1_gej_t *a) { +static void secp256k1_gej_double_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, secp256k1_fe_t *rzr) { /* Operations: 3 mul, 4 sqr, 0 normalize, 12 mul_int/add/negate */ secp256k1_fe_t t1,t2,t3,t4; /** For secp256k1, 2Q is infinity if and only if Q is infinity. This is because if 2Q = infinity, @@ -219,9 +267,18 @@ static void secp256k1_gej_double_var(secp256k1_gej_t *r, const secp256k1_gej_t * */ r->infinity = a->infinity; if (r->infinity) { + if (rzr) { + secp256k1_fe_set_int(rzr, 1); + } return; } + if (rzr) { + *rzr = a->y; + secp256k1_fe_normalize_weak(rzr); + secp256k1_fe_mul_int(rzr, 2); + } + secp256k1_fe_mul(&r->z, &a->z, &a->y); secp256k1_fe_mul_int(&r->z, 2); /* Z' = 2*Y*Z (2) */ secp256k1_fe_sqr(&t1, &a->x); @@ -244,17 +301,24 @@ static void secp256k1_gej_double_var(secp256k1_gej_t *r, const secp256k1_gej_t * secp256k1_fe_add(&r->y, &t2); /* Y' = 36*X^3*Y^2 - 27*X^6 - 8*Y^4 (4) */ } -static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_gej_t *b) { +static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_gej_t *b, secp256k1_fe_t *rzr) { /* Operations: 12 mul, 4 sqr, 2 normalize, 12 mul_int/add/negate */ secp256k1_fe_t z22, z12, u1, u2, s1, s2, h, i, i2, h2, h3, t; + if (a->infinity) { + VERIFY_CHECK(rzr == NULL); *r = *b; return; } + if (b->infinity) { + if (rzr) { + secp256k1_fe_set_int(rzr, 1); + } *r = *a; return; } + r->infinity = 0; secp256k1_fe_sqr(&z22, &b->z); secp256k1_fe_sqr(&z12, &a->z); @@ -266,8 +330,11 @@ static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, secp256k1_fe_negate(&i, &s1, 1); secp256k1_fe_add(&i, &s2); if (secp256k1_fe_normalizes_to_zero_var(&h)) { if (secp256k1_fe_normalizes_to_zero_var(&i)) { - secp256k1_gej_double_var(r, a); + secp256k1_gej_double_var(r, a, rzr); } else { + if (rzr) { + secp256k1_fe_set_int(rzr, 0); + } r->infinity = 1; } return; @@ -275,7 +342,11 @@ static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, secp256k1_fe_sqr(&i2, &i); secp256k1_fe_sqr(&h2, &h); secp256k1_fe_mul(&h3, &h, &h2); - secp256k1_fe_mul(&r->z, &a->z, &b->z); secp256k1_fe_mul(&r->z, &r->z, &h); + secp256k1_fe_mul(&h, &h, &b->z); + if (rzr) { + *rzr = h; + } + secp256k1_fe_mul(&r->z, &a->z, &h); secp256k1_fe_mul(&t, &u1, &h2); r->x = t; secp256k1_fe_mul_int(&r->x, 2); secp256k1_fe_add(&r->x, &h3); secp256k1_fe_negate(&r->x, &r->x, 3); secp256k1_fe_add(&r->x, &i2); secp256k1_fe_negate(&r->y, &r->x, 5); secp256k1_fe_add(&r->y, &t); secp256k1_fe_mul(&r->y, &r->y, &i); @@ -283,21 +354,23 @@ static void secp256k1_gej_add_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, secp256k1_fe_add(&r->y, &h3); } -static void secp256k1_gej_add_ge_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b) { +static void secp256k1_gej_add_ge_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b, secp256k1_fe_t *rzr) { /* 8 mul, 3 sqr, 4 normalize, 12 mul_int/add/negate */ secp256k1_fe_t z12, u1, u2, s1, s2, h, i, i2, h2, h3, t; if (a->infinity) { - r->infinity = b->infinity; - r->x = b->x; - r->y = b->y; - secp256k1_fe_set_int(&r->z, 1); + VERIFY_CHECK(rzr == NULL); + secp256k1_gej_set_ge(r, b); return; } if (b->infinity) { + if (rzr) { + secp256k1_fe_set_int(rzr, 1); + } *r = *a; return; } r->infinity = 0; + secp256k1_fe_sqr(&z12, &a->z); u1 = a->x; secp256k1_fe_normalize_weak(&u1); secp256k1_fe_mul(&u2, &b->x, &z12); @@ -307,7 +380,69 @@ static void secp256k1_gej_add_ge_var(secp256k1_gej_t *r, const secp256k1_gej_t * secp256k1_fe_negate(&i, &s1, 1); secp256k1_fe_add(&i, &s2); if (secp256k1_fe_normalizes_to_zero_var(&h)) { if (secp256k1_fe_normalizes_to_zero_var(&i)) { - secp256k1_gej_double_var(r, a); + secp256k1_gej_double_var(r, a, rzr); + } else { + if (rzr) { + secp256k1_fe_set_int(rzr, 0); + } + r->infinity = 1; + } + return; + } + secp256k1_fe_sqr(&i2, &i); + secp256k1_fe_sqr(&h2, &h); + secp256k1_fe_mul(&h3, &h, &h2); + if (rzr) { + *rzr = h; + } + secp256k1_fe_mul(&r->z, &a->z, &h); + secp256k1_fe_mul(&t, &u1, &h2); + r->x = t; secp256k1_fe_mul_int(&r->x, 2); secp256k1_fe_add(&r->x, &h3); secp256k1_fe_negate(&r->x, &r->x, 3); secp256k1_fe_add(&r->x, &i2); + secp256k1_fe_negate(&r->y, &r->x, 5); secp256k1_fe_add(&r->y, &t); secp256k1_fe_mul(&r->y, &r->y, &i); + secp256k1_fe_mul(&h3, &h3, &s1); secp256k1_fe_negate(&h3, &h3, 1); + secp256k1_fe_add(&r->y, &h3); +} + +static void secp256k1_gej_add_zinv_var(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b, const secp256k1_fe_t *bzinv) { + /* 9 mul, 3 sqr, 4 normalize, 12 mul_int/add/negate */ + secp256k1_fe_t az, z12, u1, u2, s1, s2, h, i, i2, h2, h3, t; + + if (b->infinity) { + *r = *a; + return; + } + if (a->infinity) { + secp256k1_fe_t bzinv2, bzinv3; + r->infinity = b->infinity; + secp256k1_fe_sqr(&bzinv2, bzinv); + secp256k1_fe_mul(&bzinv3, &bzinv2, bzinv); + secp256k1_fe_mul(&r->x, &b->x, &bzinv2); + secp256k1_fe_mul(&r->y, &b->y, &bzinv3); + secp256k1_fe_set_int(&r->z, 1); + return; + } + r->infinity = 0; + + /** We need to calculate (rx,ry,rz) = (ax,ay,az) + (bx,by,1/bzinv). Due to + * secp256k1's isomorphism we can multiply the Z coordinates on both sides + * by bzinv, and get: (rx,ry,rz*bzinv) = (ax,ay,az*bzinv) + (bx,by,1). + * This means that (rx,ry,rz) can be calculated as + * (ax,ay,az*bzinv) + (bx,by,1), when not applying the bzinv factor to rz. + * The variable az below holds the modified Z coordinate for a, which is used + * for the computation of rx and ry, but not for rz. + */ + secp256k1_fe_mul(&az, &a->z, bzinv); + + secp256k1_fe_sqr(&z12, &az); + u1 = a->x; secp256k1_fe_normalize_weak(&u1); + secp256k1_fe_mul(&u2, &b->x, &z12); + s1 = a->y; secp256k1_fe_normalize_weak(&s1); + secp256k1_fe_mul(&s2, &b->y, &z12); secp256k1_fe_mul(&s2, &s2, &az); + secp256k1_fe_negate(&h, &u1, 1); secp256k1_fe_add(&h, &u2); + secp256k1_fe_negate(&i, &s1, 1); secp256k1_fe_add(&i, &s2); + if (secp256k1_fe_normalizes_to_zero_var(&h)) { + if (secp256k1_fe_normalizes_to_zero_var(&i)) { + secp256k1_gej_double_var(r, a, NULL); } else { r->infinity = 1; } @@ -324,6 +459,7 @@ static void secp256k1_gej_add_ge_var(secp256k1_gej_t *r, const secp256k1_gej_t * secp256k1_fe_add(&r->y, &h3); } + static void secp256k1_gej_add_ge(secp256k1_gej_t *r, const secp256k1_gej_t *a, const secp256k1_ge_t *b) { /* Operations: 7 mul, 5 sqr, 5 normalize, 17 mul_int/add/negate/cmov */ static const secp256k1_fe_t fe_1 = SECP256K1_FE_CONST(0, 0, 0, 0, 0, 0, 0, 1); @@ -430,7 +566,7 @@ static SECP256K1_INLINE void secp256k1_ge_storage_cmov(secp256k1_ge_storage_t *r } #ifdef USE_ENDOMORPHISM -static void secp256k1_gej_mul_lambda(secp256k1_gej_t *r, const secp256k1_gej_t *a) { +static void secp256k1_ge_mul_lambda(secp256k1_ge_t *r, const secp256k1_ge_t *a) { static const secp256k1_fe_t beta = SECP256K1_FE_CONST( 0x7ae96a2bul, 0x657c0710ul, 0x6e64479eul, 0xac3434e9ul, 0x9cf04975ul, 0x12f58995ul, 0xc1396c28ul, 0x719501eeul diff --git a/src/tests.c b/src/tests.c index d0e05057f29..fb8ce76924a 100644 --- a/src/tests.c +++ b/src/tests.c @@ -966,6 +966,10 @@ void test_ge(void) { */ secp256k1_ge_t *ge = (secp256k1_ge_t *)malloc(sizeof(secp256k1_ge_t) * (1 + 4 * runs)); secp256k1_gej_t *gej = (secp256k1_gej_t *)malloc(sizeof(secp256k1_gej_t) * (1 + 4 * runs)); + secp256k1_fe_t *zinv = (secp256k1_fe_t *)malloc(sizeof(secp256k1_fe_t) * (1 + 4 * runs)); + secp256k1_fe_t zf; + secp256k1_fe_t zfi2, zfi3; + secp256k1_gej_set_infinity(&gej[0]); secp256k1_ge_clear(&ge[0]); secp256k1_ge_set_gej_var(&ge[0], &gej[0]); @@ -990,18 +994,65 @@ void test_ge(void) { } } + /* Compute z inverses. */ + { + secp256k1_fe_t *zs = malloc(sizeof(secp256k1_fe_t) * (1 + 4 * runs)); + for (i = 0; i < 4 * runs + 1; i++) { + if (i == 0) { + /* The point at infinity does not have a meaningful z inverse. Any should do. */ + do { + random_field_element_test(&zs[i]); + } while(secp256k1_fe_is_zero(&zs[i])); + } else { + zs[i] = gej[i].z; + } + } + secp256k1_fe_inv_all_var(4 * runs + 1, zinv, zs); + free(zs); + } + + /* Generate random zf, and zfi2 = 1/zf^2, zfi3 = 1/zf^3 */ + do { + random_field_element_test(&zf); + } while(secp256k1_fe_is_zero(&zf)); + random_field_element_magnitude(&zf); + secp256k1_fe_inv_var(&zfi3, &zf); + secp256k1_fe_sqr(&zfi2, &zfi3); + secp256k1_fe_mul(&zfi3, &zfi3, &zfi2); + for (i1 = 0; i1 < 1 + 4 * runs; i1++) { int i2; for (i2 = 0; i2 < 1 + 4 * runs; i2++) { /* Compute reference result using gej + gej (var). */ secp256k1_gej_t refj, resj; secp256k1_ge_t ref; - secp256k1_gej_add_var(&refj, &gej[i1], &gej[i2]); + secp256k1_fe_t zr; + secp256k1_gej_add_var(&refj, &gej[i1], &gej[i2], secp256k1_gej_is_infinity(&gej[i1]) ? NULL : &zr); + /* Check Z ratio. */ + if (!secp256k1_gej_is_infinity(&gej[i1]) && !secp256k1_gej_is_infinity(&refj)) { + secp256k1_fe_t zrz; secp256k1_fe_mul(&zrz, &zr, &gej[i1].z); + CHECK(secp256k1_fe_equal_var(&zrz, &refj.z)); + } secp256k1_ge_set_gej_var(&ref, &refj); - /* Test gej + ge (var). */ - secp256k1_gej_add_ge_var(&resj, &gej[i1], &ge[i2]); + /* Test gej + ge with Z ratio result (var). */ + secp256k1_gej_add_ge_var(&resj, &gej[i1], &ge[i2], secp256k1_gej_is_infinity(&gej[i1]) ? NULL : &zr); ge_equals_gej(&ref, &resj); + if (!secp256k1_gej_is_infinity(&gej[i1]) && !secp256k1_gej_is_infinity(&resj)) { + secp256k1_fe_t zrz; secp256k1_fe_mul(&zrz, &zr, &gej[i1].z); + CHECK(secp256k1_fe_equal_var(&zrz, &resj.z)); + } + + /* Test gej + ge (var, with additional Z factor). */ + { + secp256k1_ge_t ge2_zfi = ge[i2]; /* the second term with x and y rescaled for z = 1/zf */ + secp256k1_fe_mul(&ge2_zfi.x, &ge2_zfi.x, &zfi2); + secp256k1_fe_mul(&ge2_zfi.y, &ge2_zfi.y, &zfi3); + random_field_element_magnitude(&ge2_zfi.x); + random_field_element_magnitude(&ge2_zfi.y); + secp256k1_gej_add_zinv_var(&resj, &gej[i1], &ge2_zfi, &zf); + ge_equals_gej(&ref, &resj); + } /* Test gej + ge (const). */ if (i2 != 0) { @@ -1012,10 +1063,15 @@ void test_ge(void) { /* Test doubling (var). */ if ((i1 == 0 && i2 == 0) || ((i1 + 3)/4 == (i2 + 3)/4 && ((i1 + 3)%4)/2 == ((i2 + 3)%4)/2)) { - /* Normal doubling. */ - secp256k1_gej_double_var(&resj, &gej[i1]); + secp256k1_fe_t zr2; + /* Normal doubling with Z ratio result. */ + secp256k1_gej_double_var(&resj, &gej[i1], &zr2); ge_equals_gej(&ref, &resj); - secp256k1_gej_double_var(&resj, &gej[i2]); + /* Check Z ratio. */ + secp256k1_fe_mul(&zr2, &zr2, &gej[i1].z); + CHECK(secp256k1_fe_equal_var(&zr2, &resj.z)); + /* Normal doubling. */ + secp256k1_gej_double_var(&resj, &gej[i2], NULL); ge_equals_gej(&ref, &resj); } @@ -1054,27 +1110,40 @@ void test_ge(void) { } } for (i = 0; i < 4 * runs + 1; i++) { - secp256k1_gej_add_var(&sum, &sum, &gej_shuffled[i]); + secp256k1_gej_add_var(&sum, &sum, &gej_shuffled[i], NULL); } CHECK(secp256k1_gej_is_infinity(&sum)); free(gej_shuffled); } - /* Test batch gej -> ge conversion. */ + /* Test batch gej -> ge conversion with and without known z ratios. */ { + secp256k1_fe_t *zr = (secp256k1_fe_t *)malloc((4 * runs + 1) * sizeof(secp256k1_fe_t)); + secp256k1_ge_t *ge_set_table = (secp256k1_ge_t *)malloc((4 * runs + 1) * sizeof(secp256k1_ge_t)); secp256k1_ge_t *ge_set_all = (secp256k1_ge_t *)malloc((4 * runs + 1) * sizeof(secp256k1_ge_t)); + for (i = 0; i < 4 * runs + 1; i++) { + /* Compute gej[i + 1].z / gez[i].z (with gej[n].z taken to be 1). */ + if (i < 4 * runs) { + secp256k1_fe_mul(&zr[i + 1], &zinv[i], &gej[i + 1].z); + } + } + secp256k1_ge_set_table_gej_var(4 * runs + 1, ge_set_table, gej, zr); secp256k1_ge_set_all_gej_var(4 * runs + 1, ge_set_all, gej); for (i = 0; i < 4 * runs + 1; i++) { secp256k1_fe_t s; random_fe_non_zero(&s); secp256k1_gej_rescale(&gej[i], &s); + ge_equals_gej(&ge_set_table[i], &gej[i]); ge_equals_gej(&ge_set_all[i], &gej[i]); } + free(ge_set_table); free(ge_set_all); + free(zr); } free(ge); free(gej); + free(zinv); } void run_ge(void) { @@ -1139,14 +1208,14 @@ void run_ecmult_chain(void) { ); secp256k1_gej_neg(&rp, &rp); - secp256k1_gej_add_var(&rp, &rp, &x); + secp256k1_gej_add_var(&rp, &rp, &x, NULL); CHECK(secp256k1_gej_is_infinity(&rp)); } } /* redo the computation, but directly with the resulting ae and ge coefficients: */ secp256k1_ecmult(&ctx->ecmult_ctx, &x2, &a, &ae, &ge); secp256k1_gej_neg(&x2, &x2); - secp256k1_gej_add_var(&x2, &x2, &x); + secp256k1_gej_add_var(&x2, &x2, &x, NULL); CHECK(secp256k1_gej_is_infinity(&x2)); } @@ -1162,7 +1231,7 @@ void test_point_times_order(const secp256k1_gej_t *point) { secp256k1_scalar_negate(&nx, &x); secp256k1_ecmult(&ctx->ecmult_ctx, &res1, point, &x, &x); /* calc res1 = x * point + x * G; */ secp256k1_ecmult(&ctx->ecmult_ctx, &res2, point, &nx, &nx); /* calc res2 = (order - x) * point + (order - x) * G; */ - secp256k1_gej_add_var(&res1, &res1, &res2); + secp256k1_gej_add_var(&res1, &res1, &res2, NULL); CHECK(secp256k1_gej_is_infinity(&res1)); CHECK(secp256k1_gej_is_valid_var(&res1) == 0); secp256k1_ge_set_gej(&res3, &res1);