group: Save a normalize_to_zero in gej_add_ge

The code currently switches to the alternative formula for lambda only if (R,M)
= (0,0) but the alternative formula works whenever M = 0: Specifically, M = 0
implies y1 = -y2. If x1 = x2, then a = -b this is the r = infinity case that we
handle separately. If x1 != x2, then the denominator in the alternative formula
is non-zero, so this formula is well-defined.

One needs to carefully check that the infinity assignment is still correct
because now the definition of m_alt at this point in the code has changed. But
this is true:

Case y1 = -y2:
  Then degenerate = true and infinity = ((x1 - x2)Z == 0) & ~a->infinity .
  a->infinity is handled separately.
  And if ~a->infinity, then Z = Z1 != 0,
  so infinity = (x1 - x2 == 0) = (a == -b) by case condition.

Case y1 != -y2:
  Then degenerate = false and infinity = ((y1 + y2)Z == 0) & ~a->infinity .
  a->infinity is handled separately.
  And if ~a->infinity, then Z = Z1 != 0,
  so infinity = (y1 + y2 == 0) = false by case condition.

Co-Authored-By: Pieter Wuille <pieter@wuille.net>

sage/prove_group_implementations.sage

@@ -157,7 +157,7 @@ def formula_secp256k1_gej_add_ge(branch, a, b):
   zeroes = {}
   nonzeroes = {}
   a_infinity = False
-  if (branch & 4) != 0:
+  if (branch & 2) != 0:
     nonzeroes.update({a.Infinity : 'a_infinite'})
     a_infinity = True
   else:
@@ -176,15 +176,11 @@ def formula_secp256k1_gej_add_ge(branch, a, b):
   m_alt = -u2
   tt = u1 * m_alt
   rr = rr + tt
-  degenerate = (branch & 3) == 3
-  if (branch & 1) != 0:
+  degenerate = (branch & 1) != 0
+  if degenerate:
     zeroes.update({m : 'm_zero'})
   else:
     nonzeroes.update({m : 'm_nonzero'})
-  if (branch & 2) != 0:
-    zeroes.update({rr : 'rr_zero'})
-  else:
-    nonzeroes.update({rr : 'rr_nonzero'})
   rr_alt = s1
   rr_alt = rr_alt * 2
   m_alt = m_alt + u1
@@ -200,8 +196,7 @@ def formula_secp256k1_gej_add_ge(branch, a, b):
   t = rr_alt^2
   rz = a.Z * m_alt
   infinity = False
-  if (branch & 8) != 0:
-    if not a_infinity:
+  if (branch & 4) != 0:
     infinity = True
     zeroes.update({rz : 'r.z = 0'})
   else:
@@ -289,14 +284,14 @@ if __name__ == "__main__":
   success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_var", 0, 7, 5, formula_secp256k1_gej_add_var)
   success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge_var", 0, 7, 5, formula_secp256k1_gej_add_ge_var)
   success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_zinv_var", 0, 7, 5, formula_secp256k1_gej_add_zinv_var)
-  success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 16, formula_secp256k1_gej_add_ge)
+  success = success & check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 8, formula_secp256k1_gej_add_ge)
   success = success & (not check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge_old [should fail]", 0, 7, 4, formula_secp256k1_gej_add_ge_old))
 
   if len(sys.argv) >= 2 and sys.argv[1] == "--exhaustive":
     success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_var", 0, 7, 5, formula_secp256k1_gej_add_var, 43)
     success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge_var", 0, 7, 5, formula_secp256k1_gej_add_ge_var, 43)
     success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_zinv_var", 0, 7, 5, formula_secp256k1_gej_add_zinv_var, 43)
-    success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 16, formula_secp256k1_gej_add_ge, 43)
+    success = success & check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 8, formula_secp256k1_gej_add_ge, 43)
     success = success & (not check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge_old [should fail]", 0, 7, 4, formula_secp256k1_gej_add_ge_old, 43))
 
   sys.exit(int(not success))

src/group_impl.h

@@ -558,10 +558,9 @@ static void secp256k1_gej_add_ge(secp256k1_gej *r, const secp256k1_gej *a, const
     secp256k1_fe_negate(&m_alt, &u2, 1);                /* Malt = -X2*Z1^2 */
     secp256k1_fe_mul(&tt, &u1, &m_alt);                 /* tt = -U1*U2 (2) */
     secp256k1_fe_add(&rr, &tt);                         /* rr = R = T^2-U1*U2 (3) */
-    /** If lambda = R/M = 0/0 we have a problem (except in the "trivial"
+    /** If lambda = R/M = R/0 we have a problem (except in the "trivial"
      *  case that Z = z1z2 = 0, and this is special-cased later on). */
-    degenerate = secp256k1_fe_normalizes_to_zero(&m) &
-                 secp256k1_fe_normalizes_to_zero(&rr);
+    degenerate = secp256k1_fe_normalizes_to_zero(&m);
     /* This only occurs when y1 == -y2 and x1^3 == x2^3, but x1 != x2.
      * This means either x1 == beta*x2 or beta*x1 == x2, where beta is
      * a nontrivial cube root of one. In either case, an alternate
@@ -573,7 +572,7 @@ static void secp256k1_gej_add_ge(secp256k1_gej *r, const secp256k1_gej *a, const
 
     secp256k1_fe_cmov(&rr_alt, &rr, !degenerate);
     secp256k1_fe_cmov(&m_alt, &m, !degenerate);
-    /* Now Ralt / Malt = lambda and is guaranteed not to be 0/0.
+    /* Now Ralt / Malt = lambda and is guaranteed not to be Ralt / 0.
      * From here on out Ralt and Malt represent the numerator
      * and denominator of lambda; R and M represent the explicit
      * expressions x1^2 + x2^2 + x1x2 and y1 + y2. */