Effective affine addition in EC multiplication

* Make secp256k1_gej_add_var and secp256k1_gej_double return the
  Z ratio to go from a.z to r.z.
* Use these Z ratios to speed up batch point conversion to affine
  coordinates, and to speed up batch conversion of points to a
  common Z coordinate.
* Add a point addition function that takes a point with a known
  Z inverse.
* Due to secp256k1's endomorphism, all additions in the EC
  multiplication code can work on affine coordinate (with an
  implicit common Z coordinate), correcting the Z coordinate of
  the result afterwards.

Refactoring by Pieter Wuille:
* Move more global-z logic into the group code.
* Separate code for computing the odd multiples from the code to bring it
  to either storage or globalz format.
* Rename functions.
* Make all addition operations return Z ratios, and test them.
* Make the zr table format compatible with future batch chaining
  (the first entry in zr becomes the ratio between the input and the
  first output).

Original idea and code by Peter Dettman.

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,19 +994,62 @@ 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]);
             ge_equals_gej(&ref, &resj);
 
+            /* 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) {
                 /* secp256k1_gej_add_ge does not support its second argument being infinity. */
@@ -1012,10 +1059,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 +1106,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 +1204,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 +1227,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);