Avoid division when decomposing scalars

- In secp256k1_gej_split_exp, there are two divisions used. Since the denominator is a constant known at compile-time, each can be replaced by a multiplication followed by a right-shift (and rounding).
- Add the constants g1, g2 for this purpose and rewrite secp256k1_scalar_split_lambda_var accordingly.
- Remove secp256k1_num_div since no longer used

Rebased-by: Pieter Wuille

src/ecmult_impl.h

@@ -170,6 +170,8 @@ static void secp256k1_ecmult(secp256k1_gej_t *r, const secp256k1_gej_t *a, const
     /* build wnaf representation for na_1 and na_lam. */
     int wnaf_na_1[129];   int bits_na_1   = secp256k1_ecmult_wnaf(wnaf_na_1,   &na_1,   WINDOW_A);
     int wnaf_na_lam[129]; int bits_na_lam = secp256k1_ecmult_wnaf(wnaf_na_lam, &na_lam, WINDOW_A);
+    VERIFY_CHECK(bits_na_1 <= 129);
+    VERIFY_CHECK(bits_na_lam <= 129);
     int bits = bits_na_1;
     if (bits_na_lam > bits) bits = bits_na_lam;
 #else

src/num.h

@@ -47,9 +47,6 @@ static void secp256k1_num_sub(secp256k1_num_t *r, const secp256k1_num_t *a, cons
 /** Multiply two (signed) numbers. */
 static void secp256k1_num_mul(secp256k1_num_t *r, const secp256k1_num_t *a, const secp256k1_num_t *b);
 
-/** Divide two (signed) numbers. */
-static void secp256k1_num_div(secp256k1_num_t *r, const secp256k1_num_t *a, const secp256k1_num_t *b);
-
 /** Replace a number by its remainder modulo m. M's sign is ignored. The result is a number between 0 and m-1,
     even if r was negative. */
 static void secp256k1_num_mod(secp256k1_num_t *r, const secp256k1_num_t *m);

src/num_gmp_impl.h

@@ -206,25 +206,6 @@ static void secp256k1_num_mul(secp256k1_num_t *r, const secp256k1_num_t *a, cons
     memset(tmp, 0, sizeof(tmp));
 }
 
-static void secp256k1_num_div(secp256k1_num_t *r, const secp256k1_num_t *a, const secp256k1_num_t *b) {
-    secp256k1_num_sanity(a);
-    secp256k1_num_sanity(b);
-    if (b->limbs > a->limbs) {
-        r->limbs = 1;
-        r->data[0] = 0;
-        r->neg = 0;
-        return;
-    }
-
-    mp_limb_t quo[2*NUM_LIMBS+1];
-    mp_limb_t rem[2*NUM_LIMBS+1];
-    mpn_tdiv_qr(quo, rem, 0, a->data, a->limbs, b->data, b->limbs);
-    mpn_copyi(r->data, quo, a->limbs - b->limbs + 1);
-    r->limbs = a->limbs - b->limbs + 1;
-    while (r->limbs > 1 && r->data[r->limbs - 1]==0) r->limbs--;
-    r->neg = a->neg ^ b->neg;
-}
-
 static void secp256k1_num_shift(secp256k1_num_t *r, int bits) {
     if (bits % GMP_NUMB_BITS) {
         // Shift within limbs.

src/scalar_impl.h

@@ -29,7 +29,7 @@ typedef struct {
     secp256k1_num_t order;
 #endif
 #ifdef USE_ENDOMORPHISM
-    secp256k1_num_t a1b2, b1, a2;
+    secp256k1_num_t a1b2, b1, a2, g1, g2;
 #endif
 } secp256k1_scalar_consts_t;
 
@@ -65,10 +65,38 @@ static void secp256k1_scalar_start(void) {
         0x14,0xca,0x50,0xf7,0xa8,0xe2,0xf3,0xf6,
         0x57,0xc1,0x10,0x8d,0x9d,0x44,0xcf,0xd8
     };
-
+    /**
+     * g1, g2 are precomputed constants used to replace division with a rounded multiplication
+     * when decomposing the scalar for an endomorphism-based point multiplication.
+     *
+     * The possibility of using precomputed estimates is mentioned in "Guide to Elliptic Curve
+     * Cryptography" (Hankerson, Menezes, Vanstone) in section 3.5.
+     *
+     * The derivation is described in the paper "Efficient Software Implementation of Public-Key
+     * Cryptography on Sensor Networks Using the MSP430X Microcontroller" (Gouvea, Oliveira, Lopez),
+     * Section 4.3 (here we use a somewhat higher-precision estimate):
+     * d = a1*b2 - b1*a2
+     * g1 = round((2^272)*b2/d)
+     * g2 = round((2^272)*b1/d)
+     *
+     * (Note that 'd' is also equal to the curve order here because [a1,b1] and [a2,b2] are found
+     * as outputs of the Extended Euclidean Algorithm on inputs 'order' and 'lambda').
+     */
+    static const unsigned char secp256k1_scalar_consts_g1[] = {
+        0x30,0x86,
+        0xd2,0x21,0xa7,0xd4,0x6b,0xcd,0xe8,0x6c,
+        0x90,0xe4,0x92,0x84,0xeb,0x15,0x3d,0xab
+    };
+    static const unsigned char secp256k1_scalar_consts_g2[] = {
+        0xe4,0x43,
+        0x7e,0xd6,0x01,0x0e,0x88,0x28,0x6f,0x54,
+        0x7f,0xa9,0x0a,0xbf,0xe4,0xc4,0x22,0x12
+    };
     secp256k1_num_set_bin(&ret->a1b2, secp256k1_scalar_consts_a1b2, sizeof(secp256k1_scalar_consts_a1b2));
     secp256k1_num_set_bin(&ret->a2, secp256k1_scalar_consts_a2, sizeof(secp256k1_scalar_consts_a2));
     secp256k1_num_set_bin(&ret->b1, secp256k1_scalar_consts_b1, sizeof(secp256k1_scalar_consts_b1));
+    secp256k1_num_set_bin(&ret->g1, secp256k1_scalar_consts_g1, sizeof(secp256k1_scalar_consts_g1));
+    secp256k1_num_set_bin(&ret->g2, secp256k1_scalar_consts_g2, sizeof(secp256k1_scalar_consts_g2));
 #endif
 
     /* Set the global pointer. */
@@ -273,26 +301,28 @@ static void secp256k1_scalar_split_lambda_var(secp256k1_scalar_t *r1, secp256k1_
     secp256k1_num_t rn1, rn2;
 
     const secp256k1_scalar_consts_t *c = secp256k1_scalar_consts;
-    secp256k1_num_t bnc1, bnc2, bnt1, bnt2, bnn2;
-
-    secp256k1_num_copy(&bnn2, &c->order);
-    secp256k1_num_shift(&bnn2, 1);
-
-    secp256k1_num_mul(&bnc1, &na, &c->a1b2);
-    secp256k1_num_add(&bnc1, &bnc1, &bnn2);
-    secp256k1_num_div(&bnc1, &bnc1, &c->order);
-
-    secp256k1_num_mul(&bnc2, &na, &c->b1);
-    secp256k1_num_add(&bnc2, &bnc2, &bnn2);
-    secp256k1_num_div(&bnc2, &bnc2, &c->order);
-
-    secp256k1_num_mul(&bnt1, &bnc1, &c->a1b2);
-    secp256k1_num_mul(&bnt2, &bnc2, &c->a2);
-    secp256k1_num_add(&bnt1, &bnt1, &bnt2);
-    secp256k1_num_sub(&rn1, &na, &bnt1);
-    secp256k1_num_mul(&bnt1, &bnc1, &c->b1);
-    secp256k1_num_mul(&bnt2, &bnc2, &c->a1b2);
-    secp256k1_num_sub(&rn2, &bnt1, &bnt2);
+    secp256k1_num_t d1, d2, t, one;
+    unsigned char cone[1] = {0x01};
+    secp256k1_num_set_bin(&one, cone, 1);
+
+    secp256k1_num_mul(&d1, &na, &c->g1);
+    secp256k1_num_shift(&d1, 271);
+    secp256k1_num_add(&d1, &d1, &one);
+    secp256k1_num_shift(&d1, 1);
+
+    secp256k1_num_mul(&d2, &na, &c->g2);
+    secp256k1_num_shift(&d2, 271);
+    secp256k1_num_add(&d2, &d2, &one);
+    secp256k1_num_shift(&d2, 1);
+
+    secp256k1_num_mul(&t, &d1, &c->a1b2);
+    secp256k1_num_sub(&rn1, &na, &t);
+    secp256k1_num_mul(&t, &d2, &c->a2);
+    secp256k1_num_sub(&rn1, &rn1, &t);
+
+    secp256k1_num_mul(&rn2, &d1, &c->b1);
+    secp256k1_num_mul(&t, &d2, &c->a1b2);
+    secp256k1_num_sub(&rn2, &rn2, &t);
 
     secp256k1_num_get_bin(b, 32, &rn1);
     secp256k1_scalar_set_b32(r1, b, NULL);