Merge pull request #210

2d5a186 Apply effective-affine trick to precomp (Peter Dettman)
4f9791a Effective affine addition in EC multiplication (Peter Dettman)
This commit is contained in:
Pieter Wuille 2015-04-30 21:52:24 +02:00
commit 729badff14
No known key found for this signature in database
GPG Key ID: 57896D2FF8F0B657
6 changed files with 383 additions and 106 deletions

View File

@ -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);
}
}

View File

@ -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);

View File

@ -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

View File

@ -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. */

View File

@ -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

View File

@ -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);