Merge bitcoin-core/secp256k1#1033: Add _fe_half and use in _gej_add_ge and _gej_double

e848c3799c Update sage files for new formulae (Peter Dettman)
d64bb5d4f3 Add fe_half tests for worst-case inputs (Peter Dettman)
4eb8b932ff Further improve doubling formula using fe_half (Peter Dettman)
557b31fac3 Doubling formula using fe_half (Pieter Wuille)
2cbb4b1a42 Run more iterations of run_field_misc (Pieter Wuille)
9cc5c257ed Add test for secp256k1_fe_half (Pieter Wuille)
925f78d55e Add _fe_half and use in _gej_add_ge (Peter Dettman)

Pull request description:

  - Trades 1 _half for 3 _mul_int and 2 _normalize_weak

  Gives around 2-3% faster signing and ECDH, depending on compiler/platform.

ACKs for top commit:
  sipa:
    utACK e848c3799c
  jonasnick:
    ACK e848c3799c
  real-or-random:
    ACK e848c3799c

Tree-SHA512: 81a6c93b3d983f1b48ec8e8b6f262ba914215045a95415147f41ee6e85296aa4d0cbbad9f370cdf475571447baad861d2cc8e0b04a71202d48959cb8a098f584

sage/prove_group_implementations.sage

@@ -8,25 +8,20 @@ load("weierstrass_prover.sage")
 def formula_secp256k1_gej_double_var(a):
   """libsecp256k1's secp256k1_gej_double_var, used by various addition functions"""
   rz = a.Z * a.Y
-  rz = rz * 2
-  t1 = a.X^2
-  t1 = t1 * 3
-  t2 = t1^2
-  t3 = a.Y^2
-  t3 = t3 * 2
-  t4 = t3^2
-  t4 = t4 * 2
-  t3 = t3 * a.X
-  rx = t3
-  rx = rx * 4
-  rx = -rx
-  rx = rx + t2
-  t2 = -t2
-  t3 = t3 * 6
-  t3 = t3 + t2
-  ry = t1 * t3
-  t2 = -t4
-  ry = ry + t2
+  s = a.Y^2
+  l = a.X^2
+  l = l * 3
+  l = l / 2
+  t = -s
+  t = t * a.X
+  rx = l^2
+  rx = rx + t
+  rx = rx + t
+  s = s^2
+  t = t + rx
+  ry = t * l
+  ry = ry + s
+  ry = -ry
   return jacobianpoint(rx, ry, rz)
 
 def formula_secp256k1_gej_add_var(branch, a, b):
@@ -197,7 +192,8 @@ def formula_secp256k1_gej_add_ge(branch, a, b):
     rr_alt = rr
     m_alt = m
   n = m_alt^2
-  q = n * t
+  q = -t
+  q = q * n
   n = n^2
   if degenerate:
     n = m
@@ -210,8 +206,6 @@ def formula_secp256k1_gej_add_ge(branch, a, b):
     zeroes.update({rz : 'r.z=0'})
   else:
     nonzeroes.update({rz : 'r.z!=0'})
-  rz = rz * 2
-  q = -q
   t = t + q
   rx = t
   t = t * 2
@@ -219,8 +213,7 @@ def formula_secp256k1_gej_add_ge(branch, a, b):
   t = t * rr_alt
   t = t + n
   ry = -t
-  rx = rx * 4
-  ry = ry * 4
+  ry = ry / 2
   if a_infinity:
     rx = b.X
     ry = b.Y