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_gej_split_exp accordingly.
- Remove secp256k1_num_div since no longer used

Author: peterdettman

src/group.h

@@ -32,7 +32,7 @@ typedef struct {
 #ifdef USE_ENDOMORPHISM
     // constants related to secp256k1's efficiently computable endomorphism
     secp256k1_fe_t beta;
-    secp256k1_num_t lambda, a1b2, b1, a2;
+    secp256k1_num_t lambda, a1b2, b1, a2, g1, g2;
 #endif
 } secp256k1_ge_consts_t;
 

src/group_impl.h

@@ -337,26 +337,27 @@ void static secp256k1_gej_mul_lambda(secp256k1_gej_t *r, const secp256k1_gej_t *
 
 void static secp256k1_gej_split_exp_var(secp256k1_num_t *r1, secp256k1_num_t *r2, const secp256k1_num_t *a) {
     const secp256k1_ge_consts_t *c = secp256k1_ge_consts;
-    secp256k1_num_t bnc1, bnc2, bnt1, bnt2, bnn2;
-
-    secp256k1_num_copy(&bnn2, &c->order);
-    secp256k1_num_shift(&bnn2, 1);
-
-    secp256k1_num_mul(&bnc1, a, &c->a1b2);
-    secp256k1_num_add(&bnc1, &bnc1, &bnn2);
-    secp256k1_num_div(&bnc1, &bnc1, &c->order);
-
-    secp256k1_num_mul(&bnc2, a, &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(r1, a, &bnt1);
-    secp256k1_num_mul(&bnt1, &bnc1, &c->b1);
-    secp256k1_num_mul(&bnt2, &bnc2, &c->a1b2);
-    secp256k1_num_sub(r2, &bnt1, &bnt2);
+
+    secp256k1_num_t d1, d2, t;
+
+    secp256k1_num_mul(&t, a, &c->g1);
+    secp256k1_num_trunc(&d1, &t, 272);
+    if (secp256k1_num_get_bit(&t, 271))
+        secp256k1_num_inc(&d1);
+
+    secp256k1_num_mul(&t, a, &c->g2);
+    secp256k1_num_trunc(&d2, &t, 272);
+    if (secp256k1_num_get_bit(&t, 271))
+        secp256k1_num_inc(&d2);
+
+    secp256k1_num_mul(&t, &d1, &c->a1b2);
+    secp256k1_num_sub(r1, a, &t);
+    secp256k1_num_mul(&t, &d2, &c->a2);
+    secp256k1_num_sub(r1, r1, &t);
+
+    secp256k1_num_mul(r2, &d1, &c->b1);
+    secp256k1_num_mul(&t, &d2, &c->a1b2);
+    secp256k1_num_sub(r2, r2, &t);
 }
 #endif
 
@@ -407,6 +408,33 @@ void static secp256k1_ge_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_ge_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_ge_consts_g2[] = {
+        0xe4,0x43,
+        0x7e,0xd6,0x01,0x0e,0x88,0x28,0x6f,0x54,
+        0x7f,0xa9,0x0a,0xbf,0xe4,0xc4,0x22,0x12
+    };
 #endif
     if (secp256k1_ge_consts == NULL) {
         secp256k1_ge_consts_t *ret = (secp256k1_ge_consts_t*)malloc(sizeof(secp256k1_ge_consts_t));
@@ -418,6 +446,8 @@ void static secp256k1_ge_start(void) {
         secp256k1_num_set_bin(&ret->a1b2,   secp256k1_ge_consts_a1b2,   sizeof(secp256k1_ge_consts_a1b2));
         secp256k1_num_set_bin(&ret->a2,     secp256k1_ge_consts_a2,     sizeof(secp256k1_ge_consts_a2));
         secp256k1_num_set_bin(&ret->b1,     secp256k1_ge_consts_b1,     sizeof(secp256k1_ge_consts_b1));
+        secp256k1_num_set_bin(&ret->g1,     secp256k1_ge_consts_g1,     sizeof(secp256k1_ge_consts_g1));
+        secp256k1_num_set_bin(&ret->g2,     secp256k1_ge_consts_g2,     sizeof(secp256k1_ge_consts_g2));
         secp256k1_fe_set_b32(&ret->beta, secp256k1_ge_consts_beta);
 #endif
         secp256k1_fe_t g_x, g_y;

src/num.h

@@ -52,9 +52,6 @@ void static secp256k1_num_sub(secp256k1_num_t *r, const secp256k1_num_t *a, cons
 /** Multiply two (signed) numbers. */
 void static secp256k1_num_mul(secp256k1_num_t *r, const secp256k1_num_t *a, const secp256k1_num_t *b);
 
-/** Divide two (signed) numbers. */
-void static 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. */
 void static secp256k1_num_mod(secp256k1_num_t *r, const secp256k1_num_t *m);
@@ -89,6 +86,9 @@ void static secp256k1_num_get_hex(char *r, int rlen, const secp256k1_num_t *a);
 /** Split a number into a low and high part. */
 void static secp256k1_num_split(secp256k1_num_t *rl, secp256k1_num_t *rh, const secp256k1_num_t *a, int bits);
 
+/** Truncate a number, returning the high part only. */
+void static secp256k1_num_trunc(secp256k1_num_t *r, const secp256k1_num_t *a, int bits);
+
 /** Change a number's sign. */
 void static secp256k1_num_negate(secp256k1_num_t *r);
 

src/num_gmp_impl.h

@@ -234,25 +234,6 @@ void static secp256k1_num_mul(secp256k1_num_t *r, const secp256k1_num_t *a, cons
     memset(tmp, 0, sizeof(tmp));
 }
 
-void static 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;
-}
-
 void static secp256k1_num_mod_mul(secp256k1_num_t *r, const secp256k1_num_t *a, const secp256k1_num_t *b, const secp256k1_num_t *m) {
     secp256k1_num_mul(r, a, b);
     secp256k1_num_mod(r, m);
@@ -359,6 +340,24 @@ void static secp256k1_num_split(secp256k1_num_t *rl, secp256k1_num_t *rh, const
     while (rl->limbs>1 && rl->data[rl->limbs-1]==0) rl->limbs--;
 }
 
+void static secp256k1_num_trunc(secp256k1_num_t *r, const secp256k1_num_t *a, int bits) {
+    assert(bits > 0);
+    r->neg = a->neg;
+    if (bits >= a->limbs * GMP_NUMB_BITS) {
+        r->limbs = 1;
+        r->data[0] = 0;
+        return;
+    }
+    int limbs = bits / GMP_NUMB_BITS, left = bits % GMP_NUMB_BITS;
+    r->limbs = a->limbs - limbs;
+    if (left == 0) {
+        mpn_copyi(r->data, a->data + limbs, r->limbs);
+    } else {
+        mpn_rshift(r->data, a->data + limbs, r->limbs, left);
+        while (r->limbs>1 && r->data[r->limbs-1]==0) r->limbs--;
+    }
+}
+
 void static secp256k1_num_negate(secp256k1_num_t *r) {
     r->neg ^= 1;
 }