Use batch inversion in G precomputation

src/ecmult_impl.h

@@ -36,14 +36,15 @@ void static secp256k1_ecmult_table_precomp_gej(secp256k1_gej_t *pre, const secp2
         secp256k1_gej_add(&pre[i], &d, &pre[i-1]);
 }
 
-void static secp256k1_ecmult_table_precomp_ge(secp256k1_ge_t *pre, const secp256k1_ge_t *a, int w) {
-    pre[0] = *a;
-    secp256k1_gej_t x; secp256k1_gej_set_ge(&x, a);
-    secp256k1_gej_t d; secp256k1_gej_double(&d, &x);
-    for (int i=1; i<(1 << (w-2)); i++) {
-        secp256k1_gej_add_ge(&x, &d, &pre[i-1]);
-        secp256k1_ge_set_gej(&pre[i], &x);
+void static secp256k1_ecmult_table_precomp_ge(secp256k1_ge_t *pre, const secp256k1_gej_t *a, int w) {
+    const int table_size = 1 << (w-2);
+    secp256k1_gej_t prej[table_size];
+    prej[0] = *a;
+    secp256k1_gej_t d; secp256k1_gej_double(&d, a);
+    for (int i=1; i<table_size; i++) {
+        secp256k1_gej_add(&prej[i], &d, &prej[i-1]);
     }
+    secp256k1_ge_set_all_gej(table_size, pre, prej);
 }
 
 /** The number of entries a table with precomputed multiples needs to have. */
@@ -93,37 +94,45 @@ static void secp256k1_ecmult_start(void) {
 
     // get the generator
     const secp256k1_ge_t *g = &secp256k1_ge_consts->g;
+    secp256k1_gej_t gj; secp256k1_gej_set_ge(&gj, g);
 
     // calculate 2^128*generator
-    secp256k1_gej_t g_128j; secp256k1_gej_set_ge(&g_128j, g);
+    secp256k1_gej_t g_128j = gj;
     for (int i=0; i<128; i++)
         secp256k1_gej_double(&g_128j, &g_128j);
-    secp256k1_ge_t g_128; secp256k1_ge_set_gej(&g_128, &g_128j);
 
     // precompute the tables with odd multiples
-    secp256k1_ecmult_table_precomp_ge(ret->pre_g, g, WINDOW_G);
-    secp256k1_ecmult_table_precomp_ge(ret->pre_g_128, &g_128, WINDOW_G);
+    secp256k1_ecmult_table_precomp_ge(ret->pre_g, &gj, WINDOW_G);
+    secp256k1_ecmult_table_precomp_ge(ret->pre_g_128, &g_128j, WINDOW_G);
 
     // compute prec and fin
-    secp256k1_gej_t gg; secp256k1_gej_set_ge(&gg, g);
-    secp256k1_ge_t ggn; ggn = *g;
-    secp256k1_ge_t ad = *g;
+    secp256k1_gej_t tj[961];
+    int pos = 0;
     secp256k1_gej_t fn; secp256k1_gej_set_infinity(&fn);
     for (int j=0; j<64; j++) {
-        for (int k=0; k<sizeof(secp256k1_ge_t); k++)
-            ret->prec[j][k][0] = ((unsigned char*)(&ggn))[k];
-        secp256k1_gej_add(&fn, &fn, &gg);
+        secp256k1_gej_add(&fn, &fn, &gj);
+        secp256k1_gej_t adj = gj;
         for (int i=1; i<16; i++) {
-            secp256k1_gej_add_ge(&gg, &gg, &ad);
-            secp256k1_ge_set_gej(&ggn, &gg);
-            if (i == 15)
-              ad = ggn;
-            for (int k=0; k<sizeof(secp256k1_ge_t); k++)
-              ret->prec[j][k][i] = ((unsigned char*)(&ggn))[k];
+            secp256k1_gej_add(&gj, &gj, &adj);
+            tj[pos++] = gj;
         }
     }
-    secp256k1_ge_set_gej(&ret->fin, &fn);
-    secp256k1_ge_neg(&ret->fin, &ret->fin);
+    assert(pos == 960);
+    tj[pos] = fn;
+    secp256k1_ge_t t[961]; secp256k1_ge_set_all_gej(961, t, tj);
+    pos = 0;
+    const unsigned char *raw = (const unsigned char*)g;
+    for (int j=0; j<64; j++) {
+        for (int k=0; k<sizeof(secp256k1_ge_t); k++)
+            ret->prec[j][k][0] = raw[k];
+        for (int i=1; i<16; i++) {
+            raw = (const unsigned char*)(&t[pos++]);
+            for (int k=0; k<sizeof(secp256k1_ge_t); k++)
+                ret->prec[j][k][i] = raw[k];
+        }
+    }
+    assert(pos == 960);
+    secp256k1_ge_neg(&ret->fin, &t[pos]);
 }
 
 static void secp256k1_ecmult_stop(void) {

src/field.h

@@ -94,6 +94,14 @@ void static secp256k1_fe_inv(secp256k1_fe_t *r, const secp256k1_fe_t *a);
 /** Potentially faster version of secp256k1_fe_inv, without constant-time guarantee. */
 void static secp256k1_fe_inv_var(secp256k1_fe_t *r, const secp256k1_fe_t *a);
 
+/** Calculate the (modular) inverses of a batch of field elements. Requires the inputs' magnitudes to be
+ *  at most 8. The output magnitudes are 1 (but not guaranteed to be normalized). The inputs and
+ *  outputs must not overlap in memory. */
+void static secp256k1_fe_inv_all(size_t len, secp256k1_fe_t r[len], const secp256k1_fe_t a[len]);
+
+/** Potentially faster version of secp256k1_fe_inv_all, without constant-time guarantee. */
+void static secp256k1_fe_inv_all_var(size_t len, secp256k1_fe_t r[len], const secp256k1_fe_t a[len]);
+
 
 /** Convert a field element to a hexadecimal string. */
 void static secp256k1_fe_get_hex(char *r, int *rlen, const secp256k1_fe_t *a);

src/field_impl.h

@@ -214,6 +214,54 @@ void static secp256k1_fe_inv_var(secp256k1_fe_t *r, const secp256k1_fe_t *a) {
 #endif
 }
 
+void static secp256k1_fe_inv_all(size_t len, secp256k1_fe_t r[len], const secp256k1_fe_t a[len]) {
+    if (len < 1)
+        return;
+
+    assert((r + len <= a) || (a + len <= r));
+
+    r[0] = a[0];
+
+    int i = 0;
+    while (++i < len) {
+        secp256k1_fe_mul(&r[i], &r[i - 1], &a[i]);
+    }
+
+    secp256k1_fe_t u; secp256k1_fe_inv(&u, &r[--i]);
+
+    while (i > 0) {
+        int j = i--;
+        secp256k1_fe_mul(&r[j], &r[i], &u);
+        secp256k1_fe_mul(&u, &u, &a[j]);
+    }
+
+    r[0] = u;
+}
+
+void static secp256k1_fe_inv_all_var(size_t len, secp256k1_fe_t r[len], const secp256k1_fe_t a[len]) {
+    if (len < 1)
+        return;
+
+    assert((r + len <= a) || (a + len <= r));
+
+    r[0] = a[0];
+
+    int i = 0;
+    while (++i < len) {
+        secp256k1_fe_mul(&r[i], &r[i - 1], &a[i]);
+    }
+
+    secp256k1_fe_t u; secp256k1_fe_inv_var(&u, &r[--i]);
+
+    while (i > 0) {
+        int j = i--;
+        secp256k1_fe_mul(&r[j], &r[i], &u);
+        secp256k1_fe_mul(&u, &u, &a[j]);
+    }
+
+    r[0] = u;
+}
+
 void static secp256k1_fe_start(void) {
     static const unsigned char secp256k1_fe_consts_p[] = {
         0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,

src/group.h

@@ -68,6 +68,9 @@ void static secp256k1_ge_get_hex(char *r, int *rlen, const secp256k1_ge_t *a);
 /** Set a group element equal to another which is given in jacobian coordinates */
 void static 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 */
+void static secp256k1_ge_set_all_gej(size_t len, secp256k1_ge_t r[len], const secp256k1_gej_t a[len]);
+
 
 /** Set a group element (jacobian) equal to the point at infinity. */
 void static secp256k1_gej_set_infinity(secp256k1_gej_t *r);

src/group_impl.h

@@ -66,6 +66,31 @@ void static secp256k1_ge_set_gej(secp256k1_ge_t *r, secp256k1_gej_t *a) {
     r->y = a->y;
 }
 
+void static secp256k1_ge_set_all_gej(size_t len, secp256k1_ge_t r[len], const secp256k1_gej_t a[len]) {
+    int count = 0;
+    secp256k1_fe_t az[len];
+    for (int i=0; i<len; i++) {
+        if (!a[i].infinity) {
+            az[count++] = a[i].z;
+        }
+    }
+
+    secp256k1_fe_t azi[count];
+    secp256k1_fe_inv_all_var(count, azi, az);
+
+    count = 0;
+    for (int i=0; i<len; i++) {
+        r[i].infinity = a[i].infinity;
+        if (!a[i].infinity) {
+            secp256k1_fe_t *zi = &azi[count++];
+            secp256k1_fe_t zi2; secp256k1_fe_sqr(&zi2, zi);
+            secp256k1_fe_t zi3; 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);
+        }
+    }
+}
+
 void static secp256k1_gej_set_infinity(secp256k1_gej_t *r) {
     r->infinity = 1;
 }

src/tests.c

@@ -244,6 +244,74 @@ void random_fe_non_square(secp256k1_fe_t *ns) {
     }
 }
 
+int check_fe_equal(const secp256k1_fe_t *a, const secp256k1_fe_t *b) {
+    secp256k1_fe_t an = *a; secp256k1_fe_normalize(&an);
+    secp256k1_fe_t bn = *b; secp256k1_fe_normalize(&bn);
+    return secp256k1_fe_equal(&an, &bn);
+}
+
+int check_fe_inverse(const secp256k1_fe_t *a, const secp256k1_fe_t *ai) {
+    secp256k1_fe_t x; secp256k1_fe_mul(&x, a, ai);
+    secp256k1_fe_t one; secp256k1_fe_set_int(&one, 1);
+    return check_fe_equal(&x, &one);
+}
+
+void run_field_inv() {
+    secp256k1_fe_t x, xi, xii;
+    for (int i=0; i<10*count; i++) {
+        random_fe_non_zero(&x);
+        secp256k1_fe_inv(&xi, &x);
+        CHECK(check_fe_inverse(&x, &xi));
+        secp256k1_fe_inv(&xii, &xi);
+        CHECK(check_fe_equal(&x, &xii));
+    }
+}
+
+void run_field_inv_var() {
+    secp256k1_fe_t x, xi, xii;
+    for (int i=0; i<10*count; i++) {
+        random_fe_non_zero(&x);
+        secp256k1_fe_inv_var(&xi, &x);
+        CHECK(check_fe_inverse(&x, &xi));
+        secp256k1_fe_inv_var(&xii, &xi);
+        CHECK(check_fe_equal(&x, &xii));
+    }
+}
+
+void run_field_inv_all() {
+    secp256k1_fe_t x[16], xi[16], xii[16];
+    // Check it's safe to call for 0 elements
+    secp256k1_fe_inv_all(0, xi, x);
+    for (int i=0; i<count; i++) {
+        size_t len = (secp256k1_rand32() & 15) + 1;
+        for (int j=0; j<len; j++)
+            random_fe_non_zero(&x[j]);
+        secp256k1_fe_inv_all(len, xi, x);
+        for (int j=0; j<len; j++)
+            CHECK(check_fe_inverse(&x[j], &xi[j]));
+        secp256k1_fe_inv_all(len, xii, xi);
+        for (int j=0; j<len; j++)
+            CHECK(check_fe_equal(&x[j], &xii[j]));
+    }
+}
+
+void run_field_inv_all_var() {
+    secp256k1_fe_t x[16], xi[16], xii[16];
+    // Check it's safe to call for 0 elements
+    secp256k1_fe_inv_all_var(0, xi, x);
+    for (int i=0; i<count; i++) {
+        size_t len = (secp256k1_rand32() & 15) + 1;
+        for (int j=0; j<len; j++)
+            random_fe_non_zero(&x[j]);
+        secp256k1_fe_inv_all_var(len, xi, x);
+        for (int j=0; j<len; j++)
+            CHECK(check_fe_inverse(&x[j], &xi[j]));
+        secp256k1_fe_inv_all_var(len, xii, xi);
+        for (int j=0; j<len; j++)
+            CHECK(check_fe_equal(&x[j], &xii[j]));
+    }
+}
+
 void test_sqrt(const secp256k1_fe_t *a, const secp256k1_fe_t *k) {
     secp256k1_fe_t r1, r2;
     int v = secp256k1_fe_sqrt(&r1, a);
@@ -537,6 +605,10 @@ int main(int argc, char **argv) {
     run_num_smalltests();
 
     // field tests
+    run_field_inv();
+    run_field_inv_var();
+    run_field_inv_all();
+    run_field_inv_all_var();
     run_sqrt();
 
     // ecmult tests