diff --git a/Makefile.am b/Makefile.am index 6c7a9cad99..58c9635e53 100644 --- a/Makefile.am +++ b/Makefile.am @@ -30,6 +30,10 @@ noinst_HEADERS += src/field_5x52.h noinst_HEADERS += src/field_5x52_impl.h noinst_HEADERS += src/field_5x52_int128_impl.h noinst_HEADERS += src/field_5x52_asm_impl.h +noinst_HEADERS += src/modinv32.h +noinst_HEADERS += src/modinv32_impl.h +noinst_HEADERS += src/modinv64.h +noinst_HEADERS += src/modinv64_impl.h noinst_HEADERS += src/assumptions.h noinst_HEADERS += src/util.h noinst_HEADERS += src/scratch.h diff --git a/src/field_10x26_impl.h b/src/field_10x26_impl.h index a2c16a651a..1d2d916d7a 100644 --- a/src/field_10x26_impl.h +++ b/src/field_10x26_impl.h @@ -9,6 +9,7 @@ #include "util.h" #include "field.h" +#include "modinv32_impl.h" #ifdef VERIFY static void secp256k1_fe_verify(const secp256k1_fe *a) { @@ -1164,25 +1165,15 @@ static SECP256K1_INLINE void secp256k1_fe_from_storage(secp256k1_fe *r, const se #endif } -static void secp256k1_fe_decode_30(secp256k1_fe *r, const int32_t *a) { - +static void secp256k1_fe_from_signed30(secp256k1_fe *r, const secp256k1_modinv32_signed30 *a) { const uint32_t M26 = UINT32_MAX >> 6; - const uint32_t a0 = a[0], a1 = a[1], a2 = a[2], a3 = a[3], a4 = a[4], - a5 = a[5], a6 = a[6], a7 = a[7], a8 = a[8]; + const uint32_t a0 = a->v[0], a1 = a->v[1], a2 = a->v[2], a3 = a->v[3], a4 = a->v[4], + a5 = a->v[5], a6 = a->v[6], a7 = a->v[7], a8 = a->v[8]; uint32_t r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, t; t = (int32_t)a8 >> 16; - /* a must be in the range [-2^256, 2^256). */ - VERIFY_CHECK(a0 >> 30 == 0); - VERIFY_CHECK(a1 >> 30 == 0); - VERIFY_CHECK(a2 >> 30 == 0); - VERIFY_CHECK(a3 >> 30 == 0); - VERIFY_CHECK(a4 >> 30 == 0); - VERIFY_CHECK(a5 >> 30 == 0); - VERIFY_CHECK(a6 >> 30 == 0); - VERIFY_CHECK(a7 >> 30 == 0); - VERIFY_CHECK(t == 0 || t == -(uint32_t)1); + secp256k1_modinv32_signed30_verify(a); /* Add 2P if a8 is "negative". */ r0 = 0x3FFF85EUL & t; @@ -1225,8 +1216,7 @@ static void secp256k1_fe_decode_30(secp256k1_fe *r, const int32_t *a) { #endif } -static void secp256k1_fe_encode_30(int32_t *r, const secp256k1_fe *a) { - +static void secp256k1_fe_to_signed30(secp256k1_modinv32_signed30 *r, const secp256k1_fe *a) { const uint32_t M30 = UINT32_MAX >> 2; const uint32_t *n = &a->n[0]; const uint64_t a0 = n[0], a1 = n[1], a2 = n[2], a3 = n[3], a4 = n[4], @@ -1236,399 +1226,72 @@ static void secp256k1_fe_encode_30(int32_t *r, const secp256k1_fe *a) { VERIFY_CHECK(a->normalized); #endif - r[0] = (a0 | a1 << 26) & M30; - r[1] = (a1 >> 4 | a2 << 22) & M30; - r[2] = (a2 >> 8 | a3 << 18) & M30; - r[3] = (a3 >> 12 | a4 << 14) & M30; - r[4] = (a4 >> 16 | a5 << 10) & M30; - r[5] = (a5 >> 20 | a6 << 6) & M30; - r[6] = (a6 >> 24 | a7 << 2 - | a8 << 28) & M30; - r[7] = (a8 >> 2 | a9 << 24) & M30; - r[8] = a9 >> 6; -} - -static uint32_t secp256k1_fe_divsteps_30(uint32_t eta, uint32_t f0, uint32_t g0, int32_t *t) { - - uint32_t u = 1, v = 0, q = 0, r = 1; - uint32_t c1, c2, f = f0, g = g0, x, y, z; - int i; - - for (i = 0; i < 30; ++i) { - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << i); - VERIFY_CHECK((q * f0 + r * g0) == g << i); - - c1 = (int32_t)eta >> 31; - c2 = -(g & 1); - - x = (f ^ c1) - c1; - y = (u ^ c1) - c1; - z = (v ^ c1) - c1; - - g += x & c2; - q += y & c2; - r += z & c2; - - c1 &= c2; - eta = (eta ^ c1) - (c1 + 1); - - f += g & c1; - u += q & c1; - v += r & c1; - - g >>= 1; - u <<= 1; - v <<= 1; - } - - t[0] = (int32_t)u; - t[1] = (int32_t)v; - t[2] = (int32_t)q; - t[3] = (int32_t)r; - - return eta; -} - -static uint32_t secp256k1_fe_divsteps_30_var(uint32_t eta, uint32_t f0, uint32_t g0, int32_t *t) { - static const uint8_t inv256[128] = { - 0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59, - 0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31, - 0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89, - 0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61, - 0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9, - 0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91, - 0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9, - 0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1, - 0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19, - 0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1, - 0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01 - }; - - uint32_t u = 1, v = 0, q = 0, r = 1; - uint32_t f = f0, g = g0, m, w, x, y, z; - int i = 30, limit, zeros; - - for (;;) { - /* Use a sentinel bit to count zeros only up to i. */ - zeros = secp256k1_ctz32_var(g | (UINT32_MAX << i)); - - g >>= zeros; - u <<= zeros; - v <<= zeros; - eta -= zeros; - i -= zeros; - - if (i <= 0) { - break; - } - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((g & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << (30 - i)); - VERIFY_CHECK((q * f0 + r * g0) == g << (30 - i)); - - if ((int32_t)eta < 0) { - eta = -eta; - x = f; f = g; g = -x; - y = u; u = q; q = -y; - z = v; v = r; r = -z; - } - - /* Handle up to 8 divsteps at once, subject to eta and i. */ - limit = ((int)eta + 1) > i ? i : ((int)eta + 1); - m = (UINT32_MAX >> (32 - limit)) & 255U; - - w = (g * inv256[(f >> 1) & 127]) & m; - - g += f * w; - q += u * w; - r += v * w; - - VERIFY_CHECK((g & m) == 0); - } - - t[0] = (int32_t)u; - t[1] = (int32_t)v; - t[2] = (int32_t)q; - t[3] = (int32_t)r; - - return eta; + r->v[0] = (a0 | a1 << 26) & M30; + r->v[1] = (a1 >> 4 | a2 << 22) & M30; + r->v[2] = (a2 >> 8 | a3 << 18) & M30; + r->v[3] = (a3 >> 12 | a4 << 14) & M30; + r->v[4] = (a4 >> 16 | a5 << 10) & M30; + r->v[5] = (a5 >> 20 | a6 << 6) & M30; + r->v[6] = (a6 >> 24 | a7 << 2 + | a8 << 28) & M30; + r->v[7] = (a8 >> 2 | a9 << 24) & M30; + r->v[8] = a9 >> 6; } -static void secp256k1_fe_update_de_30(int32_t *d, int32_t *e, const int32_t *t) { - - /* P == 2^256 - 2^32 - C30 */ - const int64_t C30 = 0x3D1L; - /* I30 == -P^-1 mod 2^30 */ - const uint32_t I30 = 0x12253531L; - const int32_t M30 = (int32_t)(UINT32_MAX >> 2); - const int32_t u = t[0], v = t[1], q = t[2], r = t[3]; - int32_t di, ei, md, me; - int64_t cd, ce; - int i; - - di = d[0]; - ei = e[0]; - - cd = (int64_t)u * di + (int64_t)v * ei; - ce = (int64_t)q * di + (int64_t)r * ei; - - /* Calculate the multiples of P to add, to zero the 30 bottom bits. We choose md, me - * from the centred range [-2^29, 2^29) to keep d, e within [-2^256, 2^256). */ - md = ((int32_t)(I30 * 4 * (uint32_t)cd)) >> 2; - me = ((int32_t)(I30 * 4 * (uint32_t)ce)) >> 2; - - cd -= (int64_t)C30 * md; - ce -= (int64_t)C30 * me; - - VERIFY_CHECK(((int32_t)cd & M30) == 0); cd >>= 30; - VERIFY_CHECK(((int32_t)ce & M30) == 0); ce >>= 30; - - cd -= (int64_t)4 * md; - ce -= (int64_t)4 * me; - - for (i = 1; i < 8; ++i) { - - di = d[i]; - ei = e[i]; - - cd += (int64_t)u * di + (int64_t)v * ei; - ce += (int64_t)q * di + (int64_t)r * ei; - - d[i - 1] = (int32_t)cd & M30; cd >>= 30; - e[i - 1] = (int32_t)ce & M30; ce >>= 30; - } - - { - di = d[8]; - ei = e[8]; - - cd += (int64_t)u * di + (int64_t)v * ei; - ce += (int64_t)q * di + (int64_t)r * ei; - - cd += (int64_t)65536 * md; - ce += (int64_t)65536 * me; - - d[7] = (int32_t)cd & M30; cd >>= 30; - e[7] = (int32_t)ce & M30; ce >>= 30; - } - - d[8] = (int32_t)cd; - e[8] = (int32_t)ce; -} - -static void secp256k1_fe_update_fg_30(int32_t *f, int32_t *g, const int32_t *t) { - - const int32_t M30 = (int32_t)(UINT32_MAX >> 2); - const int32_t u = t[0], v = t[1], q = t[2], r = t[3]; - int32_t fi, gi; - int64_t cf, cg; - int i; - - fi = f[0]; - gi = g[0]; +static const secp256k1_modinv32_modinfo secp256k1_const_modinfo_fe = { + {{-0x3D1, -4, 0, 0, 0, 0, 0, 0, 65536}}, + 0x12253531L * 4 +}; - cf = (int64_t)u * fi + (int64_t)v * gi; - cg = (int64_t)q * fi + (int64_t)r * gi; - - VERIFY_CHECK(((int32_t)cf & M30) == 0); cf >>= 30; - VERIFY_CHECK(((int32_t)cg & M30) == 0); cg >>= 30; - - for (i = 1; i < 9; ++i) { - - fi = f[i]; - gi = g[i]; - - cf += (int64_t)u * fi + (int64_t)v * gi; - cg += (int64_t)q * fi + (int64_t)r * gi; - - f[i - 1] = (int32_t)cf & M30; cf >>= 30; - g[i - 1] = (int32_t)cg & M30; cg >>= 30; - } - - f[8] = (int32_t)cf; - g[8] = (int32_t)cg; -} - -static void secp256k1_fe_update_fg_30_var(int len, int32_t *f, int32_t *g, const int32_t *t) { - - const int32_t M30 = (int32_t)(UINT32_MAX >> 2); - const int32_t u = t[0], v = t[1], q = t[2], r = t[3]; - int32_t fi, gi; - int64_t cf, cg; - int i; - - VERIFY_CHECK(len > 0); - - fi = f[0]; - gi = g[0]; - - cf = (int64_t)u * fi + (int64_t)v * gi; - cg = (int64_t)q * fi + (int64_t)r * gi; - - VERIFY_CHECK(((int32_t)cf & M30) == 0); cf >>= 30; - VERIFY_CHECK(((int32_t)cg & M30) == 0); cg >>= 30; - - for (i = 1; i < len; ++i) { - - fi = f[i]; - gi = g[i]; - - cf += (int64_t)u * fi + (int64_t)v * gi; - cg += (int64_t)q * fi + (int64_t)r * gi; - - f[i - 1] = (int32_t)cf & M30; cf >>= 30; - g[i - 1] = (int32_t)cg & M30; cg >>= 30; - } - - f[len - 1] = (int32_t)cf; - g[len - 1] = (int32_t)cg; -} - -static void secp256k1_fe_inv(secp256k1_fe *r, const secp256k1_fe *a) { - - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int32_t t[4]; - int32_t d[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t e[9] = { 1, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t f[9] = { 0x3FFFFC2F, 0x3FFFFFFB, 0x3FFFFFFF, 0x3FFFFFFF, 0x3FFFFFFF, - 0x3FFFFFFF, 0x3FFFFFFF, 0x3FFFFFFF, 0xFFFF }; - int32_t g[9]; - secp256k1_fe b0, b1; - int i, sign; - uint32_t eta; +static void secp256k1_fe_inv(secp256k1_fe *r, const secp256k1_fe *x) { + int sign; + secp256k1_fe tmp; + secp256k1_modinv32_signed30 s; #ifdef VERIFY int zero_in; #endif - b0 = *a; - secp256k1_fe_normalize(&b0); - secp256k1_fe_encode_30(g, &b0); - + tmp = *x; + secp256k1_fe_normalize(&tmp); #ifdef VERIFY - zero_in = secp256k1_fe_is_zero(&b0); + zero_in = secp256k1_fe_normalizes_to_zero(&tmp); #endif - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If the maximum bitlength of g is known to be less than 256, then eta can be set - * initially to -(1 + (256 - maxlen(g))), and only (741 - (256 - maxlen(g))) total - * divsteps are needed. */ - eta = -(uint32_t)1; - - for (i = 0; i < 25; ++i) { - eta = secp256k1_fe_divsteps_30(eta, f[0], g[0], t); - secp256k1_fe_update_de_30(d, e, t); - secp256k1_fe_update_fg_30(f, g, t); - } - - /* At this point sufficient iterations have been performed that g must have reached 0 - * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g - * values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - VERIFY_CHECK((g[0] | g[1] | g[2] | g[3] | g[4] | g[5] | g[6] | g[7] | g[8]) == 0); - - sign = (f[0] >> 1) & 1; - - secp256k1_fe_decode_30(&b0, d); - - secp256k1_fe_negate(&b1, &b0, 2); - secp256k1_fe_cmov(&b0, &b1, sign); - secp256k1_fe_normalize_weak(&b0); + secp256k1_fe_to_signed30(&s, &tmp); + secp256k1_modinv32(&sign, &s, &secp256k1_const_modinfo_fe); + secp256k1_fe_from_signed30(r, &s); + secp256k1_fe_negate(&tmp, r, 1); + secp256k1_fe_cmov(r, &tmp, sign); + secp256k1_fe_normalize_weak(r); #ifdef VERIFY - VERIFY_CHECK(!secp256k1_fe_normalizes_to_zero(&b0) == !zero_in); + VERIFY_CHECK(secp256k1_fe_normalizes_to_zero(r) == zero_in); #endif - - *r = b0; } -static void secp256k1_fe_inv_var(secp256k1_fe *r, const secp256k1_fe *a) { - - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int32_t t[4]; - int32_t d[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t e[9] = { 1, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t f[9] = { 0x3FFFFC2F, 0x3FFFFFFB, 0x3FFFFFFF, 0x3FFFFFFF, 0x3FFFFFFF, - 0x3FFFFFFF, 0x3FFFFFFF, 0x3FFFFFFF, 0xFFFF }; - int32_t g[9]; - secp256k1_fe b; - int i, j, len = 9, sign; - uint32_t eta; - int32_t cond, fn, gn; +static void secp256k1_fe_inv_var(secp256k1_fe *r, const secp256k1_fe *x) { + int sign; + secp256k1_fe tmp; + secp256k1_modinv32_signed30 s; #ifdef VERIFY int zero_in; #endif - b = *a; - secp256k1_fe_normalize(&b); - secp256k1_fe_encode_30(g, &b); - + tmp = *x; + secp256k1_fe_normalize_var(&tmp); #ifdef VERIFY - zero_in = secp256k1_fe_is_zero(&b); + zero_in = secp256k1_fe_normalizes_to_zero(&tmp); #endif - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to - * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */ - eta = -(uint32_t)1; - - for (i = 0; i < 25; ++i) { - - eta = secp256k1_fe_divsteps_30_var(eta, f[0], g[0], t); - secp256k1_fe_update_de_30(d, e, t); - secp256k1_fe_update_fg_30_var(len, f, g, t); - - if (g[0] == 0) { - cond = 0; - for (j = 1; j < len; ++j) { - cond |= g[j]; - } - if (cond == 0) { - break; - } - } - - fn = f[len - 1]; - gn = g[len - 1]; - - cond = ((int32_t)len - 2) >> 31; - cond |= fn ^ (fn >> 31); - cond |= gn ^ (gn >> 31); - - if (cond == 0) { - f[len - 2] |= (uint32_t)fn << 30; - g[len - 2] |= (uint32_t)gn << 30; - --len; - } - } - - VERIFY_CHECK(i < 25); - - /* At this point g is 0 and (if g was not originally 0) f must now equal +/- GCD of - * the initial f, g values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - sign = (f[0] >> 1) & 1; - - secp256k1_fe_decode_30(&b, d); - + secp256k1_fe_to_signed30(&s, &tmp); + secp256k1_modinv32_var(&sign, &s, &secp256k1_const_modinfo_fe); + secp256k1_fe_from_signed30(r, &s); if (sign) { - secp256k1_fe_negate(&b, &b, 1); - secp256k1_fe_normalize_weak(&b); + secp256k1_fe_negate(r, r, 1); + secp256k1_fe_normalize_weak(r); } #ifdef VERIFY - VERIFY_CHECK(!secp256k1_fe_normalizes_to_zero(&b) == !zero_in); + VERIFY_CHECK(secp256k1_fe_normalizes_to_zero(r) == zero_in); #endif - - *r = b; } #endif /* SECP256K1_FIELD_REPR_IMPL_H */ diff --git a/src/field_5x52_impl.h b/src/field_5x52_impl.h index 8f7d2907f1..ee834a0a33 100644 --- a/src/field_5x52_impl.h +++ b/src/field_5x52_impl.h @@ -13,6 +13,7 @@ #include "util.h" #include "field.h" +#include "modinv64_impl.h" #if defined(USE_ASM_X86_64) #include "field_5x52_asm_impl.h" @@ -498,20 +499,14 @@ static SECP256K1_INLINE void secp256k1_fe_from_storage(secp256k1_fe *r, const se #endif } -static void secp256k1_fe_decode_62(secp256k1_fe *r, const int64_t *a) { - +static void secp256k1_fe_from_signed62(secp256k1_fe *r, const secp256k1_modinv64_signed62 *a) { const uint64_t M52 = UINT64_MAX >> 12; - const uint64_t a0 = a[0], a1 = a[1], a2 = a[2], a3 = a[3], a4 = a[4]; + const uint64_t a0 = a->v[0], a1 = a->v[1], a2 = a->v[2], a3 = a->v[3], a4 = a->v[4]; uint64_t r0, r1, r2, r3, r4, t; - t = (int64_t)a4 >> 8; + secp256k1_modinv64_signed62_verify(a); - /* a must be in the range [-2^256, 2^256). */ - VERIFY_CHECK(a0 >> 62 == 0); - VERIFY_CHECK(a1 >> 62 == 0); - VERIFY_CHECK(a2 >> 62 == 0); - VERIFY_CHECK(a3 >> 62 == 0); - VERIFY_CHECK(t == 0 || t == -(uint64_t)1); + t = (int64_t)a4 >> 8; /* Add 2P if a4 is "negative". */ r0 = 0xFFFFDFFFFF85EULL & t; @@ -539,8 +534,7 @@ static void secp256k1_fe_decode_62(secp256k1_fe *r, const int64_t *a) { #endif } -static void secp256k1_fe_encode_62(int64_t *r, const secp256k1_fe *a) { - +static void secp256k1_fe_to_signed62(secp256k1_modinv64_signed62 *r, const secp256k1_fe *a) { const uint64_t M62 = UINT64_MAX >> 2; const uint64_t *n = &a->n[0]; const uint64_t a0 = n[0], a1 = n[1], a2 = n[2], a3 = n[3], a4 = n[4]; @@ -549,392 +543,67 @@ static void secp256k1_fe_encode_62(int64_t *r, const secp256k1_fe *a) { VERIFY_CHECK(a->normalized); #endif - r[0] = (a0 | a1 << 52) & M62; - r[1] = (a1 >> 10 | a2 << 42) & M62; - r[2] = (a2 >> 20 | a3 << 32) & M62; - r[3] = (a3 >> 30 | a4 << 22) & M62; - r[4] = a4 >> 40; + r->v[0] = (a0 | a1 << 52) & M62; + r->v[1] = (a1 >> 10 | a2 << 42) & M62; + r->v[2] = (a2 >> 20 | a3 << 32) & M62; + r->v[3] = (a3 >> 30 | a4 << 22) & M62; + r->v[4] = a4 >> 40; } -static uint64_t secp256k1_fe_divsteps_62(uint64_t eta, uint64_t f0, uint64_t g0, int64_t *t) { - - uint64_t u = 1, v = 0, q = 0, r = 1; - uint64_t c1, c2, f = f0, g = g0, x, y, z; - int i; - - for (i = 0; i < 62; ++i) { - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << i); - VERIFY_CHECK((q * f0 + r * g0) == g << i); - - c1 = (int64_t)eta >> 63; - c2 = -(g & 1); - - x = (f ^ c1) - c1; - y = (u ^ c1) - c1; - z = (v ^ c1) - c1; - - g += x & c2; - q += y & c2; - r += z & c2; - - c1 &= c2; - eta = (eta ^ c1) - (c1 + 1); - - f += g & c1; - u += q & c1; - v += r & c1; - - g >>= 1; - u <<= 1; - v <<= 1; - } - - t[0] = (int64_t)u; - t[1] = (int64_t)v; - t[2] = (int64_t)q; - t[3] = (int64_t)r; - - return eta; -} - -static uint64_t secp256k1_fe_divsteps_62_var(uint64_t eta, uint64_t f0, uint64_t g0, int64_t *t) { - uint64_t u = 1, v = 0, q = 0, r = 1; - uint64_t f = f0, g = g0, m, w, x, y, z; - int i = 62, limit, zeros; - - for (;;) { - /* Use a sentinel bit to count zeros only up to i. */ - zeros = secp256k1_ctz64_var(g | (UINT64_MAX << i)); - - g >>= zeros; - u <<= zeros; - v <<= zeros; - eta -= zeros; - i -= zeros; - - if (i <= 0) { - break; - } - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((g & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << (62 - i)); - VERIFY_CHECK((q * f0 + r * g0) == g << (62 - i)); - - if ((int64_t)eta < 0) { - eta = -eta; - x = f; f = g; g = -x; - y = u; u = q; q = -y; - z = v; v = r; r = -z; - - /* Handle up to 6 divsteps at once, subject to eta and i. */ - limit = ((int)eta + 1) > i ? i : ((int)eta + 1); - m = (UINT64_MAX >> (64 - limit)) & 63U; - - w = (f * g * (f * f - 2)) & m; - } else { - /* Handle up to 4 divsteps at once, subject to eta and i. */ - limit = ((int)eta + 1) > i ? i : ((int)eta + 1); - m = (UINT64_MAX >> (64 - limit)) & 15U; - - w = f + (((f + 1) & 4) << 1); - w = (-w * g) & m; - } - - g += f * w; - q += u * w; - r += v * w; - - VERIFY_CHECK((g & m) == 0); - } - - t[0] = (int64_t)u; - t[1] = (int64_t)v; - t[2] = (int64_t)q; - t[3] = (int64_t)r; - - return eta; -} - -static void secp256k1_fe_update_de_62(int64_t *d, int64_t *e, const int64_t *t) { - - /* P == 2^256 - C62 */ - const int64_t C62 = 0x1000003D1LL; - /* I62 == -P^-1 mod 2^62 */ - const uint64_t I62 = 0x1838091DD2253531LL; - const int64_t M62 = (int64_t)(UINT64_MAX >> 2); - const int64_t d0 = d[0], d1 = d[1], d2 = d[2], d3 = d[3], d4 = d[4]; - const int64_t e0 = e[0], e1 = e[1], e2 = e[2], e3 = e[3], e4 = e[4]; - const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; - int64_t md, me; - int128_t cd, ce; - - cd = (int128_t)u * d0 + (int128_t)v * e0; - ce = (int128_t)q * d0 + (int128_t)r * e0; - - /* Calculate the multiples of P to add, to zero the 62 bottom bits. We choose md, me - * from the centred range [-2^61, 2^61) to keep d, e within [-2^256, 2^256). */ - md = ((int64_t)(I62 * 4 * (uint64_t)cd)) >> 2; - me = ((int64_t)(I62 * 4 * (uint64_t)ce)) >> 2; - - cd -= (int128_t)C62 * md; - ce -= (int128_t)C62 * me; - - VERIFY_CHECK(((int64_t)cd & M62) == 0); cd >>= 62; - VERIFY_CHECK(((int64_t)ce & M62) == 0); ce >>= 62; - - cd += (int128_t)u * d1 + (int128_t)v * e1; - ce += (int128_t)q * d1 + (int128_t)r * e1; - - d[0] = (int64_t)cd & M62; cd >>= 62; - e[0] = (int64_t)ce & M62; ce >>= 62; - - cd += (int128_t)u * d2 + (int128_t)v * e2; - ce += (int128_t)q * d2 + (int128_t)r * e2; - - d[1] = (int64_t)cd & M62; cd >>= 62; - e[1] = (int64_t)ce & M62; ce >>= 62; - - cd += (int128_t)u * d3 + (int128_t)v * e3; - ce += (int128_t)q * d3 + (int128_t)r * e3; - - d[2] = (int64_t)cd & M62; cd >>= 62; - e[2] = (int64_t)ce & M62; ce >>= 62; - - cd += (int128_t)u * d4 + (int128_t)v * e4; - ce += (int128_t)q * d4 + (int128_t)r * e4; - - cd += (int128_t)256 * md; - ce += (int128_t)256 * me; - - d[3] = (int64_t)cd & M62; cd >>= 62; - e[3] = (int64_t)ce & M62; ce >>= 62; - - d[4] = (int64_t)cd; - e[4] = (int64_t)ce; -} - -static void secp256k1_fe_update_fg_62(int64_t *f, int64_t *g, const int64_t *t) { - - const int64_t M62 = (int64_t)(UINT64_MAX >> 2); - const int64_t f0 = f[0], f1 = f[1], f2 = f[2], f3 = f[3], f4 = f[4]; - const int64_t g0 = g[0], g1 = g[1], g2 = g[2], g3 = g[3], g4 = g[4]; - const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; - int128_t cf, cg; - - cf = (int128_t)u * f0 + (int128_t)v * g0; - cg = (int128_t)q * f0 + (int128_t)r * g0; - - VERIFY_CHECK(((int64_t)cf & M62) == 0); cf >>= 62; - VERIFY_CHECK(((int64_t)cg & M62) == 0); cg >>= 62; - - cf += (int128_t)u * f1 + (int128_t)v * g1; - cg += (int128_t)q * f1 + (int128_t)r * g1; +static const secp256k1_modinv64_modinfo secp256k1_const_modinfo_fe = { + {{-0x1000003D1LL, 0, 0, 0, 256}}, + 0x1838091DD2253531LL * 4 +}; - f[0] = (int64_t)cf & M62; cf >>= 62; - g[0] = (int64_t)cg & M62; cg >>= 62; - - cf += (int128_t)u * f2 + (int128_t)v * g2; - cg += (int128_t)q * f2 + (int128_t)r * g2; - - f[1] = (int64_t)cf & M62; cf >>= 62; - g[1] = (int64_t)cg & M62; cg >>= 62; - - cf += (int128_t)u * f3 + (int128_t)v * g3; - cg += (int128_t)q * f3 + (int128_t)r * g3; - - f[2] = (int64_t)cf & M62; cf >>= 62; - g[2] = (int64_t)cg & M62; cg >>= 62; - - cf += (int128_t)u * f4 + (int128_t)v * g4; - cg += (int128_t)q * f4 + (int128_t)r * g4; - - f[3] = (int64_t)cf & M62; cf >>= 62; - g[3] = (int64_t)cg & M62; cg >>= 62; - - f[4] = (int64_t)cf; - g[4] = (int64_t)cg; -} - -static void secp256k1_fe_update_fg_62_var(int len, int64_t *f, int64_t *g, const int64_t *t) { - - const int64_t M62 = (int64_t)(UINT64_MAX >> 2); - const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; - int64_t fi, gi; - int128_t cf, cg; - int i; - - VERIFY_CHECK(len > 0); - - fi = f[0]; - gi = g[0]; - - cf = (int128_t)u * fi + (int128_t)v * gi; - cg = (int128_t)q * fi + (int128_t)r * gi; - - VERIFY_CHECK(((int64_t)cf & M62) == 0); cf >>= 62; - VERIFY_CHECK(((int64_t)cg & M62) == 0); cg >>= 62; - - for (i = 1; i < len; ++i) { - - fi = f[i]; - gi = g[i]; - - cf += (int128_t)u * fi + (int128_t)v * gi; - cg += (int128_t)q * fi + (int128_t)r * gi; - - f[i - 1] = (int64_t)cf & M62; cf >>= 62; - g[i - 1] = (int64_t)cg & M62; cg >>= 62; - } - - f[len - 1] = (int64_t)cf; - g[len - 1] = (int64_t)cg; -} - -static void secp256k1_fe_inv(secp256k1_fe *r, const secp256k1_fe *a) { - - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int64_t t[4]; - int64_t d[5] = { 0, 0, 0, 0, 0 }; - int64_t e[5] = { 1, 0, 0, 0, 0 }; - int64_t f[5] = { 0x3FFFFFFEFFFFFC2FLL, 0x3FFFFFFFFFFFFFFFLL, 0x3FFFFFFFFFFFFFFFLL, - 0x3FFFFFFFFFFFFFFFLL, 0xFFLL }; - int64_t g[5]; - secp256k1_fe b0, b1; - int i, sign; - uint64_t eta; +static void secp256k1_fe_inv(secp256k1_fe *r, const secp256k1_fe *x) { + int sign; + secp256k1_fe tmp; + secp256k1_modinv64_signed62 s; #ifdef VERIFY int zero_in; #endif - b0 = *a; - secp256k1_fe_normalize(&b0); - secp256k1_fe_encode_62(g, &b0); - + tmp = *x; + secp256k1_fe_normalize(&tmp); #ifdef VERIFY - zero_in = secp256k1_fe_is_zero(&b0); + zero_in = secp256k1_fe_normalizes_to_zero(&tmp); #endif - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If the maximum bitlength of g is known to be less than 256, then eta can be set - * initially to -(1 + (256 - maxlen(g))), and only (741 - (256 - maxlen(g))) total - * divsteps are needed. */ - eta = -(uint64_t)1; - - for (i = 0; i < 12; ++i) { - eta = secp256k1_fe_divsteps_62(eta, f[0], g[0], t); - secp256k1_fe_update_de_62(d, e, t); - secp256k1_fe_update_fg_62(f, g, t); - } - - /* At this point sufficient iterations have been performed that g must have reached 0 - * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g - * values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - VERIFY_CHECK((g[0] | g[1] | g[2] | g[3] | g[4]) == 0); - - sign = (f[0] >> 1) & 1; - - secp256k1_fe_decode_62(&b0, d); - - secp256k1_fe_negate(&b1, &b0, 1); - secp256k1_fe_cmov(&b0, &b1, sign); - secp256k1_fe_normalize_weak(&b0); + secp256k1_fe_to_signed62(&s, &tmp); + secp256k1_modinv64(&sign, &s, &secp256k1_const_modinfo_fe); + secp256k1_fe_from_signed62(r, &s); + secp256k1_fe_negate(&tmp, r, 1); + secp256k1_fe_cmov(r, &tmp, sign); + secp256k1_fe_normalize_weak(r); #ifdef VERIFY - VERIFY_CHECK(!secp256k1_fe_normalizes_to_zero(&b0) == !zero_in); + VERIFY_CHECK(secp256k1_fe_normalizes_to_zero(r) == zero_in); #endif - - *r = b0; } -static void secp256k1_fe_inv_var(secp256k1_fe *r, const secp256k1_fe *a) { - - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int64_t t[4]; - int64_t d[5] = { 0, 0, 0, 0, 0 }; - int64_t e[5] = { 1, 0, 0, 0, 0 }; - int64_t f[5] = { 0x3FFFFFFEFFFFFC2FLL, 0x3FFFFFFFFFFFFFFFLL, 0x3FFFFFFFFFFFFFFFLL, - 0x3FFFFFFFFFFFFFFFLL, 0xFFLL }; - int64_t g[5]; - secp256k1_fe b; - int i, j, len = 5, sign; - uint64_t eta; - int64_t cond, fn, gn; +static void secp256k1_fe_inv_var(secp256k1_fe *r, const secp256k1_fe *x) { + int sign; + secp256k1_fe tmp; + secp256k1_modinv64_signed62 s; #ifdef VERIFY int zero_in; #endif - b = *a; - secp256k1_fe_normalize(&b); - secp256k1_fe_encode_62(g, &b); - + tmp = *x; + secp256k1_fe_normalize_var(&tmp); #ifdef VERIFY - zero_in = secp256k1_fe_is_zero(&b); + zero_in = secp256k1_fe_normalizes_to_zero(&tmp); #endif - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to - * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */ - eta = -(uint64_t)1; - - for (i = 0; i < 12; ++i) { - - eta = secp256k1_fe_divsteps_62_var(eta, f[0], g[0], t); - secp256k1_fe_update_de_62(d, e, t); - secp256k1_fe_update_fg_62_var(len, f, g, t); - - if (g[0] == 0) { - cond = 0; - for (j = 1; j < len; ++j) { - cond |= g[j]; - } - if (cond == 0) { - break; - } - } - - fn = f[len - 1]; - gn = g[len - 1]; - - cond = ((int64_t)len - 2) >> 63; - cond |= fn ^ (fn >> 63); - cond |= gn ^ (gn >> 63); - - if (cond == 0) { - f[len - 2] |= fn << 62; - g[len - 2] |= gn << 62; - --len; - } - } - - VERIFY_CHECK(i < 12); - - /* At this point g is 0 and (if g was not originally 0) f must now equal +/- GCD of - * the initial f, g values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - sign = (f[0] >> 1) & 1; - - secp256k1_fe_decode_62(&b, d); - + secp256k1_fe_to_signed62(&s, &tmp); + secp256k1_modinv64_var(&sign, &s, &secp256k1_const_modinfo_fe); + secp256k1_fe_from_signed62(r, &s); if (sign) { - secp256k1_fe_negate(&b, &b, 1); - secp256k1_fe_normalize_weak(&b); + secp256k1_fe_negate(r, r, 1); + secp256k1_fe_normalize_weak(r); } #ifdef VERIFY - VERIFY_CHECK(!secp256k1_fe_normalizes_to_zero(&b) == !zero_in); + VERIFY_CHECK(secp256k1_fe_normalizes_to_zero(r) == zero_in); #endif - - *r = b; } #endif /* SECP256K1_FIELD_REPR_IMPL_H */ diff --git a/src/modinv32.h b/src/modinv32.h new file mode 100644 index 0000000000..a0f6f00467 --- /dev/null +++ b/src/modinv32.h @@ -0,0 +1,32 @@ +/********************************************************************** + * Copyright (c) 2020 Peter Dettman * + * Distributed under the MIT software license, see the accompanying * + * file COPYING or http://www.opensource.org/licenses/mit-license.php.* + **********************************************************************/ + +#ifndef SECP256K1_MODINV32_H +#define SECP256K1_MODINV32_H + +#if defined HAVE_CONFIG_H +#include "libsecp256k1-config.h" +#endif + +#include "util.h" + +typedef struct { + int32_t v[9]; +} secp256k1_modinv32_signed30; + +typedef struct { + /* The modulus in signed30 notation. */ + secp256k1_modinv32_signed30 modulus; + + /* (-modulus^{-1} mod 2^30) * 4 */ + uint32_t montmul4; +} secp256k1_modinv32_modinfo; + +static void secp256k1_modinv32_signed30_verify(const secp256k1_modinv32_signed30* a); +static void secp256k1_modinv32(int* sign, secp256k1_modinv32_signed30 *x, const secp256k1_modinv32_modinfo* modinfo); +static void secp256k1_modinv32_var(int* sign, secp256k1_modinv32_signed30 *x, const secp256k1_modinv32_modinfo* modinfo); + +#endif /* SECP256K1_MODINV32_H */ diff --git a/src/modinv32_impl.h b/src/modinv32_impl.h new file mode 100644 index 0000000000..23c150780e --- /dev/null +++ b/src/modinv32_impl.h @@ -0,0 +1,335 @@ +/********************************************************************** + * Copyright (c) 2020 Peter Dettman * + * Distributed under the MIT software license, see the accompanying * + * file COPYING or http://www.opensource.org/licenses/mit-license.php.* + **********************************************************************/ + +#ifndef SECP256K1_MODINV32_IMPL_H +#define SECP256K1_MODINV32_IMPL_H + +#include "modinv32.h" + +#include "util.h" + +static void secp256k1_modinv32_signed30_verify(const secp256k1_modinv32_signed30* a) { + /* a must be in the range [-2^256, 2^256). */ + VERIFY_CHECK(a->v[0] >> 30 == 0); + VERIFY_CHECK(a->v[1] >> 30 == 0); + VERIFY_CHECK(a->v[2] >> 30 == 0); + VERIFY_CHECK(a->v[3] >> 30 == 0); + VERIFY_CHECK(a->v[4] >> 30 == 0); + VERIFY_CHECK(a->v[5] >> 30 == 0); + VERIFY_CHECK(a->v[6] >> 30 == 0); + VERIFY_CHECK(a->v[7] >> 30 == 0); + VERIFY_CHECK(a->v[8] >> 16 == 0 || a->v[8] >> 16 == -(int32_t)1); +} + +static uint32_t secp256k1_modinv32_divsteps_30(uint32_t eta, uint32_t f0, uint32_t g0, int32_t *t) { + uint32_t u = 1, v = 0, q = 0, r = 1; + uint32_t c1, c2, f = f0, g = g0, x, y, z; + int i; + + for (i = 0; i < 30; ++i) { + VERIFY_CHECK((f & 1) == 1); + VERIFY_CHECK((u * f0 + v * g0) == f << i); + VERIFY_CHECK((q * f0 + r * g0) == g << i); + + c1 = (int32_t)eta >> 31; + c2 = -(g & 1); + + x = (f ^ c1) - c1; + y = (u ^ c1) - c1; + z = (v ^ c1) - c1; + + g += x & c2; + q += y & c2; + r += z & c2; + + c1 &= c2; + eta = (eta ^ c1) - (c1 + 1); + + f += g & c1; + u += q & c1; + v += r & c1; + + g >>= 1; + u <<= 1; + v <<= 1; + } + + t[0] = (int32_t)u; + t[1] = (int32_t)v; + t[2] = (int32_t)q; + t[3] = (int32_t)r; + + return eta; +} + +static uint32_t secp256k1_modinv32_divsteps_30_var(uint32_t eta, uint32_t f0, uint32_t g0, int32_t *t) { + /* inv256[i] = -(2*i+1)^-1 (mod 256) */ + static const uint8_t inv256[128] = { + 0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59, + 0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31, + 0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89, + 0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61, + 0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9, + 0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91, + 0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9, + 0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1, + 0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19, + 0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1, + 0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01 + }; + + uint32_t u = 1, v = 0, q = 0, r = 1; + uint32_t f = f0, g = g0, m, w, x, y, z; + int i = 30, limit, zeros; + + for (;;) { + /* Use a sentinel bit to count zeros only up to i. */ + zeros = secp256k1_ctz32_var(g | (UINT32_MAX << i)); + + g >>= zeros; + u <<= zeros; + v <<= zeros; + eta -= zeros; + i -= zeros; + + if (i <= 0) { + break; + } + + VERIFY_CHECK((f & 1) == 1); + VERIFY_CHECK((g & 1) == 1); + VERIFY_CHECK((u * f0 + v * g0) == f << (30 - i)); + VERIFY_CHECK((q * f0 + r * g0) == g << (30 - i)); + + if ((int32_t)eta < 0) { + eta = -eta; + x = f; f = g; g = -x; + y = u; u = q; q = -y; + z = v; v = r; r = -z; + } + + /* Handle up to 8 divsteps at once, subject to eta and i. */ + limit = ((int)eta + 1) > i ? i : ((int)eta + 1); + m = (UINT32_MAX >> (32 - limit)) & 255U; + + w = (g * inv256[(f >> 1) & 127]) & m; + + g += f * w; + q += u * w; + r += v * w; + + VERIFY_CHECK((g & m) == 0); + } + + t[0] = (int32_t)u; + t[1] = (int32_t)v; + t[2] = (int32_t)q; + t[3] = (int32_t)r; + + return eta; +} + +static void secp256k1_modinv32_update_de_30(secp256k1_modinv32_signed30 *d, secp256k1_modinv32_signed30 *e, const int32_t *t, const secp256k1_modinv32_modinfo* modinfo) { + const int32_t M30 = (int32_t)(UINT32_MAX >> 2); + int32_t u = t[0], v = t[1], q = t[2], r = t[3], di, ei, md, me; + int64_t cd, ce; + int i; + + di = d->v[0]; + ei = e->v[0]; + + cd = (int64_t)u * di + (int64_t)v * ei; + ce = (int64_t)q * di + (int64_t)r * ei; + + /* Calculate the multiples of P to add, to zero the 30 bottom bits. We choose md, me + * from the centred range [-2^29, 2^29) to keep d, e within [-2^256, 2^256). */ + md = ((int32_t)(modinfo->montmul4 * (uint32_t)cd)) >> 2; + me = ((int32_t)(modinfo->montmul4 * (uint32_t)ce)) >> 2; + + if (modinfo->modulus.v[0]) { + cd += (int64_t)modinfo->modulus.v[0] * md; + ce += (int64_t)modinfo->modulus.v[0] * me; + } + + VERIFY_CHECK(((int32_t)cd & M30) == 0); cd >>= 30; + VERIFY_CHECK(((int32_t)ce & M30) == 0); ce >>= 30; + + for (i = 1; i < 9; ++i) { + di = d->v[i]; + ei = e->v[i]; + + cd += (int64_t)u * di + (int64_t)v * ei; + ce += (int64_t)q * di + (int64_t)r * ei; + + if (modinfo->modulus.v[i]) { + cd += (int64_t)modinfo->modulus.v[i] * md; + ce += (int64_t)modinfo->modulus.v[i] * me; + } + + d->v[i - 1] = (int32_t)cd & M30; cd >>= 30; + e->v[i - 1] = (int32_t)ce & M30; ce >>= 30; + } + + d->v[8] = (int32_t)cd; + e->v[8] = (int32_t)ce; +} + +static void secp256k1_modinv32_update_fg_30(secp256k1_modinv32_signed30 *f, secp256k1_modinv32_signed30 *g, const int32_t *t) { + const int32_t M30 = (int32_t)(UINT32_MAX >> 2); + int32_t u = t[0], v = t[1], q = t[2], r = t[3], fi, gi; + int64_t cf, cg; + int i; + + fi = f->v[0]; + gi = g->v[0]; + + cf = (int64_t)u * fi + (int64_t)v * gi; + cg = (int64_t)q * fi + (int64_t)r * gi; + + VERIFY_CHECK(((int32_t)cf & M30) == 0); + VERIFY_CHECK(((int32_t)cg & M30) == 0); + + cf >>= 30; + cg >>= 30; + + for (i = 1; i < 9; ++i) { + fi = f->v[i]; + gi = g->v[i]; + + cf += (int64_t)u * fi + (int64_t)v * gi; + cg += (int64_t)q * fi + (int64_t)r * gi; + + f->v[i - 1] = (int32_t)cf & M30; cf >>= 30; + g->v[i - 1] = (int32_t)cg & M30; cg >>= 30; + } + + f->v[8] = (int32_t)cf; + g->v[8] = (int32_t)cg; +} + +static void secp256k1_modinv32_update_fg_30_var(int len, secp256k1_modinv32_signed30 *f, secp256k1_modinv32_signed30 *g, const int32_t *t) { + const int32_t M30 = (int32_t)(UINT32_MAX >> 2); + const int32_t u = t[0], v = t[1], q = t[2], r = t[3]; + int32_t fi, gi; + int64_t cf, cg; + int i; + + VERIFY_CHECK(len > 0); + + fi = f->v[0]; + gi = g->v[0]; + + cf = (int64_t)u * fi + (int64_t)v * gi; + cg = (int64_t)q * fi + (int64_t)r * gi; + + VERIFY_CHECK(((int32_t)cf & M30) == 0); cf >>= 30; + VERIFY_CHECK(((int32_t)cg & M30) == 0); cg >>= 30; + + for (i = 1; i < len; ++i) { + fi = f->v[i]; + gi = g->v[i]; + + cf += (int64_t)u * fi + (int64_t)v * gi; + cg += (int64_t)q * fi + (int64_t)r * gi; + + f->v[i - 1] = (int32_t)cf & M30; cf >>= 30; + g->v[i - 1] = (int32_t)cg & M30; cg >>= 30; + } + + f->v[len - 1] = (int32_t)cf; + g->v[len - 1] = (int32_t)cg; +} + +static void secp256k1_modinv32(int* sign, secp256k1_modinv32_signed30* x, const secp256k1_modinv32_modinfo* modinfo) { + /* Modular inversion based on the paper "Fast constant-time gcd computation and + * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ + int32_t t[4]; + secp256k1_modinv32_signed30 d = {{0}}; + secp256k1_modinv32_signed30 e = {{1}}; + secp256k1_modinv32_signed30 f = modinfo->modulus; + secp256k1_modinv32_signed30 g = *x; + int i; + uint32_t eta; + + /* The paper uses 'delta'; eta == -delta (a performance tweak). + * + * If the maximum bitlength of g is known to be less than 256, then eta can be set + * initially to -(1 + (256 - maxlen(g))), and only (741 - (256 - maxlen(g))) total + * divsteps are needed. */ + eta = -(uint32_t)1; + + for (i = 0; i < 25; ++i) { + eta = secp256k1_modinv32_divsteps_30(eta, f.v[0], g.v[0], t); + secp256k1_modinv32_update_de_30(&d, &e, t, modinfo); + secp256k1_modinv32_update_fg_30(&f, &g, t); + } + + /* At this point sufficient iterations have been performed that g must have reached 0 + * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g + * values i.e. +/- 1, and d now contains +/- the modular inverse. */ + VERIFY_CHECK((g.v[0] | g.v[1] | g.v[2] | g.v[3] | g.v[4] | g.v[5] | g.v[6] | g.v[7] | g.v[8]) == 0); + + *sign = (f.v[0] >> 1) & 1; + *x = d; +} + +static void secp256k1_modinv32_var(int* sign, secp256k1_modinv32_signed30 *x, const secp256k1_modinv32_modinfo* modinfo) { + /* Modular inversion based on the paper "Fast constant-time gcd computation and + * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ + int32_t t[4]; + secp256k1_modinv32_signed30 d = {{0, 0, 0, 0, 0, 0, 0, 0, 0}}; + secp256k1_modinv32_signed30 e = {{1, 0, 0, 0, 0, 0, 0, 0, 0}}; + secp256k1_modinv32_signed30 f = modinfo->modulus; + secp256k1_modinv32_signed30 g = *x; + int i, j, len = 9; + uint32_t eta; + int32_t cond, fn, gn; + + /* The paper uses 'delta'; eta == -delta (a performance tweak). + * + * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to + * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */ + eta = -(uint32_t)1; + + for (i = 0; i < 25; ++i) { + eta = secp256k1_modinv32_divsteps_30_var(eta, f.v[0], g.v[0], t); + secp256k1_modinv32_update_de_30(&d, &e, t, modinfo); + secp256k1_modinv32_update_fg_30_var(len, &f, &g, t); + + if (g.v[0] == 0) { + cond = 0; + for (j = 1; j < len; ++j) { + cond |= g.v[j]; + } + if (cond == 0) { + break; + } + } + + fn = f.v[len - 1]; + gn = g.v[len - 1]; + + cond = ((int32_t)len - 2) >> 31; + cond |= fn ^ (fn >> 31); + cond |= gn ^ (gn >> 31); + + if (cond == 0) { + f.v[len - 2] |= fn << 30; + g.v[len - 2] |= gn << 30; + --len; + } + } + + /* At this point g is 0 and (if g was not originally 0) f must now equal +/- GCD of + * the initial f, g values i.e. +/- 1, and d now contains +/- the modular inverse. */ + + VERIFY_CHECK(i < 25); + + *sign = (f.v[0] >> 1) & 1; + *x = d; +} + +#endif /* SECP256K1_MODINV32_IMPL_H */ diff --git a/src/modinv64.h b/src/modinv64.h new file mode 100644 index 0000000000..c48ca61fd9 --- /dev/null +++ b/src/modinv64.h @@ -0,0 +1,36 @@ +/********************************************************************** + * Copyright (c) 2020 Peter Dettman * + * Distributed under the MIT software license, see the accompanying * + * file COPYING or http://www.opensource.org/licenses/mit-license.php.* + **********************************************************************/ + +#ifndef SECP256K1_MODINV64_H +#define SECP256K1_MODINV64_H + +#if defined HAVE_CONFIG_H +#include "libsecp256k1-config.h" +#endif + +#include "util.h" + +#ifndef SECP256K1_WIDEMUL_INT128 +#error "modinv64 requires 128-bit wide multiplication support" +#endif + +typedef struct { + int64_t v[5]; +} secp256k1_modinv64_signed62; + +typedef struct { + /* The modulus in signed62 notation. */ + secp256k1_modinv64_signed62 modulus; + + /* (-modulus^{-1} mod 2^62) * 4 */ + uint64_t montmul4; +} secp256k1_modinv64_modinfo; + +static void secp256k1_modinv64_signed62_verify(const secp256k1_modinv64_signed62* a); +static void secp256k1_modinv64(int* sign, secp256k1_modinv64_signed62 *x, const secp256k1_modinv64_modinfo* modinfo); +static void secp256k1_modinv64_var(int* sign, secp256k1_modinv64_signed62 *x, const secp256k1_modinv64_modinfo* modinfo); + +#endif /* SECP256K1_MODINV64_H */ diff --git a/src/modinv64_impl.h b/src/modinv64_impl.h new file mode 100644 index 0000000000..cc6def7bfa --- /dev/null +++ b/src/modinv64_impl.h @@ -0,0 +1,360 @@ +/********************************************************************** + * Copyright (c) 2020 Peter Dettman * + * Distributed under the MIT software license, see the accompanying * + * file COPYING or http://www.opensource.org/licenses/mit-license.php.* + **********************************************************************/ + + +#ifndef SECP256K1_MODINV64_IMPL_H +#define SECP256K1_MODINV64_IMPL_H + +#include "modinv64.h" + +#include "util.h" + +static void secp256k1_modinv64_signed62_verify(const secp256k1_modinv64_signed62* a) { + /* a must be in the range [-2^256, 2^256). */ + VERIFY_CHECK(a->v[0] >> 62 == 0); + VERIFY_CHECK(a->v[1] >> 62 == 0); + VERIFY_CHECK(a->v[2] >> 62 == 0); + VERIFY_CHECK(a->v[3] >> 62 == 0); + VERIFY_CHECK(a->v[4] >> 8 == 0 || a->v[4] >> 8 == -(int64_t)1); +} + +static uint64_t secp256k1_modinv64_divsteps_62(uint64_t eta, uint64_t f0, uint64_t g0, int64_t *t) { + + uint64_t u = 1, v = 0, q = 0, r = 1; + uint64_t c1, c2, f = f0, g = g0, x, y, z; + int i; + + for (i = 0; i < 62; ++i) { + + VERIFY_CHECK((f & 1) == 1); + VERIFY_CHECK((u * f0 + v * g0) == f << i); + VERIFY_CHECK((q * f0 + r * g0) == g << i); + + c1 = (int64_t)eta >> 63; + c2 = -(g & 1); + + x = (f ^ c1) - c1; + y = (u ^ c1) - c1; + z = (v ^ c1) - c1; + + g += x & c2; + q += y & c2; + r += z & c2; + + c1 &= c2; + eta = (eta ^ c1) - (c1 + 1); + + f += g & c1; + u += q & c1; + v += r & c1; + + g >>= 1; + u <<= 1; + v <<= 1; + } + + t[0] = (int64_t)u; + t[1] = (int64_t)v; + t[2] = (int64_t)q; + t[3] = (int64_t)r; + + return eta; +} + +static uint64_t secp256k1_modinv64_divsteps_62_var(uint64_t eta, uint64_t f0, uint64_t g0, int64_t *t) { + uint64_t u = 1, v = 0, q = 0, r = 1; + uint64_t f = f0, g = g0, m, w, x, y, z; + int i = 62, limit, zeros; + + for (;;) { + /* Use a sentinel bit to count zeros only up to i. */ + zeros = secp256k1_ctz64_var(g | (UINT64_MAX << i)); + + g >>= zeros; + u <<= zeros; + v <<= zeros; + eta -= zeros; + i -= zeros; + + if (i <= 0) { + break; + } + + VERIFY_CHECK((f & 1) == 1); + VERIFY_CHECK((g & 1) == 1); + VERIFY_CHECK((u * f0 + v * g0) == f << (62 - i)); + VERIFY_CHECK((q * f0 + r * g0) == g << (62 - i)); + + if ((int64_t)eta < 0) { + eta = -eta; + x = f; f = g; g = -x; + y = u; u = q; q = -y; + z = v; v = r; r = -z; + + /* Handle up to 6 divsteps at once, subject to eta and i. */ + limit = ((int)eta + 1) > i ? i : ((int)eta + 1); + m = (UINT64_MAX >> (64 - limit)) & 63U; + + w = (f * g * (f * f - 2)) & m; + } else { + /* Handle up to 4 divsteps at once, subject to eta and i. */ + limit = ((int)eta + 1) > i ? i : ((int)eta + 1); + m = (UINT64_MAX >> (64 - limit)) & 15U; + + w = f + (((f + 1) & 4) << 1); + w = (-w * g) & m; + } + + g += f * w; + q += u * w; + r += v * w; + + VERIFY_CHECK((g & m) == 0); + } + + t[0] = (int64_t)u; + t[1] = (int64_t)v; + t[2] = (int64_t)q; + t[3] = (int64_t)r; + + return eta; +} + +static void secp256k1_modinv64_update_de_62(secp256k1_modinv64_signed62 *d, secp256k1_modinv64_signed62 *e, const int64_t *t, const secp256k1_modinv64_modinfo* modinfo) { + const int64_t M62 = (int64_t)(UINT64_MAX >> 2); + const int64_t d0 = d->v[0], d1 = d->v[1], d2 = d->v[2], d3 = d->v[3], d4 = d->v[4]; + const int64_t e0 = e->v[0], e1 = e->v[1], e2 = e->v[2], e3 = e->v[3], e4 = e->v[4]; + const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; + int64_t md, me; + int128_t cd, ce; + + cd = (int128_t)u * d0 + (int128_t)v * e0; + ce = (int128_t)q * d0 + (int128_t)r * e0; + + /* Calculate the multiples of P to add, to zero the 62 bottom bits. We choose md, me + * from the centred range [-2^61, 2^61) to keep d, e within [-2^256, 2^256). */ + md = ((int64_t)(modinfo->montmul4 * (uint64_t)cd)) >> 2; + me = ((int64_t)(modinfo->montmul4 * (uint64_t)ce)) >> 2; + + if (modinfo->modulus.v[0]) { + cd += (int128_t)modinfo->modulus.v[0] * md; + ce += (int128_t)modinfo->modulus.v[0] * me; + } + + VERIFY_CHECK(((int64_t)cd & M62) == 0); cd >>= 62; + VERIFY_CHECK(((int64_t)ce & M62) == 0); ce >>= 62; + + cd += (int128_t)u * d1 + (int128_t)v * e1; + ce += (int128_t)q * d1 + (int128_t)r * e1; + + if (modinfo->modulus.v[1]) { + cd += (int128_t)modinfo->modulus.v[1] * md; + ce += (int128_t)modinfo->modulus.v[1] * me; + } + + d->v[0] = (int64_t)cd & M62; cd >>= 62; + e->v[0] = (int64_t)ce & M62; ce >>= 62; + + cd += (int128_t)u * d2 + (int128_t)v * e2; + ce += (int128_t)q * d2 + (int128_t)r * e2; + + if (modinfo->modulus.v[2]) { + cd += (int128_t)modinfo->modulus.v[2] * md; + ce += (int128_t)modinfo->modulus.v[2] * me; + } + + d->v[1] = (int64_t)cd & M62; cd >>= 62; + e->v[1] = (int64_t)ce & M62; ce >>= 62; + + cd += (int128_t)u * d3 + (int128_t)v * e3; + ce += (int128_t)q * d3 + (int128_t)r * e3; + + if (modinfo->modulus.v[3]) { + cd += (int128_t)modinfo->modulus.v[3] * md; + ce += (int128_t)modinfo->modulus.v[3] * me; + } + + d->v[2] = (int64_t)cd & M62; cd >>= 62; + e->v[2] = (int64_t)ce & M62; ce >>= 62; + + cd += (int128_t)u * d4 + (int128_t)v * e4; + ce += (int128_t)q * d4 + (int128_t)r * e4; + + if (modinfo->modulus.v[4]) { + cd += (int128_t)modinfo->modulus.v[4] * md; + ce += (int128_t)modinfo->modulus.v[4] * me; + } + + d->v[3] = (int64_t)cd & M62; cd >>= 62; + e->v[3] = (int64_t)ce & M62; ce >>= 62; + + d->v[4] = (int64_t)cd; + e->v[4] = (int64_t)ce; +} + +static void secp256k1_modinv64_update_fg_62(secp256k1_modinv64_signed62 *f, secp256k1_modinv64_signed62 *g, const int64_t *t) { + const int64_t M62 = (int64_t)(UINT64_MAX >> 2); + const int64_t f0 = f->v[0], f1 = f->v[1], f2 = f->v[2], f3 = f->v[3], f4 = f->v[4]; + const int64_t g0 = g->v[0], g1 = g->v[1], g2 = g->v[2], g3 = g->v[3], g4 = g->v[4]; + const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; + int128_t cf, cg; + + cf = (int128_t)u * f0 + (int128_t)v * g0; + cg = (int128_t)q * f0 + (int128_t)r * g0; + + VERIFY_CHECK(((int64_t)cf & M62) == 0); cf >>= 62; + VERIFY_CHECK(((int64_t)cg & M62) == 0); cg >>= 62; + + cf += (int128_t)u * f1 + (int128_t)v * g1; + cg += (int128_t)q * f1 + (int128_t)r * g1; + + f->v[0] = (int64_t)cf & M62; cf >>= 62; + g->v[0] = (int64_t)cg & M62; cg >>= 62; + + cf += (int128_t)u * f2 + (int128_t)v * g2; + cg += (int128_t)q * f2 + (int128_t)r * g2; + + f->v[1] = (int64_t)cf & M62; cf >>= 62; + g->v[1] = (int64_t)cg & M62; cg >>= 62; + + cf += (int128_t)u * f3 + (int128_t)v * g3; + cg += (int128_t)q * f3 + (int128_t)r * g3; + + f->v[2] = (int64_t)cf & M62; cf >>= 62; + g->v[2] = (int64_t)cg & M62; cg >>= 62; + + cf += (int128_t)u * f4 + (int128_t)v * g4; + cg += (int128_t)q * f4 + (int128_t)r * g4; + + f->v[3] = (int64_t)cf & M62; cf >>= 62; + g->v[3] = (int64_t)cg & M62; cg >>= 62; + + f->v[4] = (int64_t)cf; + g->v[4] = (int64_t)cg; +} + +static void secp256k1_modinv64_update_fg_62_var(int len, secp256k1_modinv64_signed62 *f, secp256k1_modinv64_signed62 *g, const int64_t *t) { + const int64_t M62 = (int64_t)(UINT64_MAX >> 2); + const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; + int64_t fi, gi; + int128_t cf, cg; + int i; + + VERIFY_CHECK(len > 0); + + fi = f->v[0]; + gi = g->v[0]; + + cf = (int128_t)u * fi + (int128_t)v * gi; + cg = (int128_t)q * fi + (int128_t)r * gi; + + VERIFY_CHECK(((int64_t)cf & M62) == 0); cf >>= 62; + VERIFY_CHECK(((int64_t)cg & M62) == 0); cg >>= 62; + + for (i = 1; i < len; ++i) { + fi = f->v[i]; + gi = g->v[i]; + + cf += (int128_t)u * fi + (int128_t)v * gi; + cg += (int128_t)q * fi + (int128_t)r * gi; + + f->v[i - 1] = (int64_t)cf & M62; cf >>= 62; + g->v[i - 1] = (int64_t)cg & M62; cg >>= 62; + } + + f->v[len - 1] = (int64_t)cf; + g->v[len - 1] = (int64_t)cg; +} + +static void secp256k1_modinv64(int* sign, secp256k1_modinv64_signed62* x, const secp256k1_modinv64_modinfo* modinfo) { + /* Modular inversion based on the paper "Fast constant-time gcd computation and + * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ + + int64_t t[4]; + secp256k1_modinv64_signed62 d = {{0, 0, 0, 0, 0}}; + secp256k1_modinv64_signed62 e = {{1, 0, 0, 0, 0}}; + secp256k1_modinv64_signed62 f = modinfo->modulus; + secp256k1_modinv64_signed62 g = *x; + int i; + uint64_t eta; + + /* The paper uses 'delta'; eta == -delta (a performance tweak). + * + * If the maximum bitlength of g is known to be less than 256, then eta can be set + * initially to -(1 + (256 - maxlen(g))), and only (741 - (256 - maxlen(g))) total + * divsteps are needed. */ + eta = -(uint64_t)1; + + for (i = 0; i < 12; ++i) { + eta = secp256k1_modinv64_divsteps_62(eta, f.v[0], g.v[0], t); + secp256k1_modinv64_update_de_62(&d, &e, t, modinfo); + secp256k1_modinv64_update_fg_62(&f, &g, t); + } + + /* At this point sufficient iterations have been performed that g must have reached 0 + * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g + * values i.e. +/- 1, and d now contains +/- the modular inverse. */ + VERIFY_CHECK((g.v[0] | g.v[1] | g.v[2] | g.v[3] | g.v[4]) == 0); + + *sign = (f.v[0] >> 1) & 1; + *x = d; +} + +static void secp256k1_modinv64_var(int* sign, secp256k1_modinv64_signed62 *x, const secp256k1_modinv64_modinfo* modinfo) { + /* Modular inversion based on the paper "Fast constant-time gcd computation and + * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ + + int64_t t[4]; + secp256k1_modinv64_signed62 d = {{0, 0, 0, 0, 0}}; + secp256k1_modinv64_signed62 e = {{1, 0, 0, 0, 0}}; + secp256k1_modinv64_signed62 f = modinfo->modulus; + secp256k1_modinv64_signed62 g = *x; + int i, j, len = 5; + uint64_t eta; + int64_t cond, fn, gn; + + /* The paper uses 'delta'; eta == -delta (a performance tweak). + * + * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to + * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */ + eta = -(uint64_t)1; + + for (i = 0; i < 12; ++i) { + eta = secp256k1_modinv64_divsteps_62_var(eta, f.v[0], g.v[0], t); + secp256k1_modinv64_update_de_62(&d, &e, t, modinfo); + secp256k1_modinv64_update_fg_62_var(len, &f, &g, t); + + if (g.v[0] == 0) { + cond = 0; + for (j = 1; j < len; ++j) { + cond |= g.v[j]; + } + if (cond == 0) { + break; + } + } + + fn = f.v[len - 1]; + gn = g.v[len - 1]; + + cond = ((int64_t)len - 2) >> 63; + cond |= fn ^ (fn >> 63); + cond |= gn ^ (gn >> 63); + + if (cond == 0) { + f.v[len - 2] |= fn << 62; + g.v[len - 2] |= gn << 62; + --len; + } + } + + VERIFY_CHECK(i < 12); + + *sign = (f.v[0] >> 1) & 1; + *x = d; +} + +#endif /* SECP256K1_MODINV64_IMPL_H */ diff --git a/src/scalar_4x64_impl.h b/src/scalar_4x64_impl.h index b649928d0c..96c0fe76cb 100644 --- a/src/scalar_4x64_impl.h +++ b/src/scalar_4x64_impl.h @@ -7,6 +7,8 @@ #ifndef SECP256K1_SCALAR_REPR_IMPL_H #define SECP256K1_SCALAR_REPR_IMPL_H +#include "modinv64_impl.h" + /* Limbs of the secp256k1 order. */ #define SECP256K1_N_0 ((uint64_t)0xBFD25E8CD0364141ULL) #define SECP256K1_N_1 ((uint64_t)0xBAAEDCE6AF48A03BULL) @@ -962,18 +964,12 @@ static const secp256k1_scalar SECP256K1_SCALAR_NEG_TWO_POW_256 = SECP256K1_SCALA 0x755DB9CDUL, 0x5E914077UL, 0x7FA4BD19UL, 0xA06C8282UL ); -static void secp256k1_scalar_decode_62(secp256k1_scalar *r, const int64_t *a) { - - const uint64_t a0 = a[0], a1 = a[1], a2 = a[2], a3 = a[3], a4 = a[4]; +static void secp256k1_scalar_from_signed62(secp256k1_scalar *r, const secp256k1_modinv64_signed62 *a) { + const uint64_t a0 = a->v[0], a1 = a->v[1], a2 = a->v[2], a3 = a->v[3], a4 = a->v[4]; uint64_t r0, r1, r2, r3; secp256k1_scalar u; - /* a must be in the range [-2^256, 2^256). */ - VERIFY_CHECK(a0 >> 62 == 0); - VERIFY_CHECK(a1 >> 62 == 0); - VERIFY_CHECK(a2 >> 62 == 0); - VERIFY_CHECK(a3 >> 62 == 0); - VERIFY_CHECK((int64_t)a4 >> 8 == 0 || (int64_t)a4 >> 8 == -(int64_t)1); + secp256k1_modinv64_signed62_verify(a); r0 = a0 | a1 << 62; r1 = a1 >> 2 | a2 << 60; @@ -991,392 +987,53 @@ static void secp256k1_scalar_decode_62(secp256k1_scalar *r, const int64_t *a) { secp256k1_scalar_cmov(r, &u, a4 >> 63); } -static void secp256k1_scalar_encode_62(int64_t *r, const secp256k1_scalar *a) { - +static void secp256k1_scalar_to_signed62(secp256k1_modinv64_signed62 *r, const secp256k1_scalar *a) { const uint64_t M62 = UINT64_MAX >> 2; const uint64_t *d = &a->d[0]; const uint64_t a0 = d[0], a1 = d[1], a2 = d[2], a3 = d[3]; -#ifdef VERIFY - VERIFY_CHECK(secp256k1_scalar_check_overflow(a) == 0); -#endif - - r[0] = a0 & M62; - r[1] = (a0 >> 62 | a1 << 2) & M62; - r[2] = (a1 >> 60 | a2 << 4) & M62; - r[3] = (a2 >> 58 | a3 << 6) & M62; - r[4] = a3 >> 56; + r->v[0] = a0 & M62; + r->v[1] = (a0 >> 62 | a1 << 2) & M62; + r->v[2] = (a1 >> 60 | a2 << 4) & M62; + r->v[3] = (a2 >> 58 | a3 << 6) & M62; + r->v[4] = a3 >> 56; } -static uint64_t secp256k1_scalar_divsteps_62(uint64_t eta, uint64_t f0, uint64_t g0, int64_t *t) { - - uint64_t u = 1, v = 0, q = 0, r = 1; - uint64_t c1, c2, f = f0, g = g0, x, y, z; - int i; - - for (i = 0; i < 62; ++i) { - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << i); - VERIFY_CHECK((q * f0 + r * g0) == g << i); - - c1 = (int64_t)eta >> 63; - c2 = -(g & 1); - - x = (f ^ c1) - c1; - y = (u ^ c1) - c1; - z = (v ^ c1) - c1; - - g += x & c2; - q += y & c2; - r += z & c2; - - c1 &= c2; - eta = (eta ^ c1) - (c1 + 1); - - f += g & c1; - u += q & c1; - v += r & c1; - - g >>= 1; - u <<= 1; - v <<= 1; - } - - t[0] = (int64_t)u; - t[1] = (int64_t)v; - t[2] = (int64_t)q; - t[3] = (int64_t)r; - - return eta; -} - -static uint64_t secp256k1_scalar_divsteps_62_var(uint64_t eta, uint64_t f0, uint64_t g0, int64_t *t) { - uint64_t u = 1, v = 0, q = 0, r = 1; - uint64_t f = f0, g = g0, m, w, x, y, z; - int i = 62, limit, zeros; - - for (;;) { - /* Use a sentinel bit to count zeros only up to i. */ - zeros = secp256k1_ctz64_var(g | (UINT64_MAX << i)); - - g >>= zeros; - u <<= zeros; - v <<= zeros; - eta -= zeros; - i -= zeros; - - if (i <= 0) { - break; - } - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((g & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << (62 - i)); - VERIFY_CHECK((q * f0 + r * g0) == g << (62 - i)); - - if ((int64_t)eta < 0) { - eta = -eta; - x = f; f = g; g = -x; - y = u; u = q; q = -y; - z = v; v = r; r = -z; - - /* Handle up to 6 divsteps at once, subject to eta and i. */ - limit = ((int)eta + 1) > i ? i : ((int)eta + 1); - m = (UINT64_MAX >> (64 - limit)) & 63U; - - w = (f * g * (f * f - 2)) & m; - } else { - /* Handle up to 4 divsteps at once, subject to eta and i. */ - limit = ((int)eta + 1) > i ? i : ((int)eta + 1); - m = (UINT64_MAX >> (64 - limit)) & 15U; - - w = f + (((f + 1) & 4) << 1); - w = (-w * g) & m; - } - - g += f * w; - q += u * w; - r += v * w; - - VERIFY_CHECK((g & m) == 0); - } - - t[0] = (int64_t)u; - t[1] = (int64_t)v; - t[2] = (int64_t)q; - t[3] = (int64_t)r; - - return eta; -} - -static void secp256k1_scalar_update_de_62(int64_t *d, int64_t *e, const int64_t *t) { - - /* I62 == -P^-1 mod 2^62 */ - const uint64_t I62 = 0x0B0DFF665588B13FLL; - const int64_t M62 = (int64_t)(UINT64_MAX >> 2); - const int64_t P[5] = { 0x3FD25E8CD0364141LL, 0x2ABB739ABD2280EELL, -0x15LL, 0, 256 }; - const int64_t d0 = d[0], d1 = d[1], d2 = d[2], d3 = d[3], d4 = d[4]; - const int64_t e0 = e[0], e1 = e[1], e2 = e[2], e3 = e[3], e4 = e[4]; - const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; - int64_t md, me; - int128_t cd, ce; - - cd = (int128_t)u * d0 + (int128_t)v * e0; - ce = (int128_t)q * d0 + (int128_t)r * e0; - - /* Calculate the multiples of P to add, to zero the 62 bottom bits. We choose md, me - * from the centred range [-2^61, 2^61) to keep d, e within [-2^256, 2^256). */ - md = ((int64_t)(I62 * 4 * (uint64_t)cd)) >> 2; - me = ((int64_t)(I62 * 4 * (uint64_t)ce)) >> 2; - - cd += (int128_t)P[0] * md; - ce += (int128_t)P[0] * me; - - VERIFY_CHECK(((int64_t)cd & M62) == 0); cd >>= 62; - VERIFY_CHECK(((int64_t)ce & M62) == 0); ce >>= 62; - - cd += (int128_t)u * d1 + (int128_t)v * e1; - ce += (int128_t)q * d1 + (int128_t)r * e1; - - cd += (int128_t)P[1] * md; - ce += (int128_t)P[1] * me; - - d[0] = (int64_t)cd & M62; cd >>= 62; - e[0] = (int64_t)ce & M62; ce >>= 62; - - cd += (int128_t)u * d2 + (int128_t)v * e2; - ce += (int128_t)q * d2 + (int128_t)r * e2; - - cd += (int128_t)P[2] * md; - ce += (int128_t)P[2] * me; - - d[1] = (int64_t)cd & M62; cd >>= 62; - e[1] = (int64_t)ce & M62; ce >>= 62; - - cd += (int128_t)u * d3 + (int128_t)v * e3; - ce += (int128_t)q * d3 + (int128_t)r * e3; - - d[2] = (int64_t)cd & M62; cd >>= 62; - e[2] = (int64_t)ce & M62; ce >>= 62; - - cd += (int128_t)u * d4 + (int128_t)v * e4; - ce += (int128_t)q * d4 + (int128_t)r * e4; - - cd += (int128_t)P[4] * md; - ce += (int128_t)P[4] * me; - - d[3] = (int64_t)cd & M62; cd >>= 62; - e[3] = (int64_t)ce & M62; ce >>= 62; - - d[4] = (int64_t)cd; - e[4] = (int64_t)ce; -} - -static void secp256k1_scalar_update_fg_62(int64_t *f, int64_t *g, const int64_t *t) { - - const int64_t M62 = (int64_t)(UINT64_MAX >> 2); - const int64_t f0 = f[0], f1 = f[1], f2 = f[2], f3 = f[3], f4 = f[4]; - const int64_t g0 = g[0], g1 = g[1], g2 = g[2], g3 = g[3], g4 = g[4]; - const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; - int128_t cf, cg; - - cf = (int128_t)u * f0 + (int128_t)v * g0; - cg = (int128_t)q * f0 + (int128_t)r * g0; - - VERIFY_CHECK(((int64_t)cf & M62) == 0); cf >>= 62; - VERIFY_CHECK(((int64_t)cg & M62) == 0); cg >>= 62; - - cf += (int128_t)u * f1 + (int128_t)v * g1; - cg += (int128_t)q * f1 + (int128_t)r * g1; - - f[0] = (int64_t)cf & M62; cf >>= 62; - g[0] = (int64_t)cg & M62; cg >>= 62; - - cf += (int128_t)u * f2 + (int128_t)v * g2; - cg += (int128_t)q * f2 + (int128_t)r * g2; - - f[1] = (int64_t)cf & M62; cf >>= 62; - g[1] = (int64_t)cg & M62; cg >>= 62; - - cf += (int128_t)u * f3 + (int128_t)v * g3; - cg += (int128_t)q * f3 + (int128_t)r * g3; - - f[2] = (int64_t)cf & M62; cf >>= 62; - g[2] = (int64_t)cg & M62; cg >>= 62; - - cf += (int128_t)u * f4 + (int128_t)v * g4; - cg += (int128_t)q * f4 + (int128_t)r * g4; - - f[3] = (int64_t)cf & M62; cf >>= 62; - g[3] = (int64_t)cg & M62; cg >>= 62; - - f[4] = (int64_t)cf; - g[4] = (int64_t)cg; -} - -static void secp256k1_scalar_update_fg_62_var(int len, int64_t *f, int64_t *g, const int64_t *t) { - - const int64_t M62 = (int64_t)(UINT64_MAX >> 2); - const int64_t u = t[0], v = t[1], q = t[2], r = t[3]; - int64_t fi, gi; - int128_t cf, cg; - int i; - - VERIFY_CHECK(len > 0); - - fi = f[0]; - gi = g[0]; - - cf = (int128_t)u * fi + (int128_t)v * gi; - cg = (int128_t)q * fi + (int128_t)r * gi; - - VERIFY_CHECK(((int64_t)cf & M62) == 0); cf >>= 62; - VERIFY_CHECK(((int64_t)cg & M62) == 0); cg >>= 62; - - for (i = 1; i < len; ++i) { - - fi = f[i]; - gi = g[i]; - - cf += (int128_t)u * fi + (int128_t)v * gi; - cg += (int128_t)q * fi + (int128_t)r * gi; - - f[i - 1] = (int64_t)cf & M62; cf >>= 62; - g[i - 1] = (int64_t)cg & M62; cg >>= 62; - } - - f[len - 1] = (int64_t)cf; - g[len - 1] = (int64_t)cg; -} +static const secp256k1_modinv64_modinfo secp256k1_const_modinfo_scalar = { + {{0x3FD25E8CD0364141LL, 0x2ABB739ABD2280EELL, -0x15LL, 0, 256}}, + 0x0B0DFF665588B13FLL * 4 +}; static void secp256k1_scalar_inverse(secp256k1_scalar *r, const secp256k1_scalar *x) { - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int64_t t[4]; - int64_t d[5] = { 0, 0, 0, 0, 0 }; - int64_t e[5] = { 1, 0, 0, 0, 0 }; - int64_t f[5] = { 0x3FD25E8CD0364141LL, 0x2ABB739ABD2280EELL, 0x3FFFFFFFFFFFFFEBLL, - 0x3FFFFFFFFFFFFFFFLL, 0xFFLL }; - int64_t g[5]; - secp256k1_scalar b0; - int i, sign; - uint64_t eta; + int sign; + secp256k1_modinv64_signed62 s; #ifdef VERIFY int zero_in = secp256k1_scalar_is_zero(x); #endif - - b0 = *x; - secp256k1_scalar_encode_62(g, &b0); - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If the maximum bitlength of g is known to be less than 256, then eta can be set - * initially to -(1 + (256 - maxlen(g))), and only (741 - (256 - maxlen(g))) total - * divsteps are needed. */ - eta = -(uint64_t)1; - - for (i = 0; i < 12; ++i) { - eta = secp256k1_scalar_divsteps_62(eta, f[0], g[0], t); - secp256k1_scalar_update_de_62(d, e, t); - secp256k1_scalar_update_fg_62(f, g, t); - } - - /* At this point sufficient iterations have been performed that g must have reached 0 - * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g - * values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - VERIFY_CHECK((g[0] | g[1] | g[2] | g[3] | g[4]) == 0); - - sign = (f[0] >> 1) & 1; - - secp256k1_scalar_decode_62(&b0, d); - secp256k1_scalar_cond_negate(&b0, sign); + secp256k1_scalar_to_signed62(&s, x); + secp256k1_modinv64(&sign, &s, &secp256k1_const_modinfo_scalar); + secp256k1_scalar_from_signed62(r, &s); + secp256k1_scalar_cond_negate(r, sign); #ifdef VERIFY - VERIFY_CHECK(!secp256k1_scalar_is_zero(&b0) == !zero_in); + VERIFY_CHECK(secp256k1_scalar_is_zero(r) == zero_in); #endif - - *r = b0; } static void secp256k1_scalar_inverse_var(secp256k1_scalar *r, const secp256k1_scalar *x) { - - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int64_t t[4]; - int64_t d[5] = { 0, 0, 0, 0, 0 }; - int64_t e[5] = { 1, 0, 0, 0, 0 }; - int64_t f[5] = { 0x3FD25E8CD0364141LL, 0x2ABB739ABD2280EELL, 0x3FFFFFFFFFFFFFEBLL, - 0x3FFFFFFFFFFFFFFFLL, 0xFFLL }; - int64_t g[5]; - secp256k1_scalar b; - int i, j, len = 5, sign; - uint64_t eta; - int64_t cond, fn, gn; + int sign; + secp256k1_modinv64_signed62 s; #ifdef VERIFY int zero_in = secp256k1_scalar_is_zero(x); #endif - - b = *x; - secp256k1_scalar_encode_62(g, &b); - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to - * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */ - eta = -(uint64_t)1; - - for (i = 0; i < 12; ++i) { - - eta = secp256k1_scalar_divsteps_62_var(eta, f[0], g[0], t); - secp256k1_scalar_update_de_62(d, e, t); - secp256k1_scalar_update_fg_62_var(len, f, g, t); - - if (g[0] == 0) { - cond = 0; - for (j = 1; j < len; ++j) { - cond |= g[j]; - } - if (cond == 0) { - break; - } - } - - fn = f[len - 1]; - gn = g[len - 1]; - - cond = ((int64_t)len - 2) >> 63; - cond |= fn ^ (fn >> 63); - cond |= gn ^ (gn >> 63); - - if (cond == 0) { - f[len - 2] |= fn << 62; - g[len - 2] |= gn << 62; - --len; - } - } - - VERIFY_CHECK(i < 12); - - /* At this point g is 0 and (if g was not originally 0) f must now equal +/- GCD of - * the initial f, g values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - sign = (f[0] >> 1) & 1; - - secp256k1_scalar_decode_62(&b, d); - - if (sign) { - secp256k1_scalar_negate(&b, &b); - } + secp256k1_scalar_to_signed62(&s, x); + secp256k1_modinv64_var(&sign, &s, &secp256k1_const_modinfo_scalar); + secp256k1_scalar_from_signed62(r, &s); + if (sign) secp256k1_scalar_negate(r, r); #ifdef VERIFY - VERIFY_CHECK(!secp256k1_scalar_is_zero(&b) == !zero_in); + VERIFY_CHECK(secp256k1_scalar_is_zero(r) == zero_in); #endif - - *r = b; } SECP256K1_INLINE static int secp256k1_scalar_is_even(const secp256k1_scalar *a) { diff --git a/src/scalar_8x32_impl.h b/src/scalar_8x32_impl.h index 9fc1c6a82c..6e14eeec0b 100644 --- a/src/scalar_8x32_impl.h +++ b/src/scalar_8x32_impl.h @@ -7,6 +7,8 @@ #ifndef SECP256K1_SCALAR_REPR_IMPL_H #define SECP256K1_SCALAR_REPR_IMPL_H +#include "modinv32_impl.h" + /* Limbs of the secp256k1 order. */ #define SECP256K1_N_0 ((uint32_t)0xD0364141UL) #define SECP256K1_N_1 ((uint32_t)0xBFD25E8CUL) @@ -738,41 +740,21 @@ static const secp256k1_scalar SECP256K1_SCALAR_NEG_TWO_POW_256 = SECP256K1_SCALA 0x755DB9CDUL, 0x5E914077UL, 0x7FA4BD19UL, 0xA06C8282UL ); -static void secp256k1_scalar_decode_30(secp256k1_scalar *r, const int32_t *a) { - - const uint32_t a0 = a[0], a1 = a[1], a2 = a[2], a3 = a[3], a4 = a[4], - a5 = a[5], a6 = a[6], a7 = a[7], a8 = a[8]; - uint32_t r0, r1, r2, r3, r4, r5, r6, r7; +static void secp256k1_scalar_from_signed30(secp256k1_scalar *r, const secp256k1_modinv32_signed30 *a) { + const uint32_t a0 = a->v[0], a1 = a->v[1], a2 = a->v[2], a3 = a->v[3], a4 = a->v[4], + a5 = a->v[5], a6 = a->v[6], a7 = a->v[7], a8 = a->v[8]; secp256k1_scalar u; - /* a must be in the range [-2^256, 2^256). */ - VERIFY_CHECK(a0 >> 30 == 0); - VERIFY_CHECK(a1 >> 30 == 0); - VERIFY_CHECK(a2 >> 30 == 0); - VERIFY_CHECK(a3 >> 30 == 0); - VERIFY_CHECK(a4 >> 30 == 0); - VERIFY_CHECK(a5 >> 30 == 0); - VERIFY_CHECK(a6 >> 30 == 0); - VERIFY_CHECK(a7 >> 30 == 0); - VERIFY_CHECK((int32_t)a8 >> 16 == 0 || (int32_t)a8 >> 16 == -(int32_t)1); - - r0 = a0 | a1 << 30; - r1 = a1 >> 2 | a2 << 28; - r2 = a2 >> 4 | a3 << 26; - r3 = a3 >> 6 | a4 << 24; - r4 = a4 >> 8 | a5 << 22; - r5 = a5 >> 10 | a6 << 20; - r6 = a6 >> 12 | a7 << 18; - r7 = a7 >> 14 | a8 << 16; - - r->d[0] = r0; - r->d[1] = r1; - r->d[2] = r2; - r->d[3] = r3; - r->d[4] = r4; - r->d[5] = r5; - r->d[6] = r6; - r->d[7] = r7; + secp256k1_modinv32_signed30_verify(a); + + r->d[0] = a0 | a1 << 30; + r->d[1] = a1 >> 2 | a2 << 28; + r->d[2] = a2 >> 4 | a3 << 26; + r->d[3] = a3 >> 6 | a4 << 24; + r->d[4] = a4 >> 8 | a5 << 22; + r->d[5] = a5 >> 10 | a6 << 20; + r->d[6] = a6 >> 12 | a7 << 18; + r->d[7] = a7 >> 14 | a8 << 16; secp256k1_scalar_reduce(r, secp256k1_scalar_check_overflow(r)); @@ -780,7 +762,7 @@ static void secp256k1_scalar_decode_30(secp256k1_scalar *r, const int32_t *a) { secp256k1_scalar_cmov(r, &u, a8 >> 31); } -static void secp256k1_scalar_encode_30(int32_t *r, const secp256k1_scalar *a) { +static void secp256k1_scalar_to_signed30(secp256k1_modinv32_signed30 *r, const secp256k1_scalar *a) { const uint32_t M30 = UINT32_MAX >> 2; const uint32_t *d = &a->d[0]; @@ -791,370 +773,52 @@ static void secp256k1_scalar_encode_30(int32_t *r, const secp256k1_scalar *a) { VERIFY_CHECK(secp256k1_scalar_check_overflow(a) == 0); #endif - r[0] = a0 & M30; - r[1] = (a0 >> 30 | a1 << 2) & M30; - r[2] = (a1 >> 28 | a2 << 4) & M30; - r[3] = (a2 >> 26 | a3 << 6) & M30; - r[4] = (a3 >> 24 | a4 << 8) & M30; - r[5] = (a4 >> 22 | a5 << 10) & M30; - r[6] = (a5 >> 20 | a6 << 12) & M30; - r[7] = (a6 >> 18 | a7 << 14) & M30; - r[8] = a7 >> 16; + r->v[0] = a0 & M30; + r->v[1] = (a0 >> 30 | a1 << 2) & M30; + r->v[2] = (a1 >> 28 | a2 << 4) & M30; + r->v[3] = (a2 >> 26 | a3 << 6) & M30; + r->v[4] = (a3 >> 24 | a4 << 8) & M30; + r->v[5] = (a4 >> 22 | a5 << 10) & M30; + r->v[6] = (a5 >> 20 | a6 << 12) & M30; + r->v[7] = (a6 >> 18 | a7 << 14) & M30; + r->v[8] = a7 >> 16; } -static uint32_t secp256k1_scalar_divsteps_30(uint32_t eta, uint32_t f0, uint32_t g0, int32_t *t) { - - uint32_t u = 1, v = 0, q = 0, r = 1; - uint32_t c1, c2, f = f0, g = g0, x, y, z; - int i; - - for (i = 0; i < 30; ++i) { - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << i); - VERIFY_CHECK((q * f0 + r * g0) == g << i); - - c1 = (int32_t)eta >> 31; - c2 = -(g & 1); - - x = (f ^ c1) - c1; - y = (u ^ c1) - c1; - z = (v ^ c1) - c1; - - g += x & c2; - q += y & c2; - r += z & c2; - - c1 &= c2; - eta = (eta ^ c1) - (c1 + 1); - - f += g & c1; - u += q & c1; - v += r & c1; - - g >>= 1; - u <<= 1; - v <<= 1; - } - - t[0] = (int32_t)u; - t[1] = (int32_t)v; - t[2] = (int32_t)q; - t[3] = (int32_t)r; - - return eta; -} - -static uint32_t secp256k1_scalar_divsteps_30_var(uint32_t eta, uint32_t f0, uint32_t g0, int32_t *t) { - static const uint8_t inv256[128] = { - 0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59, - 0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31, - 0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89, - 0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61, - 0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9, - 0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91, - 0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9, - 0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1, - 0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19, - 0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1, - 0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01 - }; - - uint32_t u = 1, v = 0, q = 0, r = 1; - uint32_t f = f0, g = g0, m, w, x, y, z; - int i = 30, limit, zeros; - - for (;;) { - /* Use a sentinel bit to count zeros only up to i. */ - zeros = secp256k1_ctz32_var(g | (UINT32_MAX << i)); - - g >>= zeros; - u <<= zeros; - v <<= zeros; - eta -= zeros; - i -= zeros; - - if (i <= 0) { - break; - } - - VERIFY_CHECK((f & 1) == 1); - VERIFY_CHECK((g & 1) == 1); - VERIFY_CHECK((u * f0 + v * g0) == f << (30 - i)); - VERIFY_CHECK((q * f0 + r * g0) == g << (30 - i)); - - if ((int32_t)eta < 0) { - eta = -eta; - x = f; f = g; g = -x; - y = u; u = q; q = -y; - z = v; v = r; r = -z; - } - - /* Handle up to 8 divsteps at once, subject to eta and i. */ - limit = ((int)eta + 1) > i ? i : ((int)eta + 1); - m = (UINT32_MAX >> (32 - limit)) & 255U; - - w = (g * inv256[(f >> 1) & 127]) & m; - - g += f * w; - q += u * w; - r += v * w; - - VERIFY_CHECK((g & m) == 0); - } - - t[0] = (int32_t)u; - t[1] = (int32_t)v; - t[2] = (int32_t)q; - t[3] = (int32_t)r; - - return eta; -} - -static void secp256k1_scalar_update_de_30(int32_t *d, int32_t *e, const int32_t *t) { - - /* I30 == -P^-1 mod 2^30 */ - const uint32_t I30 = 0x1588B13FL; - const int32_t P[9] = { 0x10364141L, 0x3F497A33L, 0x348A03BBL, 0x2BB739ABL, -0x146L, - 0, 0, 0, 65536 }; - const int32_t M30 = (int32_t)(UINT32_MAX >> 2); - int32_t u = t[0], v = t[1], q = t[2], r = t[3], di, ei, md, me; - int64_t cd, ce; - int i; - - di = d[0]; - ei = e[0]; - - cd = (int64_t)u * di + (int64_t)v * ei; - ce = (int64_t)q * di + (int64_t)r * ei; - - /* Calculate the multiples of P to add, to zero the 30 bottom bits. We choose md, me - * from the centred range [-2^29, 2^29) to keep d, e within [-2^256, 2^256). */ - md = ((int32_t)(I30 * 4 * (uint32_t)cd)) >> 2; - me = ((int32_t)(I30 * 4 * (uint32_t)ce)) >> 2; - - cd += (int64_t)P[0] * md; - ce += (int64_t)P[0] * me; - - VERIFY_CHECK(((int32_t)cd & M30) == 0); cd >>= 30; - VERIFY_CHECK(((int32_t)ce & M30) == 0); ce >>= 30; - - for (i = 1; i < 9; ++i) { - - di = d[i]; - ei = e[i]; - - cd += (int64_t)u * di + (int64_t)v * ei; - ce += (int64_t)q * di + (int64_t)r * ei; - - cd += (int64_t)P[i] * md; - ce += (int64_t)P[i] * me; - - d[i - 1] = (int32_t)cd & M30; cd >>= 30; - e[i - 1] = (int32_t)ce & M30; ce >>= 30; - } - - d[8] = (int32_t)cd; - e[8] = (int32_t)ce; -} - -static void secp256k1_scalar_update_fg_30(int32_t *f, int32_t *g, int32_t *t) { - - const int32_t M30 = (int32_t)(UINT32_MAX >> 2); - int32_t u = t[0], v = t[1], q = t[2], r = t[3], fi, gi; - int64_t cf, cg; - int i; - - fi = f[0]; - gi = g[0]; - - cf = (int64_t)u * fi + (int64_t)v * gi; - cg = (int64_t)q * fi + (int64_t)r * gi; - - VERIFY_CHECK(((int32_t)cf & M30) == 0); - VERIFY_CHECK(((int32_t)cg & M30) == 0); - - cf >>= 30; - cg >>= 30; - - for (i = 1; i < 9; ++i) { - - fi = f[i]; - gi = g[i]; - - cf += (int64_t)u * fi + (int64_t)v * gi; - cg += (int64_t)q * fi + (int64_t)r * gi; - - f[i - 1] = (int32_t)cf & M30; cf >>= 30; - g[i - 1] = (int32_t)cg & M30; cg >>= 30; - } - - f[8] = (int32_t)cf; - g[8] = (int32_t)cg; -} - -static void secp256k1_scalar_update_fg_30_var(int len, int32_t *f, int32_t *g, const int32_t *t) { - - const int32_t M30 = (int32_t)(UINT32_MAX >> 2); - const int32_t u = t[0], v = t[1], q = t[2], r = t[3]; - int32_t fi, gi; - int64_t cf, cg; - int i; - - VERIFY_CHECK(len > 0); - - fi = f[0]; - gi = g[0]; - - cf = (int64_t)u * fi + (int64_t)v * gi; - cg = (int64_t)q * fi + (int64_t)r * gi; - - VERIFY_CHECK(((int32_t)cf & M30) == 0); cf >>= 30; - VERIFY_CHECK(((int32_t)cg & M30) == 0); cg >>= 30; - - for (i = 1; i < len; ++i) { - - fi = f[i]; - gi = g[i]; - - cf += (int64_t)u * fi + (int64_t)v * gi; - cg += (int64_t)q * fi + (int64_t)r * gi; - - f[i - 1] = (int32_t)cf & M30; cf >>= 30; - g[i - 1] = (int32_t)cg & M30; cg >>= 30; - } - - f[len - 1] = (int32_t)cf; - g[len - 1] = (int32_t)cg; -} +static const secp256k1_modinv32_modinfo secp256k1_const_modinfo_scalar = { + {{0x10364141L, 0x3F497A33L, 0x348A03BBL, 0x2BB739ABL, -0x146L, 0, 0, 0, 65536}}, + 0x1588B13FL * 4 +}; static void secp256k1_scalar_inverse(secp256k1_scalar *r, const secp256k1_scalar *x) { - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int32_t t[4]; - int32_t d[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t e[9] = { 1, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t f[9] = { 0x10364141L, 0x3F497A33L, 0x348A03BBL, 0x2BB739ABL, - 0x3FFFFEBAL, 0x3FFFFFFFL, 0x3FFFFFFFL, 0x3FFFFFFFL, 0xFFFFL }; - int32_t g[9]; - secp256k1_scalar b0; - int i, sign; - uint32_t eta; + int sign; + secp256k1_modinv32_signed30 s; #ifdef VERIFY int zero_in = secp256k1_scalar_is_zero(x); #endif - - b0 = *x; - secp256k1_scalar_encode_30(g, &b0); - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If the maximum bitlength of g is known to be less than 256, then eta can be set - * initially to -(1 + (256 - maxlen(g))), and only (741 - (256 - maxlen(g))) total - * divsteps are needed. */ - eta = -(uint32_t)1; - - for (i = 0; i < 25; ++i) { - eta = secp256k1_scalar_divsteps_30(eta, f[0], g[0], t); - secp256k1_scalar_update_de_30(d, e, t); - secp256k1_scalar_update_fg_30(f, g, t); - } - - /* At this point sufficient iterations have been performed that g must have reached 0 - * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g - * values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - VERIFY_CHECK((g[0] | g[1] | g[2] | g[3] | g[4] | g[5] | g[6] | g[7] | g[8]) == 0); - - sign = (f[0] >> 1) & 1; - - secp256k1_scalar_decode_30(&b0, d); - secp256k1_scalar_cond_negate(&b0, sign); + secp256k1_scalar_to_signed30(&s, x); + secp256k1_modinv32(&sign, &s, &secp256k1_const_modinfo_scalar); + secp256k1_scalar_from_signed30(r, &s); + secp256k1_scalar_cond_negate(r, sign); #ifdef VERIFY - VERIFY_CHECK(!secp256k1_scalar_is_zero(&b0) == !zero_in); + VERIFY_CHECK(secp256k1_scalar_is_zero(r) == zero_in); #endif - - *r = b0; } static void secp256k1_scalar_inverse_var(secp256k1_scalar *r, const secp256k1_scalar *x) { - - /* Modular inversion based on the paper "Fast constant-time gcd computation and - * modular inversion" by Daniel J. Bernstein and Bo-Yin Yang. */ - - int32_t t[4]; - int32_t d[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t e[9] = { 1, 0, 0, 0, 0, 0, 0, 0, 0 }; - int32_t f[9] = { 0x10364141L, 0x3F497A33L, 0x348A03BBL, 0x2BB739ABL, - 0x3FFFFEBAL, 0x3FFFFFFFL, 0x3FFFFFFFL, 0x3FFFFFFFL, 0xFFFFL }; - int32_t g[9]; - secp256k1_scalar b; - int i, j, len = 9, sign; - uint32_t eta; - int32_t cond, fn, gn; + int sign; + secp256k1_modinv32_signed30 s; #ifdef VERIFY int zero_in = secp256k1_scalar_is_zero(x); #endif - - b = *x; - secp256k1_scalar_encode_30(g, &b); - - /* The paper uses 'delta'; eta == -delta (a performance tweak). - * - * If g has leading zeros (w.r.t 256 bits), then eta can be set initially to - * -(1 + clz(g)), and the worst-case divstep count would be only (741 - clz(g)). */ - eta = -(uint32_t)1; - - for (i = 0; i < 25; ++i) { - - eta = secp256k1_scalar_divsteps_30_var(eta, f[0], g[0], t); - secp256k1_scalar_update_de_30(d, e, t); - secp256k1_scalar_update_fg_30_var(len, f, g, t); - - if (g[0] == 0) { - cond = 0; - for (j = 1; j < len; ++j) { - cond |= g[j]; - } - if (cond == 0) { - break; - } - } - - fn = f[len - 1]; - gn = g[len - 1]; - - cond = ((int32_t)len - 2) >> 31; - cond |= fn ^ (fn >> 31); - cond |= gn ^ (gn >> 31); - - if (cond == 0) { - f[len - 2] |= (uint32_t)fn << 30; - g[len - 2] |= (uint32_t)gn << 30; - --len; - } - } - - VERIFY_CHECK(i < 25); - - /* At this point g is 0 and (if g was not originally 0) f must now equal +/- GCD of - * the initial f, g values i.e. +/- 1, and d now contains +/- the modular inverse. */ - - sign = (f[0] >> 1) & 1; - - secp256k1_scalar_decode_30(&b, d); - - if (sign) { - secp256k1_scalar_negate(&b, &b); - } + secp256k1_scalar_to_signed30(&s, x); + secp256k1_modinv32_var(&sign, &s, &secp256k1_const_modinfo_scalar); + secp256k1_scalar_from_signed30(r, &s); + if (sign) secp256k1_scalar_negate(r, r); #ifdef VERIFY - VERIFY_CHECK(!secp256k1_scalar_is_zero(&b) == !zero_in); + VERIFY_CHECK(secp256k1_scalar_is_zero(r) == zero_in); #endif - - *r = b; } SECP256K1_INLINE static int secp256k1_scalar_is_even(const secp256k1_scalar *a) {