Add _fe_half and use in _gej_add_ge and _gej_double

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,17 +206,14 @@ 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
  t = t + q
  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

src/bench_internal.c

@@ -140,6 +140,15 @@ void bench_scalar_inverse_var(void* arg, int iters) {
  CHECK(j <= iters);
 }
 
+void bench_field_half(void* arg, int iters) {
+ int i;
+ bench_inv *data = (bench_inv*)arg;
+
+ for (i = 0; i < iters; i++) {
+ secp256k1_fe_half(&data->fe[0]);
+ }
+}
+
 void bench_field_normalize(void* arg, int iters) {
  int i;
  bench_inv *data = (bench_inv*)arg;
@@ -354,6 +363,7 @@ int main(int argc, char **argv) {
  if (d || have_flag(argc, argv, "scalar") || have_flag(argc, argv, "inverse")) run_benchmark("scalar_inverse", bench_scalar_inverse, bench_setup, NULL, &data, 10, iters);
  if (d || have_flag(argc, argv, "scalar") || have_flag(argc, argv, "inverse")) run_benchmark("scalar_inverse_var", bench_scalar_inverse_var, bench_setup, NULL, &data, 10, iters);
 
+ if (d || have_flag(argc, argv, "field") || have_flag(argc, argv, "half")) run_benchmark("field_half", bench_field_half, bench_setup, NULL, &data, 10, iters*100);
  if (d || have_flag(argc, argv, "field") || have_flag(argc, argv, "normalize")) run_benchmark("field_normalize", bench_field_normalize, bench_setup, NULL, &data, 10, iters*100);
  if (d || have_flag(argc, argv, "field") || have_flag(argc, argv, "normalize")) run_benchmark("field_normalize_weak", bench_field_normalize_weak, bench_setup, NULL, &data, 10, iters*100);
  if (d || have_flag(argc, argv, "field") || have_flag(argc, argv, "sqr")) run_benchmark("field_sqr", bench_field_sqr, bench_setup, NULL, &data, 10, iters*10);

src/field.h

@@ -130,4 +130,13 @@ static void secp256k1_fe_storage_cmov(secp256k1_fe_storage *r, const secp256k1_f
 /** If flag is true, set *r equal to *a; otherwise leave it. Constant-time. Both *r and *a must be initialized.*/
 static void secp256k1_fe_cmov(secp256k1_fe *r, const secp256k1_fe *a, int flag);
 
+/** Halves the value of a field element modulo the field prime. Constant-time.
+ * For an input magnitude 'm', the output magnitude is set to 'floor(m/2) + 1'.
+ * The output is not guaranteed to be normalized, regardless of the input. */
+static void secp256k1_fe_half(secp256k1_fe *r);
+
+/** Sets each limb of 'r' to its upper bound at magnitude 'm'. The output will also have its
+ * magnitude set to 'm' and is normalized if (and only if) 'm' is zero. */
+static void secp256k1_fe_get_bounds(secp256k1_fe *r, int m);
+
 #endif /* SECP256K1_FIELD_H */

src/field_10x26_impl.h

@@ -49,6 +49,26 @@ static void secp256k1_fe_verify(const secp256k1_fe *a) {
 }
 #endif
 
+static void secp256k1_fe_get_bounds(secp256k1_fe *r, int m) {
+ VERIFY_CHECK(m >= 0);
+ VERIFY_CHECK(m <= 2048);
+ r->n[0] = 0x3FFFFFFUL * 2 * m;
+ r->n[1] = 0x3FFFFFFUL * 2 * m;
+ r->n[2] = 0x3FFFFFFUL * 2 * m;
+ r->n[3] = 0x3FFFFFFUL * 2 * m;
+ r->n[4] = 0x3FFFFFFUL * 2 * m;
+ r->n[5] = 0x3FFFFFFUL * 2 * m;
+ r->n[6] = 0x3FFFFFFUL * 2 * m;
+ r->n[7] = 0x3FFFFFFUL * 2 * m;
+ r->n[8] = 0x3FFFFFFUL * 2 * m;
+ r->n[9] = 0x03FFFFFUL * 2 * m;
+#ifdef VERIFY
+ r->magnitude = m;
+ r->normalized = (m == 0);
+ secp256k1_fe_verify(r);
+#endif
+}
+
 static void secp256k1_fe_normalize(secp256k1_fe *r) {
  uint32_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4],
  t5 = r->n[5], t6 = r->n[6], t7 = r->n[7], t8 = r->n[8], t9 = r->n[9];
@@ -1133,6 +1153,82 @@ static SECP256K1_INLINE void secp256k1_fe_cmov(secp256k1_fe *r, const secp256k1_
 #endif
 }
 
+static SECP256K1_INLINE void secp256k1_fe_half(secp256k1_fe *r) {
+ uint32_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4],
+ t5 = r->n[5], t6 = r->n[6], t7 = r->n[7], t8 = r->n[8], t9 = r->n[9];
+ uint32_t one = (uint32_t)1;
+ uint32_t mask = -(t0 & one) >> 6;
+
+#ifdef VERIFY
+ secp256k1_fe_verify(r);
+ VERIFY_CHECK(r->magnitude < 32);
+#endif
+
+ /* Bounds analysis (over the rationals).
+ *
+ * Let m = r->magnitude
+ * C = 0x3FFFFFFUL * 2
+ * D = 0x03FFFFFUL * 2
+ *
+ * Initial bounds: t0..t8 <= C * m
+ * t9 <= D * m
+ */
+
+ t0 += 0x3FFFC2FUL & mask;
+ t1 += 0x3FFFFBFUL & mask;
+ t2 += mask;
+ t3 += mask;
+ t4 += mask;
+ t5 += mask;
+ t6 += mask;
+ t7 += mask;
+ t8 += mask;
+ t9 += mask >> 4;
+
+ VERIFY_CHECK((t0 & one) == 0);
+
+ /* t0..t8: added <= C/2
+ * t9: added <= D/2
+ *
+ * Current bounds: t0..t8 <= C * (m + 1/2)
+ * t9 <= D * (m + 1/2)
+ */
+
+ r->n[0] = (t0 >> 1) + ((t1 & one) << 25);
+ r->n[1] = (t1 >> 1) + ((t2 & one) << 25);
+ r->n[2] = (t2 >> 1) + ((t3 & one) << 25);
+ r->n[3] = (t3 >> 1) + ((t4 & one) << 25);
+ r->n[4] = (t4 >> 1) + ((t5 & one) << 25);
+ r->n[5] = (t5 >> 1) + ((t6 & one) << 25);
+ r->n[6] = (t6 >> 1) + ((t7 & one) << 25);
+ r->n[7] = (t7 >> 1) + ((t8 & one) << 25);
+ r->n[8] = (t8 >> 1) + ((t9 & one) << 25);
+ r->n[9] = (t9 >> 1);
+
+ /* t0..t8: shifted right and added <= C/4 + 1/2
+ * t9: shifted right
+ *
+ * Current bounds: t0..t8 <= C * (m/2 + 1/2)
+ * t9 <= D * (m/2 + 1/4)
+ */
+
+#ifdef VERIFY
+ /* Therefore the output magnitude (M) has to be set such that:
+ * t0..t8: C * M >= C * (m/2 + 1/2)
+ * t9: D * M >= D * (m/2 + 1/4)
+ *
+ * It suffices for all limbs that, for any input magnitude m:
+ * M >= m/2 + 1/2
+ *
+ * and since we want the smallest such integer value for M:
+ * M == floor(m/2) + 1
+ */
+ r->magnitude = (r->magnitude >> 1) + 1;
+ r->normalized = 0;
+ secp256k1_fe_verify(r);
+#endif
+}
+
 static SECP256K1_INLINE void secp256k1_fe_storage_cmov(secp256k1_fe_storage *r, const secp256k1_fe_storage *a, int flag) {
  uint32_t mask0, mask1;
  VG_CHECK_VERIFY(r->n, sizeof(r->n));

src/field_5x52_impl.h

@@ -58,6 +58,21 @@ static void secp256k1_fe_verify(const secp256k1_fe *a) {
 }
 #endif
 
+static void secp256k1_fe_get_bounds(secp256k1_fe *r, int m) {
+ VERIFY_CHECK(m >= 0);
+ VERIFY_CHECK(m <= 2048);
+ r->n[0] = 0xFFFFFFFFFFFFFULL * 2 * m;
+ r->n[1] = 0xFFFFFFFFFFFFFULL * 2 * m;
+ r->n[2] = 0xFFFFFFFFFFFFFULL * 2 * m;
+ r->n[3] = 0xFFFFFFFFFFFFFULL * 2 * m;
+ r->n[4] = 0x0FFFFFFFFFFFFULL * 2 * m;
+#ifdef VERIFY
+ r->magnitude = m;
+ r->normalized = (m == 0);
+ secp256k1_fe_verify(r);
+#endif
+}
+
 static void secp256k1_fe_normalize(secp256k1_fe *r) {
  uint64_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4];
 
@@ -477,6 +492,71 @@ static SECP256K1_INLINE void secp256k1_fe_cmov(secp256k1_fe *r, const secp256k1_
 #endif
 }
 
+static SECP256K1_INLINE void secp256k1_fe_half(secp256k1_fe *r) {
+ uint64_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4];
+ uint64_t one = (uint64_t)1;
+ uint64_t mask = -(t0 & one) >> 12;
+
+#ifdef VERIFY
+ secp256k1_fe_verify(r);
+ VERIFY_CHECK(r->magnitude < 32);
+#endif
+
+ /* Bounds analysis (over the rationals).
+ *
+ * Let m = r->magnitude
+ * C = 0xFFFFFFFFFFFFFULL * 2
+ * D = 0x0FFFFFFFFFFFFULL * 2
+ *
+ * Initial bounds: t0..t3 <= C * m
+ * t4 <= D * m
+ */
+
+ t0 += 0xFFFFEFFFFFC2FULL & mask;
+ t1 += mask;
+ t2 += mask;
+ t3 += mask;
+ t4 += mask >> 4;
+
+ VERIFY_CHECK((t0 & one) == 0);
+
+ /* t0..t3: added <= C/2
+ * t4: added <= D/2
+ *
+ * Current bounds: t0..t3 <= C * (m + 1/2)
+ * t4 <= D * (m + 1/2)
+ */
+
+ r->n[0] = (t0 >> 1) + ((t1 & one) << 51);
+ r->n[1] = (t1 >> 1) + ((t2 & one) << 51);
+ r->n[2] = (t2 >> 1) + ((t3 & one) << 51);
+ r->n[3] = (t3 >> 1) + ((t4 & one) << 51);
+ r->n[4] = (t4 >> 1);
+
+ /* t0..t3: shifted right and added <= C/4 + 1/2
+ * t4: shifted right
+ *
+ * Current bounds: t0..t3 <= C * (m/2 + 1/2)
+ * t4 <= D * (m/2 + 1/4)
+ */
+
+#ifdef VERIFY
+ /* Therefore the output magnitude (M) has to be set such that:
+ * t0..t3: C * M >= C * (m/2 + 1/2)
+ * t4: D * M >= D * (m/2 + 1/4)
+ *
+ * It suffices for all limbs that, for any input magnitude m:
+ * M >= m/2 + 1/2
+ *
+ * and since we want the smallest such integer value for M:
+ * M == floor(m/2) + 1
+ */
+ r->magnitude = (r->magnitude >> 1) + 1;
+ r->normalized = 0;
+ secp256k1_fe_verify(r);
+#endif
+}
+
 static SECP256K1_INLINE void secp256k1_fe_storage_cmov(secp256k1_fe_storage *r, const secp256k1_fe_storage *a, int flag) {
  uint64_t mask0, mask1;
  VG_CHECK_VERIFY(r->n, sizeof(r->n));