ec_c25519_m62.c
1 /* 2 * Copyright (c) 2018 Thomas Pornin <pornin@bolet.org> 3 * 4 * Permission is hereby granted, free of charge, to any person obtaining 5 * a copy of this software and associated documentation files (the 6 * "Software"), to deal in the Software without restriction, including 7 * without limitation the rights to use, copy, modify, merge, publish, 8 * distribute, sublicense, and/or sell copies of the Software, and to 9 * permit persons to whom the Software is furnished to do so, subject to 10 * the following conditions: 11 * 12 * The above copyright notice and this permission notice shall be 13 * included in all copies or substantial portions of the Software. 14 * 15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 16 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 17 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 18 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS 19 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN 20 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN 21 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 22 * SOFTWARE. 23 */ 24 25 #include "inner.h" 26 27 #if BR_INT128 || BR_UMUL128 28 29 #if BR_UMUL128 30 #include <intrin.h> 31 #endif 32 33 static const unsigned char GEN[] = { 34 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 35 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 36 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 37 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00 38 }; 39 40 static const unsigned char ORDER[] = { 41 0x7F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 42 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 43 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 44 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF 45 }; 46 47 static const unsigned char * 48 api_generator(int curve, size_t *len) 49 { 50 (void)curve; 51 *len = 32; 52 return GEN; 53 } 54 55 static const unsigned char * 56 api_order(int curve, size_t *len) 57 { 58 (void)curve; 59 *len = 32; 60 return ORDER; 61 } 62 63 static size_t 64 api_xoff(int curve, size_t *len) 65 { 66 (void)curve; 67 *len = 32; 68 return 0; 69 } 70 71 /* 72 * A field element is encoded as five 64-bit integers, in basis 2^51. 73 * Limbs may be occasionally larger than 2^51, to save on carry 74 * propagation costs. 75 */ 76 77 #define MASK51 (((uint64_t)1 << 51) - (uint64_t)1) 78 79 /* 80 * Swap two field elements, conditionally on a flag. 81 */ 82 static inline void 83 f255_cswap(uint64_t *a, uint64_t *b, uint32_t ctl) 84 { 85 uint64_t m, w; 86 87 m = -(uint64_t)ctl; 88 w = m & (a[0] ^ b[0]); a[0] ^= w; b[0] ^= w; 89 w = m & (a[1] ^ b[1]); a[1] ^= w; b[1] ^= w; 90 w = m & (a[2] ^ b[2]); a[2] ^= w; b[2] ^= w; 91 w = m & (a[3] ^ b[3]); a[3] ^= w; b[3] ^= w; 92 w = m & (a[4] ^ b[4]); a[4] ^= w; b[4] ^= w; 93 } 94 95 /* 96 * Addition with no carry propagation. Limbs double in size. 97 */ 98 static inline void 99 f255_add(uint64_t *d, const uint64_t *a, const uint64_t *b) 100 { 101 d[0] = a[0] + b[0]; 102 d[1] = a[1] + b[1]; 103 d[2] = a[2] + b[2]; 104 d[3] = a[3] + b[3]; 105 d[4] = a[4] + b[4]; 106 } 107 108 /* 109 * Subtraction. 110 * On input, limbs must fit on 60 bits each. On output, result is 111 * partially reduced, with max value 2^255+19456; moreover, all 112 * limbs will fit on 51 bits, except the low limb, which may have 113 * value up to 2^51+19455. 114 */ 115 static inline void 116 f255_sub(uint64_t *d, const uint64_t *a, const uint64_t *b) 117 { 118 uint64_t cc, w; 119 120 /* 121 * We compute d = (2^255-19)*1024 + a - b. Since the limbs 122 * fit on 60 bits, the maximum value of operands are slightly 123 * more than 2^264, but much less than 2^265-19456. This 124 * ensures that the result is positive. 125 */ 126 127 /* 128 * Initial carry is 19456, since we add 2^265-19456. Each 129 * individual subtraction may yield a carry up to 513. 130 */ 131 w = a[0] - b[0] - 19456; 132 d[0] = w & MASK51; 133 cc = -(w >> 51) & 0x3FF; 134 w = a[1] - b[1] - cc; 135 d[1] = w & MASK51; 136 cc = -(w >> 51) & 0x3FF; 137 w = a[2] - b[2] - cc; 138 d[2] = w & MASK51; 139 cc = -(w >> 51) & 0x3FF; 140 w = a[3] - b[3] - cc; 141 d[3] = w & MASK51; 142 cc = -(w >> 51) & 0x3FF; 143 d[4] = ((uint64_t)1 << 61) + a[4] - b[4] - cc; 144 145 /* 146 * Partial reduction. The intermediate result may be up to 147 * slightly above 2^265, but less than 2^265+2^255. When we 148 * truncate to 255 bits, the upper bits will be at most 1024. 149 */ 150 d[0] += 19 * (d[4] >> 51); 151 d[4] &= MASK51; 152 } 153 154 /* 155 * UMUL51(hi, lo, x, y) computes: 156 * 157 * hi = floor((x * y) / (2^51)) 158 * lo = x * y mod 2^51 159 * 160 * Note that lo < 2^51, but "hi" may be larger, if the input operands are 161 * larger. 162 */ 163 #if BR_INT128 164 165 #define UMUL51(hi, lo, x, y) do { \ 166 unsigned __int128 umul_tmp; \ 167 umul_tmp = (unsigned __int128)(x) * (unsigned __int128)(y); \ 168 (hi) = (uint64_t)(umul_tmp >> 51); \ 169 (lo) = (uint64_t)umul_tmp & MASK51; \ 170 } while (0) 171 172 #elif BR_UMUL128 173 174 #define UMUL51(hi, lo, x, y) do { \ 175 uint64_t umul_hi, umul_lo; \ 176 umul_lo = _umul128((x), (y), &umul_hi); \ 177 (hi) = (umul_hi << 13) | (umul_lo >> 51); \ 178 (lo) = umul_lo & MASK51; \ 179 } while (0) 180 181 #endif 182 183 /* 184 * Multiplication. 185 * On input, limbs must fit on 54 bits each. 186 * On output, limb 0 is at most 2^51 + 155647, and other limbs fit 187 * on 51 bits each. 188 */ 189 static inline void 190 f255_mul(uint64_t *d, uint64_t *a, uint64_t *b) 191 { 192 uint64_t t[10], hi, lo, w, cc; 193 194 /* 195 * Perform cross products, accumulating values without carry 196 * propagation. 197 * 198 * Since input limbs fit on 54 bits each, each individual 199 * UMUL51 will produce a "hi" of less than 2^57. The maximum 200 * sum will be at most 5*(2^57-1) + 4*(2^51-1) (for t[5]), 201 * i.e. less than 324*2^51. 202 */ 203 204 UMUL51(t[1], t[0], a[0], b[0]); 205 206 UMUL51(t[2], lo, a[1], b[0]); t[1] += lo; 207 UMUL51(hi, lo, a[0], b[1]); t[1] += lo; t[2] += hi; 208 209 UMUL51(t[3], lo, a[2], b[0]); t[2] += lo; 210 UMUL51(hi, lo, a[1], b[1]); t[2] += lo; t[3] += hi; 211 UMUL51(hi, lo, a[0], b[2]); t[2] += lo; t[3] += hi; 212 213 UMUL51(t[4], lo, a[3], b[0]); t[3] += lo; 214 UMUL51(hi, lo, a[2], b[1]); t[3] += lo; t[4] += hi; 215 UMUL51(hi, lo, a[1], b[2]); t[3] += lo; t[4] += hi; 216 UMUL51(hi, lo, a[0], b[3]); t[3] += lo; t[4] += hi; 217 218 UMUL51(t[5], lo, a[4], b[0]); t[4] += lo; 219 UMUL51(hi, lo, a[3], b[1]); t[4] += lo; t[5] += hi; 220 UMUL51(hi, lo, a[2], b[2]); t[4] += lo; t[5] += hi; 221 UMUL51(hi, lo, a[1], b[3]); t[4] += lo; t[5] += hi; 222 UMUL51(hi, lo, a[0], b[4]); t[4] += lo; t[5] += hi; 223 224 UMUL51(t[6], lo, a[4], b[1]); t[5] += lo; 225 UMUL51(hi, lo, a[3], b[2]); t[5] += lo; t[6] += hi; 226 UMUL51(hi, lo, a[2], b[3]); t[5] += lo; t[6] += hi; 227 UMUL51(hi, lo, a[1], b[4]); t[5] += lo; t[6] += hi; 228 229 UMUL51(t[7], lo, a[4], b[2]); t[6] += lo; 230 UMUL51(hi, lo, a[3], b[3]); t[6] += lo; t[7] += hi; 231 UMUL51(hi, lo, a[2], b[4]); t[6] += lo; t[7] += hi; 232 233 UMUL51(t[8], lo, a[4], b[3]); t[7] += lo; 234 UMUL51(hi, lo, a[3], b[4]); t[7] += lo; t[8] += hi; 235 236 UMUL51(t[9], lo, a[4], b[4]); t[8] += lo; 237 238 /* 239 * The upper words t[5]..t[9] are folded back into the lower 240 * words, using the rule that 2^255 = 19 in the field. 241 * 242 * Since each t[i] is less than 324*2^51, the additions below 243 * will yield less than 6480*2^51 in each limb; this fits in 244 * 64 bits (6480*2^51 < 8192*2^51 = 2^64), hence there is 245 * no overflow. 246 */ 247 t[0] += 19 * t[5]; 248 t[1] += 19 * t[6]; 249 t[2] += 19 * t[7]; 250 t[3] += 19 * t[8]; 251 t[4] += 19 * t[9]; 252 253 /* 254 * Propagate carries. 255 */ 256 w = t[0]; 257 d[0] = w & MASK51; 258 cc = w >> 51; 259 w = t[1] + cc; 260 d[1] = w & MASK51; 261 cc = w >> 51; 262 w = t[2] + cc; 263 d[2] = w & MASK51; 264 cc = w >> 51; 265 w = t[3] + cc; 266 d[3] = w & MASK51; 267 cc = w >> 51; 268 w = t[4] + cc; 269 d[4] = w & MASK51; 270 cc = w >> 51; 271 272 /* 273 * Since the limbs were 64-bit values, the top carry is at 274 * most 8192 (in practice, that cannot be reached). We simply 275 * performed a partial reduction. 276 */ 277 d[0] += 19 * cc; 278 } 279 280 /* 281 * Multiplication by A24 = 121665. 282 * Input must have limbs of 60 bits at most. 283 */ 284 static inline void 285 f255_mul_a24(uint64_t *d, const uint64_t *a) 286 { 287 uint64_t t[5], cc, w; 288 289 /* 290 * 121665 = 15 * 8111. We first multiply by 15, with carry 291 * propagation and partial reduction. 292 */ 293 w = a[0] * 15; 294 t[0] = w & MASK51; 295 cc = w >> 51; 296 w = a[1] * 15 + cc; 297 t[1] = w & MASK51; 298 cc = w >> 51; 299 w = a[2] * 15 + cc; 300 t[2] = w & MASK51; 301 cc = w >> 51; 302 w = a[3] * 15 + cc; 303 t[3] = w & MASK51; 304 cc = w >> 51; 305 w = a[4] * 15 + cc; 306 t[4] = w & MASK51; 307 t[0] += 19 * (w >> 51); 308 309 /* 310 * Then multiplication by 8111. At that point, we known that 311 * t[0] is less than 2^51 + 19*8192, and other limbs are less 312 * than 2^51; thus, there will be no overflow. 313 */ 314 w = t[0] * 8111; 315 d[0] = w & MASK51; 316 cc = w >> 51; 317 w = t[1] * 8111 + cc; 318 d[1] = w & MASK51; 319 cc = w >> 51; 320 w = t[2] * 8111 + cc; 321 d[2] = w & MASK51; 322 cc = w >> 51; 323 w = t[3] * 8111 + cc; 324 d[3] = w & MASK51; 325 cc = w >> 51; 326 w = t[4] * 8111 + cc; 327 d[4] = w & MASK51; 328 d[0] += 19 * (w >> 51); 329 } 330 331 /* 332 * Finalize reduction. 333 * On input, limbs must fit on 51 bits, except possibly the low limb, 334 * which may be slightly above 2^51. 335 */ 336 static inline void 337 f255_final_reduce(uint64_t *a) 338 { 339 uint64_t t[5], cc, w; 340 341 /* 342 * We add 19. If the result (in t[]) is below 2^255, then a[] 343 * is already less than 2^255-19, thus already reduced. 344 * Otherwise, we subtract 2^255 from t[], in which case we 345 * have t = a - (2^255-19), and that's our result. 346 */ 347 w = a[0] + 19; 348 t[0] = w & MASK51; 349 cc = w >> 51; 350 w = a[1] + cc; 351 t[1] = w & MASK51; 352 cc = w >> 51; 353 w = a[2] + cc; 354 t[2] = w & MASK51; 355 cc = w >> 51; 356 w = a[3] + cc; 357 t[3] = w & MASK51; 358 cc = w >> 51; 359 w = a[4] + cc; 360 t[4] = w & MASK51; 361 cc = w >> 51; 362 363 /* 364 * The bit 255 of t is in cc. If that bit is 0, when a[] must 365 * be unchanged; otherwise, it must be replaced with t[]. 366 */ 367 cc = -cc; 368 a[0] ^= cc & (a[0] ^ t[0]); 369 a[1] ^= cc & (a[1] ^ t[1]); 370 a[2] ^= cc & (a[2] ^ t[2]); 371 a[3] ^= cc & (a[3] ^ t[3]); 372 a[4] ^= cc & (a[4] ^ t[4]); 373 } 374 375 static uint32_t 376 api_mul(unsigned char *G, size_t Glen, 377 const unsigned char *kb, size_t kblen, int curve) 378 { 379 unsigned char k[32]; 380 uint64_t x1[5], x2[5], z2[5], x3[5], z3[5]; 381 uint32_t swap; 382 int i; 383 384 (void)curve; 385 386 /* 387 * Points are encoded over exactly 32 bytes. Multipliers must fit 388 * in 32 bytes as well. 389 */ 390 if (Glen != 32 || kblen > 32) { 391 return 0; 392 } 393 394 /* 395 * RFC 7748 mandates that the high bit of the last point byte must 396 * be ignored/cleared; the "& MASK51" in the initialization for 397 * x1[4] clears that bit. 398 */ 399 x1[0] = br_dec64le(&G[0]) & MASK51; 400 x1[1] = (br_dec64le(&G[6]) >> 3) & MASK51; 401 x1[2] = (br_dec64le(&G[12]) >> 6) & MASK51; 402 x1[3] = (br_dec64le(&G[19]) >> 1) & MASK51; 403 x1[4] = (br_dec64le(&G[24]) >> 12) & MASK51; 404 405 /* 406 * We can use memset() to clear values, because exact-width types 407 * like uint64_t are guaranteed to have no padding bits or 408 * trap representations. 409 */ 410 memset(x2, 0, sizeof x2); 411 x2[0] = 1; 412 memset(z2, 0, sizeof z2); 413 memcpy(x3, x1, sizeof x1); 414 memcpy(z3, x2, sizeof x2); 415 416 /* 417 * The multiplier is provided in big-endian notation, and 418 * possibly shorter than 32 bytes. 419 */ 420 memset(k, 0, (sizeof k) - kblen); 421 memcpy(k + (sizeof k) - kblen, kb, kblen); 422 k[31] &= 0xF8; 423 k[0] &= 0x7F; 424 k[0] |= 0x40; 425 426 swap = 0; 427 428 for (i = 254; i >= 0; i --) { 429 uint64_t a[5], aa[5], b[5], bb[5], e[5]; 430 uint64_t c[5], d[5], da[5], cb[5]; 431 uint32_t kt; 432 433 kt = (k[31 - (i >> 3)] >> (i & 7)) & 1; 434 swap ^= kt; 435 f255_cswap(x2, x3, swap); 436 f255_cswap(z2, z3, swap); 437 swap = kt; 438 439 /* 440 * At that point, limbs of x_2 and z_2 are assumed to fit 441 * on at most 52 bits each. 442 * 443 * Each f255_add() adds one bit to the maximum range of 444 * the values, but f255_sub() and f255_mul() bring back 445 * the limbs into 52 bits. All f255_add() outputs are 446 * used only as inputs for f255_mul(), which ensures 447 * that limbs remain in the proper range. 448 */ 449 450 /* A = x_2 + z_2 -- limbs fit on 53 bits each */ 451 f255_add(a, x2, z2); 452 453 /* AA = A^2 */ 454 f255_mul(aa, a, a); 455 456 /* B = x_2 - z_2 */ 457 f255_sub(b, x2, z2); 458 459 /* BB = B^2 */ 460 f255_mul(bb, b, b); 461 462 /* E = AA - BB */ 463 f255_sub(e, aa, bb); 464 465 /* C = x_3 + z_3 -- limbs fit on 53 bits each */ 466 f255_add(c, x3, z3); 467 468 /* D = x_3 - z_3 */ 469 f255_sub(d, x3, z3); 470 471 /* DA = D * A */ 472 f255_mul(da, d, a); 473 474 /* CB = C * B */ 475 f255_mul(cb, c, b); 476 477 /* x_3 = (DA + CB)^2 */ 478 f255_add(x3, da, cb); 479 f255_mul(x3, x3, x3); 480 481 /* z_3 = x_1 * (DA - CB)^2 */ 482 f255_sub(z3, da, cb); 483 f255_mul(z3, z3, z3); 484 f255_mul(z3, x1, z3); 485 486 /* x_2 = AA * BB */ 487 f255_mul(x2, aa, bb); 488 489 /* z_2 = E * (AA + a24 * E) */ 490 f255_mul_a24(z2, e); 491 f255_add(z2, aa, z2); 492 f255_mul(z2, e, z2); 493 } 494 495 f255_cswap(x2, x3, swap); 496 f255_cswap(z2, z3, swap); 497 498 /* 499 * Compute 1/z2 = z2^(p-2). Since p = 2^255-19, we can mutualize 500 * most non-squarings. We use x1 and x3, now useless, as temporaries. 501 */ 502 memcpy(x1, z2, sizeof z2); 503 for (i = 0; i < 15; i ++) { 504 f255_mul(x1, x1, x1); 505 f255_mul(x1, x1, z2); 506 } 507 memcpy(x3, x1, sizeof x1); 508 for (i = 0; i < 14; i ++) { 509 int j; 510 511 for (j = 0; j < 16; j ++) { 512 f255_mul(x3, x3, x3); 513 } 514 f255_mul(x3, x3, x1); 515 } 516 for (i = 14; i >= 0; i --) { 517 f255_mul(x3, x3, x3); 518 if ((0xFFEB >> i) & 1) { 519 f255_mul(x3, z2, x3); 520 } 521 } 522 523 /* 524 * Compute x2/z2. We have 1/z2 in x3. 525 */ 526 f255_mul(x2, x2, x3); 527 f255_final_reduce(x2); 528 529 /* 530 * Encode the final x2 value in little-endian. We first assemble 531 * the limbs into 64-bit values. 532 */ 533 x2[0] |= x2[1] << 51; 534 x2[1] = (x2[1] >> 13) | (x2[2] << 38); 535 x2[2] = (x2[2] >> 26) | (x2[3] << 25); 536 x2[3] = (x2[3] >> 39) | (x2[4] << 12); 537 br_enc64le(G, x2[0]); 538 br_enc64le(G + 8, x2[1]); 539 br_enc64le(G + 16, x2[2]); 540 br_enc64le(G + 24, x2[3]); 541 return 1; 542 } 543 544 static size_t 545 api_mulgen(unsigned char *R, 546 const unsigned char *x, size_t xlen, int curve) 547 { 548 const unsigned char *G; 549 size_t Glen; 550 551 G = api_generator(curve, &Glen); 552 memcpy(R, G, Glen); 553 api_mul(R, Glen, x, xlen, curve); 554 return Glen; 555 } 556 557 static uint32_t 558 api_muladd(unsigned char *A, const unsigned char *B, size_t len, 559 const unsigned char *x, size_t xlen, 560 const unsigned char *y, size_t ylen, int curve) 561 { 562 /* 563 * We don't implement this method, since it is used for ECDSA 564 * only, and there is no ECDSA over Curve25519 (which instead 565 * uses EdDSA). 566 */ 567 (void)A; 568 (void)B; 569 (void)len; 570 (void)x; 571 (void)xlen; 572 (void)y; 573 (void)ylen; 574 (void)curve; 575 return 0; 576 } 577 578 /* see bearssl_ec.h */ 579 const br_ec_impl br_ec_c25519_m62 = { 580 (uint32_t)0x20000000, 581 &api_generator, 582 &api_order, 583 &api_xoff, 584 &api_mul, 585 &api_mulgen, 586 &api_muladd 587 }; 588 589 /* see bearssl_ec.h */ 590 const br_ec_impl * 591 br_ec_c25519_m62_get(void) 592 { 593 return &br_ec_c25519_m62; 594 } 595 596 #else 597 598 /* see bearssl_ec.h */ 599 const br_ec_impl * 600 br_ec_c25519_m62_get(void) 601 { 602 return 0; 603 } 604 605 #endif