Skip to content

Commit

Permalink
Merge pull request #210
Browse files Browse the repository at this point in the history
2d5a186 Apply effective-affine trick to precomp (Peter Dettman)
4f9791a Effective affine addition in EC multiplication (Peter Dettman)
  • Loading branch information
sipa committed Apr 30, 2015
2 parents 22f60a6 + 2d5a186 commit 729badf
Show file tree
Hide file tree
Showing 6 changed files with 374 additions and 97 deletions.
6 changes: 3 additions & 3 deletions src/bench_internal.c
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}

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

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

Expand Down
10 changes: 5 additions & 5 deletions src/ecmult_gen_impl.h
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand All @@ -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);
Expand Down
166 changes: 112 additions & 54 deletions src/ecmult_impl.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
/** 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_table_precomp_gej_var(secp256k1_gej_t *pre, const secp256k1_gej_t *a, int w) {
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;
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]);

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

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;
/** 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;
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++) {

/* 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(prej);

free(prea);
free(prej);
free(zr);
}

/** The number of entries a table with precomputed multiples needs to have. */
#define ECMULT_TABLE_SIZE(w) (1 << ((w)-2))

/** 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)); \
Expand Down Expand Up @@ -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
{
Expand All @@ -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
}
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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) */
Expand All @@ -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
28 changes: 21 additions & 7 deletions src/group.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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. */
Expand Down
Loading

0 comments on commit 729badf

Please sign in to comment.