Split off Z factor addition and change explanation

src/ecmult_gen_impl.h

@@ -50,7 +50,7 @@ static void secp256k1_ecmult_gen_start(void) {
         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, NULL, g);
+        secp256k1_gej_add_ge_var(&nums_gej, &nums_gej, g);
     }
 
     /* compute prec. */

src/ecmult_impl.h

@@ -224,7 +224,16 @@ static void secp256k1_ecmult(secp256k1_gej_t *r, const secp256k1_gej_t *a, const
     int bits = bits_na;
 #endif
 
-    /* calculate odd multiples of 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_fe_t Z;
     secp256k1_ge_t pre_a[ECMULT_TABLE_SIZE(WINDOW_A)];
     secp256k1_ecmult_table_precomp_coz_var(pre_a, &Z, a, WINDOW_A);
@@ -259,28 +268,28 @@ static void secp256k1_ecmult(secp256k1_gej_t *r, const secp256k1_gej_t *a, const
 #ifdef USE_ENDOMORPHISM
         if (i < bits_na_1 && (n = wnaf_na_1[i])) {
             ECMULT_TABLE_GET_GE(&tmpa, pre_a, n, WINDOW_A);
-            secp256k1_gej_add_ge_var(r, r, NULL, &tmpa);
+            secp256k1_gej_add_ge_var(r, r, &tmpa);
         }
         if (i < bits_na_lam && (n = wnaf_na_lam[i])) {
             ECMULT_TABLE_GET_GE(&tmpa, pre_a_lam, n, WINDOW_A);
-            secp256k1_gej_add_ge_var(r, r, NULL, &tmpa);
+            secp256k1_gej_add_ge_var(r, r, &tmpa);
         }
         if (i < bits_ng_1 && (n = wnaf_ng_1[i])) {
             ECMULT_TABLE_GET_GE(&tmpa, c->pre_g, n, WINDOW_G);
-            secp256k1_gej_add_ge_var(r, r, &Z, &tmpa);
+            secp256k1_gej_add_zinv_var(r, r, &tmpa, &Z);
         }
         if (i < bits_ng_128 && (n = wnaf_ng_128[i])) {
             ECMULT_TABLE_GET_GE(&tmpa, c->pre_g_128, n, WINDOW_G);
-            secp256k1_gej_add_ge_var(r, r, &Z, &tmpa);
+            secp256k1_gej_add_zinv_var(r, r, &tmpa, &Z);
         }
 #else
         if (i < bits_na && (n = wnaf_na[i])) {
             ECMULT_TABLE_GET_GE(&tmpa, pre_a, n, WINDOW_A);
-            secp256k1_gej_add_ge_var(r, r, NULL, &tmpa);
+            secp256k1_gej_add_ge_var(r, r, &tmpa);
         }
         if (i < bits_ng && (n = wnaf_ng[i])) {
             ECMULT_TABLE_GET_GE(&tmpa, c->pre_g, n, WINDOW_G);
-            secp256k1_gej_add_ge_var(r, r, &Z, &tmpa);
+            secp256k1_gej_add_zinv_var(r, r, &tmpa, &Z);
         }
 #endif
     }

src/group.h

@@ -116,7 +116,10 @@ static void secp256k1_gej_add_ge(secp256k1_gej_t *r, const secp256k1_gej_t *a, c
 /** 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_fe_t *azr, 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);
+
+/** 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);
 
 /** Get a hex representation of a point. *rlen will be overwritten with the real length. */
 static void secp256k1_gej_get_hex(char *r, int *rlen, const secp256k1_gej_t *a);

src/group_impl.h

@@ -304,33 +304,72 @@ 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_fe_t *azr, 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) {
     if (a->infinity) {
         r->infinity = b->infinity;
-        if (azr == NULL) {
-            r->x = b->x;
-            r->y = b->y;
+        r->x = b->x;
+        r->y = b->y;
+        secp256k1_fe_set_int(&r->z, 1);
+        return;
+    }
+    if (b->infinity) {
+        *r = *a;
+        return;
+    }
+    r->infinity = 0;
+
+    secp256k1_fe_t z12; secp256k1_fe_sqr(&z12, &a->z);
+    secp256k1_fe_t u1 = a->x; secp256k1_fe_normalize_weak(&u1);
+    secp256k1_fe_t u2; secp256k1_fe_mul(&u2, &b->x, &z12);
+    secp256k1_fe_t s1 = a->y; secp256k1_fe_normalize_weak(&s1);
+    secp256k1_fe_t s2; secp256k1_fe_mul(&s2, &b->y, &z12); secp256k1_fe_mul(&s2, &s2, &a->z);
+    secp256k1_fe_t h; secp256k1_fe_negate(&h, &u1, 1); secp256k1_fe_add(&h, &u2);
+    secp256k1_fe_t i; 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);
         } else {
-            secp256k1_fe_t azr2; secp256k1_fe_sqr(&azr2, azr);
-            secp256k1_fe_t azr3; secp256k1_fe_mul(&azr3, &azr2, azr);
-            secp256k1_fe_mul(&r->x, &b->x, &azr2);
-            secp256k1_fe_mul(&r->y, &b->y, &azr3);
+            r->infinity = 1;
         }
-        secp256k1_fe_set_int(&r->z, 1);
         return;
     }
+    secp256k1_fe_t i2; secp256k1_fe_sqr(&i2, &i);
+    secp256k1_fe_t h2; secp256k1_fe_sqr(&h2, &h);
+    secp256k1_fe_t h3; secp256k1_fe_mul(&h3, &h, &h2);
+    r->z = a->z; secp256k1_fe_mul(&r->z, &r->z, &h);
+    secp256k1_fe_t t; 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) {
     if (b->infinity) {
         *r = *a;
         return;
     }
+    if (a->infinity) {
+        r->infinity = b->infinity;
+        secp256k1_fe_t bzinv2; secp256k1_fe_sqr(&bzinv2, bzinv);
+        secp256k1_fe_t bzinv3; 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_t az;
-    if (azr == NULL) {
-        az = a->z;
-    } else {
-        secp256k1_fe_mul(&az, &a->z, azr);
-    }
+    secp256k1_fe_mul(&az, &a->z, bzinv);
 
     secp256k1_fe_t z12; secp256k1_fe_sqr(&z12, &az);
     secp256k1_fe_t u1 = a->x; secp256k1_fe_normalize_weak(&u1);
@@ -358,6 +397,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) {
     VERIFY_CHECK(!b->infinity);
     VERIFY_CHECK(a->infinity == 0 || a->infinity == 1);

src/tests.c

@@ -872,7 +872,7 @@ void test_ge(void) {
             secp256k1_ge_set_gej_var(&ref, &refj);
 
             /* Test gej + ge (var). */
-            secp256k1_gej_add_ge_var(&resj, &gej[i1], NULL, &ge[i2]);
+            secp256k1_gej_add_ge_var(&resj, &gej[i1], &ge[i2]);
             ge_equals_gej(&ref, &resj);
 
             /* Test gej + ge (var, with additional Z factor). */
@@ -882,7 +882,7 @@ void test_ge(void) {
                 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_ge_var(&resj, &gej[i1], &zf, &ge2_zfi); /* The addition undoes the scaling by passing zf. */
+                secp256k1_gej_add_zinv_var(&resj, &gej[i1], &ge2_zfi, &zf);
                 ge_equals_gej(&ref, &resj);
             }