src/scalar_4x64_impl.h

   1 /**********************************************************************
   2  * Copyright (c) 2013, 2014 Pieter Wuille                             *
   3  * Distributed under the MIT software license, see the accompanying   *
   4  * file COPYING or http://www.opensource.org/licenses/mit-license.php.*
   5  **********************************************************************/
   6
   7 #ifndef _SECP256K1_SCALAR_REPR_IMPL_H_
   8 #define _SECP256K1_SCALAR_REPR_IMPL_H_
   9
  10 typedef unsigned __int128 uint128_t;
  11
  12 /* Limbs of the secp256k1 order. */
  13 #define SECP256K1_N_0 ((uint64_t)0xBFD25E8CD0364141ULL)
  14 #define SECP256K1_N_1 ((uint64_t)0xBAAEDCE6AF48A03BULL)
  15 #define SECP256K1_N_2 ((uint64_t)0xFFFFFFFFFFFFFFFEULL)
  16 #define SECP256K1_N_3 ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
  17
  18 /* Limbs of 2^256 minus the secp256k1 order. */
  19 #define SECP256K1_N_C_0 (~SECP256K1_N_0 + 1)
  20 #define SECP256K1_N_C_1 (~SECP256K1_N_1)
  21 #define SECP256K1_N_C_2 (1)
  22
  23 /* Limbs of half the secp256k1 order. */
  24 #define SECP256K1_N_H_0 ((uint64_t)0xDFE92F46681B20A0ULL)
  25 #define SECP256K1_N_H_1 ((uint64_t)0x5D576E7357A4501DULL)
  26 #define SECP256K1_N_H_2 ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
  27 #define SECP256K1_N_H_3 ((uint64_t)0x7FFFFFFFFFFFFFFFULL)
  28
  29 SECP256K1_INLINE static void secp256k1_scalar_clear(secp256k1_scalar_t *r) {
  30     r->d[0] = 0;
  31     r->d[1] = 0;
  32     r->d[2] = 0;
  33     r->d[3] = 0;
  34 }
  35
  36 SECP256K1_INLINE static void secp256k1_scalar_set_int(secp256k1_scalar_t *r, unsigned int v) {
  37     r->d[0] = v;
  38     r->d[1] = 0;
  39     r->d[2] = 0;
  40     r->d[3] = 0;
  41 }
  42
  43 SECP256K1_INLINE static unsigned int secp256k1_scalar_get_bits(const secp256k1_scalar_t *a, unsigned int offset, unsigned int count) {
  44     VERIFY_CHECK((offset + count - 1) >> 6 == offset >> 6);
  45     return (a->d[offset >> 6] >> (offset & 0x3F)) & ((((uint64_t)1) << count) - 1);
  46 }
  47
  48 SECP256K1_INLINE static unsigned int secp256k1_scalar_get_bits_var(const secp256k1_scalar_t *a, unsigned int offset, unsigned int count) {
  49     VERIFY_CHECK(count < 32);
  50     VERIFY_CHECK(offset + count <= 256);
  51     if ((offset + count - 1) >> 6 == offset >> 6) {
  52         return secp256k1_scalar_get_bits(a, offset, count);
  53     } else {
  54         VERIFY_CHECK((offset >> 6) + 1 < 4);
  55         return ((a->d[offset >> 6] >> (offset & 0x3F)) | (a->d[(offset >> 6) + 1] << (64 - (offset & 0x3F)))) & ((((uint64_t)1) << count) - 1);
  56     }
  57 }
  58
  59 SECP256K1_INLINE static int secp256k1_scalar_check_overflow(const secp256k1_scalar_t *a) {
  60     int yes = 0;
  61     int no = 0;
  62     no |= (a->d[3] < SECP256K1_N_3); /* No need for a > check. */
  63     no |= (a->d[2] < SECP256K1_N_2);
  64     yes |= (a->d[2] > SECP256K1_N_2) & ~no;
  65     no |= (a->d[1] < SECP256K1_N_1);
  66     yes |= (a->d[1] > SECP256K1_N_1) & ~no;
  67     yes |= (a->d[0] >= SECP256K1_N_0) & ~no;
  68     return yes;
  69 }
  70
  71 SECP256K1_INLINE static int secp256k1_scalar_reduce(secp256k1_scalar_t *r, unsigned int overflow) {
  72     VERIFY_CHECK(overflow <= 1);
  73     uint128_t t = (uint128_t)r->d[0] + overflow * SECP256K1_N_C_0;
  74     r->d[0] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  75     t += (uint128_t)r->d[1] + overflow * SECP256K1_N_C_1;
  76     r->d[1] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  77     t += (uint128_t)r->d[2] + overflow * SECP256K1_N_C_2;
  78     r->d[2] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  79     t += (uint64_t)r->d[3];
  80     r->d[3] = t & 0xFFFFFFFFFFFFFFFFULL;
  81     return overflow;
  82 }
  83
  84 static void secp256k1_scalar_add(secp256k1_scalar_t *r, const secp256k1_scalar_t *a, const secp256k1_scalar_t *b) {
  85     uint128_t t = (uint128_t)a->d[0] + b->d[0];
  86     r->d[0] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  87     t += (uint128_t)a->d[1] + b->d[1];
  88     r->d[1] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  89     t += (uint128_t)a->d[2] + b->d[2];
  90     r->d[2] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  91     t += (uint128_t)a->d[3] + b->d[3];
  92     r->d[3] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
  93     secp256k1_scalar_reduce(r, t + secp256k1_scalar_check_overflow(r));
  94 }
  95
  96 static void secp256k1_scalar_add_bit(secp256k1_scalar_t *r, unsigned int bit) {
  97     VERIFY_CHECK(bit < 256);
  98     uint128_t t = (uint128_t)r->d[0] + (((uint64_t)((bit >> 6) == 0)) << bit);
  99     r->d[0] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
 100     t += (uint128_t)r->d[1] + (((uint64_t)((bit >> 6) == 1)) << (bit & 0x3F));
 101     r->d[1] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
 102     t += (uint128_t)r->d[2] + (((uint64_t)((bit >> 6) == 2)) << (bit & 0x3F));
 103     r->d[2] = t & 0xFFFFFFFFFFFFFFFFULL; t >>= 64;
 104     t += (uint128_t)r->d[3] + (((uint64_t)((bit >> 6) == 3)) << (bit & 0x3F));
 105     r->d[3] = t & 0xFFFFFFFFFFFFFFFFULL;
 106 #ifdef VERIFY
 107     VERIFY_CHECK((t >> 64) == 0);
 108     VERIFY_CHECK(secp256k1_scalar_check_overflow(r) == 0);
 109 #endif
 110 }
 111
 112 static void secp256k1_scalar_set_b32(secp256k1_scalar_t *r, const unsigned char *b32, int *overflow) {
 113     r->d[0] = (uint64_t)b32[31] | (uint64_t)b32[30] << 8 | (uint64_t)b32[29] << 16 | (uint64_t)b32[28] << 24 | (uint64_t)b32[27] << 32 | (uint64_t)b32[26] << 40 | (uint64_t)b32[25] << 48 | (uint64_t)b32[24] << 56;
 114     r->d[1] = (uint64_t)b32[23] | (uint64_t)b32[22] << 8 | (uint64_t)b32[21] << 16 | (uint64_t)b32[20] << 24 | (uint64_t)b32[19] << 32 | (uint64_t)b32[18] << 40 | (uint64_t)b32[17] << 48 | (uint64_t)b32[16] << 56;
 115     r->d[2] = (uint64_t)b32[15] | (uint64_t)b32[14] << 8 | (uint64_t)b32[13] << 16 | (uint64_t)b32[12] << 24 | (uint64_t)b32[11] << 32 | (uint64_t)b32[10] << 40 | (uint64_t)b32[9] << 48 | (uint64_t)b32[8] << 56;
 116     r->d[3] = (uint64_t)b32[7] | (uint64_t)b32[6] << 8 | (uint64_t)b32[5] << 16 | (uint64_t)b32[4] << 24 | (uint64_t)b32[3] << 32 | (uint64_t)b32[2] << 40 | (uint64_t)b32[1] << 48 | (uint64_t)b32[0] << 56;
 117     int over = secp256k1_scalar_reduce(r, secp256k1_scalar_check_overflow(r));
 118     if (overflow) {
 119         *overflow = over;
 120     }
 121 }
 122
 123 static void secp256k1_scalar_get_b32(unsigned char *bin, const secp256k1_scalar_t* a) {
 124     bin[0] = a->d[3] >> 56; bin[1] = a->d[3] >> 48; bin[2] = a->d[3] >> 40; bin[3] = a->d[3] >> 32; bin[4] = a->d[3] >> 24; bin[5] = a->d[3] >> 16; bin[6] = a->d[3] >> 8; bin[7] = a->d[3];
 125     bin[8] = a->d[2] >> 56; bin[9] = a->d[2] >> 48; bin[10] = a->d[2] >> 40; bin[11] = a->d[2] >> 32; bin[12] = a->d[2] >> 24; bin[13] = a->d[2] >> 16; bin[14] = a->d[2] >> 8; bin[15] = a->d[2];
 126     bin[16] = a->d[1] >> 56; bin[17] = a->d[1] >> 48; bin[18] = a->d[1] >> 40; bin[19] = a->d[1] >> 32; bin[20] = a->d[1] >> 24; bin[21] = a->d[1] >> 16; bin[22] = a->d[1] >> 8; bin[23] = a->d[1];
 127     bin[24] = a->d[0] >> 56; bin[25] = a->d[0] >> 48; bin[26] = a->d[0] >> 40; bin[27] = a->d[0] >> 32; bin[28] = a->d[0] >> 24; bin[29] = a->d[0] >> 16; bin[30] = a->d[0] >> 8; bin[31] = a->d[0];
 128 }
 129
 130 SECP256K1_INLINE static int secp256k1_scalar_is_zero(const secp256k1_scalar_t *a) {
 131     return (a->d[0] | a->d[1] | a->d[2] | a->d[3]) == 0;
 132 }
 133
 134 static void secp256k1_scalar_negate(secp256k1_scalar_t *r, const secp256k1_scalar_t *a) {
 135     uint64_t nonzero = 0xFFFFFFFFFFFFFFFFULL * (secp256k1_scalar_is_zero(a) == 0);
 136     uint128_t t = (uint128_t)(~a->d[0]) + SECP256K1_N_0 + 1;
 137     r->d[0] = t & nonzero; t >>= 64;
 138     t += (uint128_t)(~a->d[1]) + SECP256K1_N_1;
 139     r->d[1] = t & nonzero; t >>= 64;
 140     t += (uint128_t)(~a->d[2]) + SECP256K1_N_2;
 141     r->d[2] = t & nonzero; t >>= 64;
 142     t += (uint128_t)(~a->d[3]) + SECP256K1_N_3;
 143     r->d[3] = t & nonzero;
 144 }
 145
 146 SECP256K1_INLINE static int secp256k1_scalar_is_one(const secp256k1_scalar_t *a) {
 147     return ((a->d[0] ^ 1) | a->d[1] | a->d[2] | a->d[3]) == 0;
 148 }
 149
 150 static int secp256k1_scalar_is_high(const secp256k1_scalar_t *a) {
 151     int yes = 0;
 152     int no = 0;
 153     no |= (a->d[3] < SECP256K1_N_H_3);
 154     yes |= (a->d[3] > SECP256K1_N_H_3) & ~no;
 155     no |= (a->d[2] < SECP256K1_N_H_2) & ~yes; /* No need for a > check. */
 156     no |= (a->d[1] < SECP256K1_N_H_1) & ~yes;
 157     yes |= (a->d[1] > SECP256K1_N_H_1) & ~no;
 158     yes |= (a->d[0] > SECP256K1_N_H_0) & ~no;
 159     return yes;
 160 }
 161
 162 /* Inspired by the macros in OpenSSL's crypto/bn/asm/x86_64-gcc.c. */
 163
 164 /** Add a*b to the number defined by (c0,c1,c2). c2 must never overflow. */
 165 #define muladd(a,b) { \
 166     uint64_t tl, th; \
 167     { \
 168         uint128_t t = (uint128_t)a * b; \
 169         th = t >> 64;         /* at most 0xFFFFFFFFFFFFFFFE */ \
 170         tl = t; \
 171     } \
 172     c0 += tl;                 /* overflow is handled on the next line */ \
 173     th += (c0 < tl) ? 1 : 0;  /* at most 0xFFFFFFFFFFFFFFFF */ \
 174     c1 += th;                 /* overflow is handled on the next line */ \
 175     c2 += (c1 < th) ? 1 : 0;  /* never overflows by contract (verified in the next line) */ \
 176     VERIFY_CHECK((c1 >= th) || (c2 != 0)); \
 177 }
 178
 179 /** Add a*b to the number defined by (c0,c1). c1 must never overflow. */
 180 #define muladd_fast(a,b) { \
 181     uint64_t tl, th; \
 182     { \
 183         uint128_t t = (uint128_t)a * b; \
 184         th = t >> 64;         /* at most 0xFFFFFFFFFFFFFFFE */ \
 185         tl = t; \
 186     } \
 187     c0 += tl;                 /* overflow is handled on the next line */ \
 188     th += (c0 < tl) ? 1 : 0;  /* at most 0xFFFFFFFFFFFFFFFF */ \
 189     c1 += th;                 /* never overflows by contract (verified in the next line) */ \
 190     VERIFY_CHECK(c1 >= th); \
 191 }
 192
 193 /** Add 2*a*b to the number defined by (c0,c1,c2). c2 must never overflow. */
 194 #define muladd2(a,b) { \
 195     uint64_t tl, th; \
 196     { \
 197         uint128_t t = (uint128_t)a * b; \
 198         th = t >> 64;               /* at most 0xFFFFFFFFFFFFFFFE */ \
 199         tl = t; \
 200     } \
 201     uint64_t th2 = th + th;         /* at most 0xFFFFFFFFFFFFFFFE (in case th was 0x7FFFFFFFFFFFFFFF) */ \
 202     c2 += (th2 < th) ? 1 : 0;       /* never overflows by contract (verified the next line) */ \
 203     VERIFY_CHECK((th2 >= th) || (c2 != 0)); \
 204     uint64_t tl2 = tl + tl;         /* at most 0xFFFFFFFFFFFFFFFE (in case the lowest 63 bits of tl were 0x7FFFFFFFFFFFFFFF) */ \
 205     th2 += (tl2 < tl) ? 1 : 0;      /* at most 0xFFFFFFFFFFFFFFFF */ \
 206     c0 += tl2;                      /* overflow is handled on the next line */ \
 207     th2 += (c0 < tl2) ? 1 : 0;      /* second overflow is handled on the next line */ \
 208     c2 += (c0 < tl2) & (th2 == 0);  /* never overflows by contract (verified the next line) */ \
 209     VERIFY_CHECK((c0 >= tl2) || (th2 != 0) || (c2 != 0)); \
 210     c1 += th2;                      /* overflow is handled on the next line */ \
 211     c2 += (c1 < th2) ? 1 : 0;       /* never overflows by contract (verified the next line) */ \
 212     VERIFY_CHECK((c1 >= th2) || (c2 != 0)); \
 213 }
 214
 215 /** Add a to the number defined by (c0,c1,c2). c2 must never overflow. */
 216 #define sumadd(a) { \
 217     c0 += (a);                  /* overflow is handled on the next line */ \
 218     unsigned int over = (c0 < (a)) ? 1 : 0; \
 219     c1 += over;                 /* overflow is handled on the next line */ \
 220     c2 += (c1 < over) ? 1 : 0;  /* never overflows by contract */ \
 221 }
 222
 223 /** Add a to the number defined by (c0,c1). c1 must never overflow, c2 must be zero. */
 224 #define sumadd_fast(a) { \
 225     c0 += (a);                 /* overflow is handled on the next line */ \
 226     c1 += (c0 < (a)) ? 1 : 0;  /* never overflows by contract (verified the next line) */ \
 227     VERIFY_CHECK((c1 != 0) | (c0 >= (a))); \
 228     VERIFY_CHECK(c2 == 0); \
 229 }
 230
 231 /** Extract the lowest 64 bits of (c0,c1,c2) into n, and left shift the number 64 bits. */
 232 #define extract(n) { \
 233     (n) = c0; \
 234     c0 = c1; \
 235     c1 = c2; \
 236     c2 = 0; \
 237 }
 238
 239 /** Extract the lowest 64 bits of (c0,c1,c2) into n, and left shift the number 64 bits. c2 is required to be zero. */
 240 #define extract_fast(n) { \
 241     (n) = c0; \
 242     c0 = c1; \
 243     c1 = 0; \
 244     VERIFY_CHECK(c2 == 0); \
 245 }
 246
 247 static void secp256k1_scalar_reduce_512(secp256k1_scalar_t *r, const uint64_t *l) {
 248     uint64_t n0 = l[4], n1 = l[5], n2 = l[6], n3 = l[7];
 249
 250     /* 160 bit accumulator. */
 251     uint64_t c0, c1;
 252     uint32_t c2;
 253
 254     /* Reduce 512 bits into 385. */
 255     /* m[0..6] = l[0..3] + n[0..3] * SECP256K1_N_C. */
 256     c0 = l[0]; c1 = 0; c2 = 0;
 257     muladd_fast(n0, SECP256K1_N_C_0);
 258     uint64_t m0; extract_fast(m0);
 259     sumadd_fast(l[1]);
 260     muladd(n1, SECP256K1_N_C_0);
 261     muladd(n0, SECP256K1_N_C_1);
 262     uint64_t m1; extract(m1);
 263     sumadd(l[2]);
 264     muladd(n2, SECP256K1_N_C_0);
 265     muladd(n1, SECP256K1_N_C_1);
 266     sumadd(n0);
 267     uint64_t m2; extract(m2);
 268     sumadd(l[3]);
 269     muladd(n3, SECP256K1_N_C_0);
 270     muladd(n2, SECP256K1_N_C_1);
 271     sumadd(n1);
 272     uint64_t m3; extract(m3);
 273     muladd(n3, SECP256K1_N_C_1);
 274     sumadd(n2);
 275     uint64_t m4; extract(m4);
 276     sumadd_fast(n3);
 277     uint64_t m5; extract_fast(m5);
 278     VERIFY_CHECK(c0 <= 1);
 279     uint32_t m6 = c0;
 280
 281     /* Reduce 385 bits into 258. */
 282     /* p[0..4] = m[0..3] + m[4..6] * SECP256K1_N_C. */
 283     c0 = m0; c1 = 0; c2 = 0;
 284     muladd_fast(m4, SECP256K1_N_C_0);
 285     uint64_t p0; extract_fast(p0);
 286     sumadd_fast(m1);
 287     muladd(m5, SECP256K1_N_C_0);
 288     muladd(m4, SECP256K1_N_C_1);
 289     uint64_t p1; extract(p1);
 290     sumadd(m2);
 291     muladd(m6, SECP256K1_N_C_0);
 292     muladd(m5, SECP256K1_N_C_1);
 293     sumadd(m4);
 294     uint64_t p2; extract(p2);
 295     sumadd_fast(m3);
 296     muladd_fast(m6, SECP256K1_N_C_1);
 297     sumadd_fast(m5);
 298     uint64_t p3; extract_fast(p3);
 299     uint32_t p4 = c0 + m6;
 300     VERIFY_CHECK(p4 <= 2);
 301
 302     /* Reduce 258 bits into 256. */
 303     /* r[0..3] = p[0..3] + p[4] * SECP256K1_N_C. */
 304     uint128_t c = p0 + (uint128_t)SECP256K1_N_C_0 * p4;
 305     r->d[0] = c & 0xFFFFFFFFFFFFFFFFULL; c >>= 64;
 306     c += p1 + (uint128_t)SECP256K1_N_C_1 * p4;
 307     r->d[1] = c & 0xFFFFFFFFFFFFFFFFULL; c >>= 64;
 308     c += p2 + (uint128_t)p4;
 309     r->d[2] = c & 0xFFFFFFFFFFFFFFFFULL; c >>= 64;
 310     c += p3;
 311     r->d[3] = c & 0xFFFFFFFFFFFFFFFFULL; c >>= 64;
 312
 313     /* Final reduction of r. */
 314     secp256k1_scalar_reduce(r, c + secp256k1_scalar_check_overflow(r));
 315 }
 316
 317 static void secp256k1_scalar_mul(secp256k1_scalar_t *r, const secp256k1_scalar_t *a, const secp256k1_scalar_t *b) {
 318     /* 160 bit accumulator. */
 319     uint64_t c0 = 0, c1 = 0;
 320     uint32_t c2 = 0;
 321
 322     uint64_t l[8];
 323
 324     /* l[0..7] = a[0..3] * b[0..3]. */
 325     muladd_fast(a->d[0], b->d[0]);
 326     extract_fast(l[0]);
 327     muladd(a->d[0], b->d[1]);
 328     muladd(a->d[1], b->d[0]);
 329     extract(l[1]);
 330     muladd(a->d[0], b->d[2]);
 331     muladd(a->d[1], b->d[1]);
 332     muladd(a->d[2], b->d[0]);
 333     extract(l[2]);
 334     muladd(a->d[0], b->d[3]);
 335     muladd(a->d[1], b->d[2]);
 336     muladd(a->d[2], b->d[1]);
 337     muladd(a->d[3], b->d[0]);
 338     extract(l[3]);
 339     muladd(a->d[1], b->d[3]);
 340     muladd(a->d[2], b->d[2]);
 341     muladd(a->d[3], b->d[1]);
 342     extract(l[4]);
 343     muladd(a->d[2], b->d[3]);
 344     muladd(a->d[3], b->d[2]);
 345     extract(l[5]);
 346     muladd_fast(a->d[3], b->d[3]);
 347     extract_fast(l[6]);
 348     VERIFY_CHECK(c1 <= 0);
 349     l[7] = c0;
 350
 351     secp256k1_scalar_reduce_512(r, l);
 352 }
 353
 354 static void secp256k1_scalar_sqr(secp256k1_scalar_t *r, const secp256k1_scalar_t *a) {
 355     /* 160 bit accumulator. */
 356     uint64_t c0 = 0, c1 = 0;
 357     uint32_t c2 = 0;
 358
 359     uint64_t l[8];
 360
 361     /* l[0..7] = a[0..3] * b[0..3]. */
 362     muladd_fast(a->d[0], a->d[0]);
 363     extract_fast(l[0]);
 364     muladd2(a->d[0], a->d[1]);
 365     extract(l[1]);
 366     muladd2(a->d[0], a->d[2]);
 367     muladd(a->d[1], a->d[1]);
 368     extract(l[2]);
 369     muladd2(a->d[0], a->d[3]);
 370     muladd2(a->d[1], a->d[2]);
 371     extract(l[3]);
 372     muladd2(a->d[1], a->d[3]);
 373     muladd(a->d[2], a->d[2]);
 374     extract(l[4]);
 375     muladd2(a->d[2], a->d[3]);
 376     extract(l[5]);
 377     muladd_fast(a->d[3], a->d[3]);
 378     extract_fast(l[6]);
 379     VERIFY_CHECK(c1 == 0);
 380     l[7] = c0;
 381
 382     secp256k1_scalar_reduce_512(r, l);
 383 }
 384
 385 #undef sumadd
 386 #undef sumadd_fast
 387 #undef muladd
 388 #undef muladd_fast
 389 #undef muladd2
 390 #undef extract
 391 #undef extract_fast
 392
 393 #endif