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 33static const unsigned char P256_G[] = { 34 0x04, 0x6B, 0x17, 0xD1, 0xF2, 0xE1, 0x2C, 0x42, 0x47, 0xF8, 35 0xBC, 0xE6, 0xE5, 0x63, 0xA4, 0x40, 0xF2, 0x77, 0x03, 0x7D, 36 0x81, 0x2D, 0xEB, 0x33, 0xA0, 0xF4, 0xA1, 0x39, 0x45, 0xD8, 37 0x98, 0xC2, 0x96, 0x4F, 0xE3, 0x42, 0xE2, 0xFE, 0x1A, 0x7F, 38 0x9B, 0x8E, 0xE7, 0xEB, 0x4A, 0x7C, 0x0F, 0x9E, 0x16, 0x2B, 39 0xCE, 0x33, 0x57, 0x6B, 0x31, 0x5E, 0xCE, 0xCB, 0xB6, 0x40, 40 0x68, 0x37, 0xBF, 0x51, 0xF5 41}; 42 43static const unsigned char P256_N[] = { 44 0xFF, 0xFF, 0xFF, 0xFF, 0x00, 0x00, 0x00, 0x00, 0xFF, 0xFF, 45 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xBC, 0xE6, 0xFA, 0xAD, 46 0xA7, 0x17, 0x9E, 0x84, 0xF3, 0xB9, 0xCA, 0xC2, 0xFC, 0x63, 47 0x25, 0x51 48}; 49 50static const unsigned char * 51api_generator(int curve, size_t *len) 52{ 53 (void)curve; 54 *len = sizeof P256_G; 55 return P256_G; 56} 57 58static const unsigned char * 59api_order(int curve, size_t *len) 60{ 61 (void)curve; 62 *len = sizeof P256_N; 63 return P256_N; 64} 65 66static size_t 67api_xoff(int curve, size_t *len) 68{ 69 (void)curve; 70 *len = 32; 71 return 1; 72} 73 74/* 75 * A field element is encoded as five 64-bit integers, in basis 2^52. 76 * Limbs may occasionally exceed 2^52. 77 * 78 * A _partially reduced_ value is such that the following hold: 79 * - top limb is less than 2^48 + 2^30 80 * - the other limbs fit on 53 bits each 81 * In particular, such a value is less than twice the modulus p. 82 */ 83 84#define BIT(n) ((uint64_t)1 << (n)) 85#define MASK48 (BIT(48) - BIT(0)) 86#define MASK52 (BIT(52) - BIT(0)) 87 88/* R = 2^260 mod p */ 89static const uint64_t F256_R[] = { 90 0x0000000000010, 0xF000000000000, 0xFFFFFFFFFFFFF, 91 0xFFEFFFFFFFFFF, 0x00000000FFFFF 92}; 93 94/* Curve equation is y^2 = x^3 - 3*x + B. This constant is B*R mod p 95 (Montgomery representation of B). */ 96static const uint64_t P256_B_MONTY[] = { 97 0xDF6229C4BDDFD, 0xCA8843090D89C, 0x212ED6ACF005C, 98 0x83415A220ABF7, 0x0C30061DD4874 99}; 100 101/* 102 * Addition in the field. Carry propagation is not performed. 103 * On input, limbs may be up to 63 bits each; on output, they will 104 * be up to one bit more than on input. 105 */ 106static inline void 107f256_add(uint64_t *d, const uint64_t *a, const uint64_t *b) 108{ 109 d[0] = a[0] + b[0]; 110 d[1] = a[1] + b[1]; 111 d[2] = a[2] + b[2]; 112 d[3] = a[3] + b[3]; 113 d[4] = a[4] + b[4]; 114} 115 116/* 117 * Partially reduce the provided value. 118 * Input: limbs can go up to 61 bits each. 119 * Output: partially reduced. 120 */ 121static inline void 122f256_partial_reduce(uint64_t *a) 123{ 124 uint64_t w, cc, s; 125 126 /* 127 * Propagate carries. 128 */ 129 w = a[0]; 130 a[0] = w & MASK52; 131 cc = w >> 52; 132 w = a[1] + cc; 133 a[1] = w & MASK52; 134 cc = w >> 52; 135 w = a[2] + cc; 136 a[2] = w & MASK52; 137 cc = w >> 52; 138 w = a[3] + cc; 139 a[3] = w & MASK52; 140 cc = w >> 52; 141 a[4] += cc; 142 143 s = a[4] >> 48; /* s < 2^14 */ 144 a[0] += s; /* a[0] < 2^52 + 2^14 */ 145 w = a[1] - (s << 44); 146 a[1] = w & MASK52; /* a[1] < 2^52 */ 147 cc = -(w >> 52) & 0xFFF; /* cc < 16 */ 148 w = a[2] - cc; 149 a[2] = w & MASK52; /* a[2] < 2^52 */ 150 cc = w >> 63; /* cc = 0 or 1 */ 151 w = a[3] - cc - (s << 36); 152 a[3] = w & MASK52; /* a[3] < 2^52 */ 153 cc = w >> 63; /* cc = 0 or 1 */ 154 w = a[4] & MASK48; 155 a[4] = w + (s << 16) - cc; /* a[4] < 2^48 + 2^30 */ 156} 157 158/* 159 * Subtraction in the field. 160 * Input: limbs must fit on 60 bits each; in particular, the complete 161 * integer will be less than 2^268 + 2^217. 162 * Output: partially reduced. 163 */ 164static inline void 165f256_sub(uint64_t *d, const uint64_t *a, const uint64_t *b) 166{ 167 uint64_t t[5], w, s, cc; 168 169 /* 170 * We compute d = 2^13*p + a - b; this ensures a positive 171 * intermediate value. 172 * 173 * Each individual addition/subtraction may yield a positive or 174 * negative result; thus, we need to handle a signed carry, thus 175 * with sign extension. We prefer not to use signed types (int64_t) 176 * because conversion from unsigned to signed is cumbersome (a 177 * direct cast with the top bit set is undefined behavior; instead, 178 * we have to use pointer aliasing, using the guaranteed properties 179 * of exact-width types, but this requires the compiler to optimize 180 * away the writes and reads from RAM), and right-shifting a 181 * signed negative value is implementation-defined. Therefore, 182 * we use a custom sign extension. 183 */ 184 185 w = a[0] - b[0] - BIT(13); 186 t[0] = w & MASK52; 187 cc = w >> 52; 188 cc |= -(cc & BIT(11)); 189 w = a[1] - b[1] + cc; 190 t[1] = w & MASK52; 191 cc = w >> 52; 192 cc |= -(cc & BIT(11)); 193 w = a[2] - b[2] + cc; 194 t[2] = (w & MASK52) + BIT(5); 195 cc = w >> 52; 196 cc |= -(cc & BIT(11)); 197 w = a[3] - b[3] + cc; 198 t[3] = (w & MASK52) + BIT(49); 199 cc = w >> 52; 200 cc |= -(cc & BIT(11)); 201 t[4] = (BIT(61) - BIT(29)) + a[4] - b[4] + cc; 202 203 /* 204 * Perform partial reduction. Rule is: 205 * 2^256 = 2^224 - 2^192 - 2^96 + 1 mod p 206 * 207 * At that point: 208 * 0 <= t[0] <= 2^52 - 1 209 * 0 <= t[1] <= 2^52 - 1 210 * 2^5 <= t[2] <= 2^52 + 2^5 - 1 211 * 2^49 <= t[3] <= 2^52 + 2^49 - 1 212 * 2^59 < t[4] <= 2^61 + 2^60 - 2^29 213 * 214 * Thus, the value 's' (t[4] / 2^48) will be necessarily 215 * greater than 2048, and less than 12288. 216 */ 217 s = t[4] >> 48; 218 219 d[0] = t[0] + s; /* d[0] <= 2^52 + 12287 */ 220 w = t[1] - (s << 44); 221 d[1] = w & MASK52; /* d[1] <= 2^52 - 1 */ 222 cc = -(w >> 52) & 0xFFF; /* cc <= 48 */ 223 w = t[2] - cc; 224 cc = w >> 63; /* cc = 0 or 1 */ 225 d[2] = w + (cc << 52); /* d[2] <= 2^52 + 31 */ 226 w = t[3] - cc - (s << 36); 227 cc = w >> 63; /* cc = 0 or 1 */ 228 d[3] = w + (cc << 52); /* t[3] <= 2^52 + 2^49 - 1 */ 229 d[4] = (t[4] & MASK48) + (s << 16) - cc; /* d[4] < 2^48 + 2^30 */ 230 231 /* 232 * If s = 0, then none of the limbs is modified, and there cannot 233 * be an overflow; if s != 0, then (s << 16) > cc, and there is 234 * no overflow either. 235 */ 236} 237 238/* 239 * Montgomery multiplication in the field. 240 * Input: limbs must fit on 56 bits each. 241 * Output: partially reduced. 242 */ 243static void 244f256_montymul(uint64_t *d, const uint64_t *a, const uint64_t *b) 245{ 246#if BR_INT128 247 248 int i; 249 uint64_t t[5]; 250 251 t[0] = 0; 252 t[1] = 0; 253 t[2] = 0; 254 t[3] = 0; 255 t[4] = 0; 256 for (i = 0; i < 5; i ++) { 257 uint64_t x, f, cc, w, s; 258 unsigned __int128 z; 259 260 /* 261 * Since limbs of a[] and b[] fit on 56 bits each, 262 * each individual product fits on 112 bits. Also, 263 * the factor f fits on 52 bits, so f<<48 fits on 264 * 112 bits too. This guarantees that carries (cc) 265 * will fit on 62 bits, thus no overflow. 266 * 267 * The operations below compute: 268 * t <- (t + x*b + f*p) / 2^64 269 */ 270 x = a[i]; 271 z = (unsigned __int128)b[0] * (unsigned __int128)x 272 + (unsigned __int128)t[0]; 273 f = (uint64_t)z & MASK52; 274 cc = (uint64_t)(z >> 52); 275 z = (unsigned __int128)b[1] * (unsigned __int128)x 276 + (unsigned __int128)t[1] + cc 277 + ((unsigned __int128)f << 44); 278 t[0] = (uint64_t)z & MASK52; 279 cc = (uint64_t)(z >> 52); 280 z = (unsigned __int128)b[2] * (unsigned __int128)x 281 + (unsigned __int128)t[2] + cc; 282 t[1] = (uint64_t)z & MASK52; 283 cc = (uint64_t)(z >> 52); 284 z = (unsigned __int128)b[3] * (unsigned __int128)x 285 + (unsigned __int128)t[3] + cc 286 + ((unsigned __int128)f << 36); 287 t[2] = (uint64_t)z & MASK52; 288 cc = (uint64_t)(z >> 52); 289 z = (unsigned __int128)b[4] * (unsigned __int128)x 290 + (unsigned __int128)t[4] + cc 291 + ((unsigned __int128)f << 48) 292 - ((unsigned __int128)f << 16); 293 t[3] = (uint64_t)z & MASK52; 294 t[4] = (uint64_t)(z >> 52); 295 296 /* 297 * t[4] may be up to 62 bits here; we need to do a 298 * partial reduction. Note that limbs t[0] to t[3] 299 * fit on 52 bits each. 300 */ 301 s = t[4] >> 48; /* s < 2^14 */ 302 t[0] += s; /* t[0] < 2^52 + 2^14 */ 303 w = t[1] - (s << 44); 304 t[1] = w & MASK52; /* t[1] < 2^52 */ 305 cc = -(w >> 52) & 0xFFF; /* cc < 16 */ 306 w = t[2] - cc; 307 t[2] = w & MASK52; /* t[2] < 2^52 */ 308 cc = w >> 63; /* cc = 0 or 1 */ 309 w = t[3] - cc - (s << 36); 310 t[3] = w & MASK52; /* t[3] < 2^52 */ 311 cc = w >> 63; /* cc = 0 or 1 */ 312 w = t[4] & MASK48; 313 t[4] = w + (s << 16) - cc; /* t[4] < 2^48 + 2^30 */ 314 315 /* 316 * The final t[4] cannot overflow because cc is 0 or 1, 317 * and cc can be 1 only if s != 0. 318 */ 319 } 320 321 d[0] = t[0]; 322 d[1] = t[1]; 323 d[2] = t[2]; 324 d[3] = t[3]; 325 d[4] = t[4]; 326 327#elif BR_UMUL128 328 329 int i; 330 uint64_t t[5]; 331 332 t[0] = 0; 333 t[1] = 0; 334 t[2] = 0; 335 t[3] = 0; 336 t[4] = 0; 337 for (i = 0; i < 5; i ++) { 338 uint64_t x, f, cc, w, s, zh, zl; 339 unsigned char k; 340 341 /* 342 * Since limbs of a[] and b[] fit on 56 bits each, 343 * each individual product fits on 112 bits. Also, 344 * the factor f fits on 52 bits, so f<<48 fits on 345 * 112 bits too. This guarantees that carries (cc) 346 * will fit on 62 bits, thus no overflow. 347 * 348 * The operations below compute: 349 * t <- (t + x*b + f*p) / 2^64 350 */ 351 x = a[i]; 352 zl = _umul128(b[0], x, &zh); 353 k = _addcarry_u64(0, t[0], zl, &zl); 354 (void)_addcarry_u64(k, 0, zh, &zh); 355 f = zl & MASK52; 356 cc = (zl >> 52) | (zh << 12); 357 358 zl = _umul128(b[1], x, &zh); 359 k = _addcarry_u64(0, t[1], zl, &zl); 360 (void)_addcarry_u64(k, 0, zh, &zh); 361 k = _addcarry_u64(0, cc, zl, &zl); 362 (void)_addcarry_u64(k, 0, zh, &zh); 363 k = _addcarry_u64(0, f << 44, zl, &zl); 364 (void)_addcarry_u64(k, f >> 20, zh, &zh); 365 t[0] = zl & MASK52; 366 cc = (zl >> 52) | (zh << 12); 367 368 zl = _umul128(b[2], x, &zh); 369 k = _addcarry_u64(0, t[2], zl, &zl); 370 (void)_addcarry_u64(k, 0, zh, &zh); 371 k = _addcarry_u64(0, cc, zl, &zl); 372 (void)_addcarry_u64(k, 0, zh, &zh); 373 t[1] = zl & MASK52; 374 cc = (zl >> 52) | (zh << 12); 375 376 zl = _umul128(b[3], x, &zh); 377 k = _addcarry_u64(0, t[3], zl, &zl); 378 (void)_addcarry_u64(k, 0, zh, &zh); 379 k = _addcarry_u64(0, cc, zl, &zl); 380 (void)_addcarry_u64(k, 0, zh, &zh); 381 k = _addcarry_u64(0, f << 36, zl, &zl); 382 (void)_addcarry_u64(k, f >> 28, zh, &zh); 383 t[2] = zl & MASK52; 384 cc = (zl >> 52) | (zh << 12); 385 386 zl = _umul128(b[4], x, &zh); 387 k = _addcarry_u64(0, t[4], zl, &zl); 388 (void)_addcarry_u64(k, 0, zh, &zh); 389 k = _addcarry_u64(0, cc, zl, &zl); 390 (void)_addcarry_u64(k, 0, zh, &zh); 391 k = _addcarry_u64(0, f << 48, zl, &zl); 392 (void)_addcarry_u64(k, f >> 16, zh, &zh); 393 k = _subborrow_u64(0, zl, f << 16, &zl); 394 (void)_subborrow_u64(k, zh, f >> 48, &zh); 395 t[3] = zl & MASK52; 396 t[4] = (zl >> 52) | (zh << 12); 397 398 /* 399 * t[4] may be up to 62 bits here; we need to do a 400 * partial reduction. Note that limbs t[0] to t[3] 401 * fit on 52 bits each. 402 */ 403 s = t[4] >> 48; /* s < 2^14 */ 404 t[0] += s; /* t[0] < 2^52 + 2^14 */ 405 w = t[1] - (s << 44); 406 t[1] = w & MASK52; /* t[1] < 2^52 */ 407 cc = -(w >> 52) & 0xFFF; /* cc < 16 */ 408 w = t[2] - cc; 409 t[2] = w & MASK52; /* t[2] < 2^52 */ 410 cc = w >> 63; /* cc = 0 or 1 */ 411 w = t[3] - cc - (s << 36); 412 t[3] = w & MASK52; /* t[3] < 2^52 */ 413 cc = w >> 63; /* cc = 0 or 1 */ 414 w = t[4] & MASK48; 415 t[4] = w + (s << 16) - cc; /* t[4] < 2^48 + 2^30 */ 416 417 /* 418 * The final t[4] cannot overflow because cc is 0 or 1, 419 * and cc can be 1 only if s != 0. 420 */ 421 } 422 423 d[0] = t[0]; 424 d[1] = t[1]; 425 d[2] = t[2]; 426 d[3] = t[3]; 427 d[4] = t[4]; 428 429#endif 430} 431 432/* 433 * Montgomery squaring in the field; currently a basic wrapper around 434 * multiplication (inline, should be optimized away). 435 * TODO: see if some extra speed can be gained here. 436 */ 437static inline void 438f256_montysquare(uint64_t *d, const uint64_t *a) 439{ 440 f256_montymul(d, a, a); 441} 442 443/* 444 * Convert to Montgomery representation. 445 */ 446static void 447f256_tomonty(uint64_t *d, const uint64_t *a) 448{ 449 /* 450 * R2 = 2^520 mod p. 451 * If R = 2^260 mod p, then R2 = R^2 mod p; and the Montgomery 452 * multiplication of a by R2 is: a*R2/R = a*R mod p, i.e. the 453 * conversion to Montgomery representation. 454 */ 455 static const uint64_t R2[] = { 456 0x0000000000300, 0xFFFFFFFF00000, 0xFFFFEFFFFFFFB, 457 0xFDFFFFFFFFFFF, 0x0000004FFFFFF 458 }; 459 460 f256_montymul(d, a, R2); 461} 462 463/* 464 * Convert from Montgomery representation. 465 */ 466static void 467f256_frommonty(uint64_t *d, const uint64_t *a) 468{ 469 /* 470 * Montgomery multiplication by 1 is division by 2^260 modulo p. 471 */ 472 static const uint64_t one[] = { 1, 0, 0, 0, 0 }; 473 474 f256_montymul(d, a, one); 475} 476 477/* 478 * Inversion in the field. If the source value is 0 modulo p, then this 479 * returns 0 or p. This function uses Montgomery representation. 480 */ 481static void 482f256_invert(uint64_t *d, const uint64_t *a) 483{ 484 /* 485 * We compute a^(p-2) mod p. The exponent pattern (from high to 486 * low) is: 487 * - 32 bits of value 1 488 * - 31 bits of value 0 489 * - 1 bit of value 1 490 * - 96 bits of value 0 491 * - 94 bits of value 1 492 * - 1 bit of value 0 493 * - 1 bit of value 1 494 * To speed up the square-and-multiply algorithm, we precompute 495 * a^(2^31-1). 496 */ 497 498 uint64_t r[5], t[5]; 499 int i; 500 501 memcpy(t, a, sizeof t); 502 for (i = 0; i < 30; i ++) { 503 f256_montysquare(t, t); 504 f256_montymul(t, t, a); 505 } 506 507 memcpy(r, t, sizeof t); 508 for (i = 224; i >= 0; i --) { 509 f256_montysquare(r, r); 510 switch (i) { 511 case 0: 512 case 2: 513 case 192: 514 case 224: 515 f256_montymul(r, r, a); 516 break; 517 case 3: 518 case 34: 519 case 65: 520 f256_montymul(r, r, t); 521 break; 522 } 523 } 524 memcpy(d, r, sizeof r); 525} 526 527/* 528 * Finalize reduction. 529 * Input value should be partially reduced. 530 * On output, limbs a[0] to a[3] fit on 52 bits each, limb a[4] fits 531 * on 48 bits, and the integer is less than p. 532 */ 533static inline void 534f256_final_reduce(uint64_t *a) 535{ 536 uint64_t r[5], t[5], w, cc; 537 int i; 538 539 /* 540 * Propagate carries to ensure that limbs 0 to 3 fit on 52 bits. 541 */ 542 cc = 0; 543 for (i = 0; i < 5; i ++) { 544 w = a[i] + cc; 545 r[i] = w & MASK52; 546 cc = w >> 52; 547 } 548 549 /* 550 * We compute t = r + (2^256 - p) = r + 2^224 - 2^192 - 2^96 + 1. 551 * If t < 2^256, then r < p, and we return r. Otherwise, we 552 * want to return r - p = t - 2^256. 553 */ 554 555 /* 556 * Add 2^224 + 1, and propagate carries to ensure that limbs 557 * t[0] to t[3] fit in 52 bits each. 558 */ 559 w = r[0] + 1; 560 t[0] = w & MASK52; 561 cc = w >> 52; 562 w = r[1] + cc; 563 t[1] = w & MASK52; 564 cc = w >> 52; 565 w = r[2] + cc; 566 t[2] = w & MASK52; 567 cc = w >> 52; 568 w = r[3] + cc; 569 t[3] = w & MASK52; 570 cc = w >> 52; 571 t[4] = r[4] + cc + BIT(16); 572 573 /* 574 * Subtract 2^192 + 2^96. Since we just added 2^224 + 1, the 575 * result cannot be negative. 576 */ 577 w = t[1] - BIT(44); 578 t[1] = w & MASK52; 579 cc = w >> 63; 580 w = t[2] - cc; 581 t[2] = w & MASK52; 582 cc = w >> 63; 583 w = t[3] - BIT(36); 584 t[3] = w & MASK52; 585 cc = w >> 63; 586 t[4] -= cc; 587 588 /* 589 * If the top limb t[4] fits on 48 bits, then r[] is already 590 * in the proper range. Otherwise, t[] is the value to return 591 * (truncated to 256 bits). 592 */ 593 cc = -(t[4] >> 48); 594 t[4] &= MASK48; 595 for (i = 0; i < 5; i ++) { 596 a[i] = r[i] ^ (cc & (r[i] ^ t[i])); 597 } 598} 599 600/* 601 * Points in affine and Jacobian coordinates. 602 * 603 * - In affine coordinates, the point-at-infinity cannot be encoded. 604 * - Jacobian coordinates (X,Y,Z) correspond to affine (X/Z^2,Y/Z^3); 605 * if Z = 0 then this is the point-at-infinity. 606 */ 607typedef struct { 608 uint64_t x[5]; 609 uint64_t y[5]; 610} p256_affine; 611 612typedef struct { 613 uint64_t x[5]; 614 uint64_t y[5]; 615 uint64_t z[5]; 616} p256_jacobian; 617 618/* 619 * Decode a field element (unsigned big endian notation). 620 */ 621static void 622f256_decode(uint64_t *a, const unsigned char *buf) 623{ 624 uint64_t w0, w1, w2, w3; 625 626 w3 = br_dec64be(buf + 0); 627 w2 = br_dec64be(buf + 8); 628 w1 = br_dec64be(buf + 16); 629 w0 = br_dec64be(buf + 24); 630 a[0] = w0 & MASK52; 631 a[1] = ((w0 >> 52) | (w1 << 12)) & MASK52; 632 a[2] = ((w1 >> 40) | (w2 << 24)) & MASK52; 633 a[3] = ((w2 >> 28) | (w3 << 36)) & MASK52; 634 a[4] = w3 >> 16; 635} 636 637/* 638 * Encode a field element (unsigned big endian notation). The field 639 * element MUST be fully reduced. 640 */ 641static void 642f256_encode(unsigned char *buf, const uint64_t *a) 643{ 644 uint64_t w0, w1, w2, w3; 645 646 w0 = a[0] | (a[1] << 52); 647 w1 = (a[1] >> 12) | (a[2] << 40); 648 w2 = (a[2] >> 24) | (a[3] << 28); 649 w3 = (a[3] >> 36) | (a[4] << 16); 650 br_enc64be(buf + 0, w3); 651 br_enc64be(buf + 8, w2); 652 br_enc64be(buf + 16, w1); 653 br_enc64be(buf + 24, w0); 654} 655 656/* 657 * Decode a point. The returned point is in Jacobian coordinates, but 658 * with z = 1. If the encoding is invalid, or encodes a point which is 659 * not on the curve, or encodes the point at infinity, then this function 660 * returns 0. Otherwise, 1 is returned. 661 * 662 * The buffer is assumed to have length exactly 65 bytes. 663 */ 664static uint32_t 665point_decode(p256_jacobian *P, const unsigned char *buf) 666{ 667 uint64_t x[5], y[5], t[5], x3[5], tt; 668 uint32_t r; 669 670 /* 671 * Header byte shall be 0x04. 672 */ 673 r = EQ(buf[0], 0x04); 674 675 /* 676 * Decode X and Y coordinates, and convert them into 677 * Montgomery representation. 678 */ 679 f256_decode(x, buf + 1); 680 f256_decode(y, buf + 33); 681 f256_tomonty(x, x); 682 f256_tomonty(y, y); 683 684 /* 685 * Verify y^2 = x^3 + A*x + B. In curve P-256, A = -3. 686 * Note that the Montgomery representation of 0 is 0. We must 687 * take care to apply the final reduction to make sure we have 688 * 0 and not p. 689 */ 690 f256_montysquare(t, y); 691 f256_montysquare(x3, x); 692 f256_montymul(x3, x3, x); 693 f256_sub(t, t, x3); 694 f256_add(t, t, x); 695 f256_add(t, t, x); 696 f256_add(t, t, x); 697 f256_sub(t, t, P256_B_MONTY); 698 f256_final_reduce(t); 699 tt = t[0] | t[1] | t[2] | t[3] | t[4]; 700 r &= EQ((uint32_t)(tt | (tt >> 32)), 0); 701 702 /* 703 * Return the point in Jacobian coordinates (and Montgomery 704 * representation). 705 */ 706 memcpy(P->x, x, sizeof x); 707 memcpy(P->y, y, sizeof y); 708 memcpy(P->z, F256_R, sizeof F256_R); 709 return r; 710} 711 712/* 713 * Final conversion for a point: 714 * - The point is converted back to affine coordinates. 715 * - Final reduction is performed. 716 * - The point is encoded into the provided buffer. 717 * 718 * If the point is the point-at-infinity, all operations are performed, 719 * but the buffer contents are indeterminate, and 0 is returned. Otherwise, 720 * the encoded point is written in the buffer, and 1 is returned. 721 */ 722static uint32_t 723point_encode(unsigned char *buf, const p256_jacobian *P) 724{ 725 uint64_t t1[5], t2[5], z; 726 727 /* Set t1 = 1/z^2 and t2 = 1/z^3. */ 728 f256_invert(t2, P->z); 729 f256_montysquare(t1, t2); 730 f256_montymul(t2, t2, t1); 731 732 /* Compute affine coordinates x (in t1) and y (in t2). */ 733 f256_montymul(t1, P->x, t1); 734 f256_montymul(t2, P->y, t2); 735 736 /* Convert back from Montgomery representation, and finalize 737 reductions. */ 738 f256_frommonty(t1, t1); 739 f256_frommonty(t2, t2); 740 f256_final_reduce(t1); 741 f256_final_reduce(t2); 742 743 /* Encode. */ 744 buf[0] = 0x04; 745 f256_encode(buf + 1, t1); 746 f256_encode(buf + 33, t2); 747 748 /* Return success if and only if P->z != 0. */ 749 z = P->z[0] | P->z[1] | P->z[2] | P->z[3] | P->z[4]; 750 return NEQ((uint32_t)(z | z >> 32), 0); 751} 752 753/* 754 * Point doubling in Jacobian coordinates: point P is doubled. 755 * Note: if the source point is the point-at-infinity, then the result is 756 * still the point-at-infinity, which is correct. Moreover, if the three 757 * coordinates were zero, then they still are zero in the returned value. 758 */ 759static void 760p256_double(p256_jacobian *P) 761{ 762 /* 763 * Doubling formulas are: 764 * 765 * s = 4*x*y^2 766 * m = 3*(x + z^2)*(x - z^2) 767 * x' = m^2 - 2*s 768 * y' = m*(s - x') - 8*y^4 769 * z' = 2*y*z 770 * 771 * These formulas work for all points, including points of order 2 772 * and points at infinity: 773 * - If y = 0 then z' = 0. But there is no such point in P-256 774 * anyway. 775 * - If z = 0 then z' = 0. 776 */ 777 uint64_t t1[5], t2[5], t3[5], t4[5]; 778 779 /* 780 * Compute z^2 in t1. 781 */ 782 f256_montysquare(t1, P->z); 783 784 /* 785 * Compute x-z^2 in t2 and x+z^2 in t1. 786 */ 787 f256_add(t2, P->x, t1); 788 f256_sub(t1, P->x, t1); 789 790 /* 791 * Compute 3*(x+z^2)*(x-z^2) in t1. 792 */ 793 f256_montymul(t3, t1, t2); 794 f256_add(t1, t3, t3); 795 f256_add(t1, t3, t1); 796 797 /* 798 * Compute 4*x*y^2 (in t2) and 2*y^2 (in t3). 799 */ 800 f256_montysquare(t3, P->y); 801 f256_add(t3, t3, t3); 802 f256_montymul(t2, P->x, t3); 803 f256_add(t2, t2, t2); 804 805 /* 806 * Compute x' = m^2 - 2*s. 807 */ 808 f256_montysquare(P->x, t1); 809 f256_sub(P->x, P->x, t2); 810 f256_sub(P->x, P->x, t2); 811 812 /* 813 * Compute z' = 2*y*z. 814 */ 815 f256_montymul(t4, P->y, P->z); 816 f256_add(P->z, t4, t4); 817 f256_partial_reduce(P->z); 818 819 /* 820 * Compute y' = m*(s - x') - 8*y^4. Note that we already have 821 * 2*y^2 in t3. 822 */ 823 f256_sub(t2, t2, P->x); 824 f256_montymul(P->y, t1, t2); 825 f256_montysquare(t4, t3); 826 f256_add(t4, t4, t4); 827 f256_sub(P->y, P->y, t4); 828} 829 830/* 831 * Point addition (Jacobian coordinates): P1 is replaced with P1+P2. 832 * This function computes the wrong result in the following cases: 833 * 834 * - If P1 == 0 but P2 != 0 835 * - If P1 != 0 but P2 == 0 836 * - If P1 == P2 837 * 838 * In all three cases, P1 is set to the point at infinity. 839 * 840 * Returned value is 0 if one of the following occurs: 841 * 842 * - P1 and P2 have the same Y coordinate. 843 * - P1 == 0 and P2 == 0. 844 * - The Y coordinate of one of the points is 0 and the other point is 845 * the point at infinity. 846 * 847 * The third case cannot actually happen with valid points, since a point 848 * with Y == 0 is a point of order 2, and there is no point of order 2 on 849 * curve P-256. 850 * 851 * Therefore, assuming that P1 != 0 and P2 != 0 on input, then the caller 852 * can apply the following: 853 * 854 * - If the result is not the point at infinity, then it is correct. 855 * - Otherwise, if the returned value is 1, then this is a case of 856 * P1+P2 == 0, so the result is indeed the point at infinity. 857 * - Otherwise, P1 == P2, so a "double" operation should have been 858 * performed. 859 * 860 * Note that you can get a returned value of 0 with a correct result, 861 * e.g. if P1 and P2 have the same Y coordinate, but distinct X coordinates. 862 */ 863static uint32_t 864p256_add(p256_jacobian *P1, const p256_jacobian *P2) 865{ 866 /* 867 * Addtions formulas are: 868 * 869 * u1 = x1 * z2^2 870 * u2 = x2 * z1^2 871 * s1 = y1 * z2^3 872 * s2 = y2 * z1^3 873 * h = u2 - u1 874 * r = s2 - s1 875 * x3 = r^2 - h^3 - 2 * u1 * h^2 876 * y3 = r * (u1 * h^2 - x3) - s1 * h^3 877 * z3 = h * z1 * z2 878 */ 879 uint64_t t1[5], t2[5], t3[5], t4[5], t5[5], t6[5], t7[5], tt; 880 uint32_t ret; 881 882 /* 883 * Compute u1 = x1*z2^2 (in t1) and s1 = y1*z2^3 (in t3). 884 */ 885 f256_montysquare(t3, P2->z); 886 f256_montymul(t1, P1->x, t3); 887 f256_montymul(t4, P2->z, t3); 888 f256_montymul(t3, P1->y, t4); 889 890 /* 891 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4). 892 */ 893 f256_montysquare(t4, P1->z); 894 f256_montymul(t2, P2->x, t4); 895 f256_montymul(t5, P1->z, t4); 896 f256_montymul(t4, P2->y, t5); 897 898 /* 899 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4). 900 * We need to test whether r is zero, so we will do some extra 901 * reduce. 902 */ 903 f256_sub(t2, t2, t1); 904 f256_sub(t4, t4, t3); 905 f256_final_reduce(t4); 906 tt = t4[0] | t4[1] | t4[2] | t4[3] | t4[4]; 907 ret = (uint32_t)(tt | (tt >> 32)); 908 ret = (ret | -ret) >> 31; 909 910 /* 911 * Compute u1*h^2 (in t6) and h^3 (in t5); 912 */ 913 f256_montysquare(t7, t2); 914 f256_montymul(t6, t1, t7); 915 f256_montymul(t5, t7, t2); 916 917 /* 918 * Compute x3 = r^2 - h^3 - 2*u1*h^2. 919 */ 920 f256_montysquare(P1->x, t4); 921 f256_sub(P1->x, P1->x, t5); 922 f256_sub(P1->x, P1->x, t6); 923 f256_sub(P1->x, P1->x, t6); 924 925 /* 926 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3. 927 */ 928 f256_sub(t6, t6, P1->x); 929 f256_montymul(P1->y, t4, t6); 930 f256_montymul(t1, t5, t3); 931 f256_sub(P1->y, P1->y, t1); 932 933 /* 934 * Compute z3 = h*z1*z2. 935 */ 936 f256_montymul(t1, P1->z, P2->z); 937 f256_montymul(P1->z, t1, t2); 938 939 return ret; 940} 941 942/* 943 * Point addition (mixed coordinates): P1 is replaced with P1+P2. 944 * This is a specialised function for the case when P2 is a non-zero point 945 * in affine coordinates. 946 * 947 * This function computes the wrong result in the following cases: 948 * 949 * - If P1 == 0 950 * - If P1 == P2 951 * 952 * In both cases, P1 is set to the point at infinity. 953 * 954 * Returned value is 0 if one of the following occurs: 955 * 956 * - P1 and P2 have the same Y (affine) coordinate. 957 * - The Y coordinate of P2 is 0 and P1 is the point at infinity. 958 * 959 * The second case cannot actually happen with valid points, since a point 960 * with Y == 0 is a point of order 2, and there is no point of order 2 on 961 * curve P-256. 962 * 963 * Therefore, assuming that P1 != 0 on input, then the caller 964 * can apply the following: 965 * 966 * - If the result is not the point at infinity, then it is correct. 967 * - Otherwise, if the returned value is 1, then this is a case of 968 * P1+P2 == 0, so the result is indeed the point at infinity. 969 * - Otherwise, P1 == P2, so a "double" operation should have been 970 * performed. 971 * 972 * Again, a value of 0 may be returned in some cases where the addition 973 * result is correct. 974 */ 975static uint32_t 976p256_add_mixed(p256_jacobian *P1, const p256_affine *P2) 977{ 978 /* 979 * Addtions formulas are: 980 * 981 * u1 = x1 982 * u2 = x2 * z1^2 983 * s1 = y1 984 * s2 = y2 * z1^3 985 * h = u2 - u1 986 * r = s2 - s1 987 * x3 = r^2 - h^3 - 2 * u1 * h^2 988 * y3 = r * (u1 * h^2 - x3) - s1 * h^3 989 * z3 = h * z1 990 */ 991 uint64_t t1[5], t2[5], t3[5], t4[5], t5[5], t6[5], t7[5], tt; 992 uint32_t ret; 993 994 /* 995 * Compute u1 = x1 (in t1) and s1 = y1 (in t3). 996 */ 997 memcpy(t1, P1->x, sizeof t1); 998 memcpy(t3, P1->y, sizeof t3); 999 1000 /* 1001 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4). 1002 */ 1003 f256_montysquare(t4, P1->z); 1004 f256_montymul(t2, P2->x, t4); 1005 f256_montymul(t5, P1->z, t4); 1006 f256_montymul(t4, P2->y, t5); 1007 1008 /* 1009 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4). 1010 * We need to test whether r is zero, so we will do some extra 1011 * reduce. 1012 */ 1013 f256_sub(t2, t2, t1); 1014 f256_sub(t4, t4, t3); 1015 f256_final_reduce(t4); 1016 tt = t4[0] | t4[1] | t4[2] | t4[3] | t4[4]; 1017 ret = (uint32_t)(tt | (tt >> 32)); 1018 ret = (ret | -ret) >> 31; 1019 1020 /* 1021 * Compute u1*h^2 (in t6) and h^3 (in t5); 1022 */ 1023 f256_montysquare(t7, t2); 1024 f256_montymul(t6, t1, t7); 1025 f256_montymul(t5, t7, t2); 1026 1027 /* 1028 * Compute x3 = r^2 - h^3 - 2*u1*h^2. 1029 */ 1030 f256_montysquare(P1->x, t4); 1031 f256_sub(P1->x, P1->x, t5); 1032 f256_sub(P1->x, P1->x, t6); 1033 f256_sub(P1->x, P1->x, t6); 1034 1035 /* 1036 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3. 1037 */ 1038 f256_sub(t6, t6, P1->x); 1039 f256_montymul(P1->y, t4, t6); 1040 f256_montymul(t1, t5, t3); 1041 f256_sub(P1->y, P1->y, t1); 1042 1043 /* 1044 * Compute z3 = h*z1*z2. 1045 */ 1046 f256_montymul(P1->z, P1->z, t2); 1047 1048 return ret; 1049} 1050 1051#if 0 1052/* unused */ 1053/* 1054 * Point addition (mixed coordinates, complete): P1 is replaced with P1+P2. 1055 * This is a specialised function for the case when P2 is a non-zero point 1056 * in affine coordinates. 1057 * 1058 * This function returns the correct result in all cases. 1059 */ 1060static uint32_t 1061p256_add_complete_mixed(p256_jacobian *P1, const p256_affine *P2) 1062{ 1063 /* 1064 * Addtions formulas, in the general case, are: 1065 * 1066 * u1 = x1 1067 * u2 = x2 * z1^2 1068 * s1 = y1 1069 * s2 = y2 * z1^3 1070 * h = u2 - u1 1071 * r = s2 - s1 1072 * x3 = r^2 - h^3 - 2 * u1 * h^2 1073 * y3 = r * (u1 * h^2 - x3) - s1 * h^3 1074 * z3 = h * z1 1075 * 1076 * These formulas mishandle the two following cases: 1077 * 1078 * - If P1 is the point-at-infinity (z1 = 0), then z3 is 1079 * incorrectly set to 0. 1080 * 1081 * - If P1 = P2, then u1 = u2 and s1 = s2, and x3, y3 and z3 1082 * are all set to 0. 1083 * 1084 * However, if P1 + P2 = 0, then u1 = u2 but s1 != s2, and then 1085 * we correctly get z3 = 0 (the point-at-infinity). 1086 * 1087 * To fix the case P1 = 0, we perform at the end a copy of P2 1088 * over P1, conditional to z1 = 0. 1089 * 1090 * For P1 = P2: in that case, both h and r are set to 0, and 1091 * we get x3, y3 and z3 equal to 0. We can test for that 1092 * occurrence to make a mask which will be all-one if P1 = P2, 1093 * or all-zero otherwise; then we can compute the double of P2 1094 * and add it, combined with the mask, to (x3,y3,z3). 1095 * 1096 * Using the doubling formulas in p256_double() on (x2,y2), 1097 * simplifying since P2 is affine (i.e. z2 = 1, implicitly), 1098 * we get: 1099 * s = 4*x2*y2^2 1100 * m = 3*(x2 + 1)*(x2 - 1) 1101 * x' = m^2 - 2*s 1102 * y' = m*(s - x') - 8*y2^4 1103 * z' = 2*y2 1104 * which requires only 6 multiplications. Added to the 11 1105 * multiplications of the normal mixed addition in Jacobian 1106 * coordinates, we get a cost of 17 multiplications in total. 1107 */ 1108 uint64_t t1[5], t2[5], t3[5], t4[5], t5[5], t6[5], t7[5], tt, zz; 1109 int i; 1110 1111 /* 1112 * Set zz to -1 if P1 is the point at infinity, 0 otherwise. 1113 */ 1114 zz = P1->z[0] | P1->z[1] | P1->z[2] | P1->z[3] | P1->z[4]; 1115 zz = ((zz | -zz) >> 63) - (uint64_t)1; 1116 1117 /* 1118 * Compute u1 = x1 (in t1) and s1 = y1 (in t3). 1119 */ 1120 memcpy(t1, P1->x, sizeof t1); 1121 memcpy(t3, P1->y, sizeof t3); 1122 1123 /* 1124 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4). 1125 */ 1126 f256_montysquare(t4, P1->z); 1127 f256_montymul(t2, P2->x, t4); 1128 f256_montymul(t5, P1->z, t4); 1129 f256_montymul(t4, P2->y, t5); 1130 1131 /* 1132 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4). 1133 * reduce. 1134 */ 1135 f256_sub(t2, t2, t1); 1136 f256_sub(t4, t4, t3); 1137 1138 /* 1139 * If both h = 0 and r = 0, then P1 = P2, and we want to set 1140 * the mask tt to -1; otherwise, the mask will be 0. 1141 */ 1142 f256_final_reduce(t2); 1143 f256_final_reduce(t4); 1144 tt = t2[0] | t2[1] | t2[2] | t2[3] | t2[4] 1145 | t4[0] | t4[1] | t4[2] | t4[3] | t4[4]; 1146 tt = ((tt | -tt) >> 63) - (uint64_t)1; 1147 1148 /* 1149 * Compute u1*h^2 (in t6) and h^3 (in t5); 1150 */ 1151 f256_montysquare(t7, t2); 1152 f256_montymul(t6, t1, t7); 1153 f256_montymul(t5, t7, t2); 1154 1155 /* 1156 * Compute x3 = r^2 - h^3 - 2*u1*h^2. 1157 */ 1158 f256_montysquare(P1->x, t4); 1159 f256_sub(P1->x, P1->x, t5); 1160 f256_sub(P1->x, P1->x, t6); 1161 f256_sub(P1->x, P1->x, t6); 1162 1163 /* 1164 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3. 1165 */ 1166 f256_sub(t6, t6, P1->x); 1167 f256_montymul(P1->y, t4, t6); 1168 f256_montymul(t1, t5, t3); 1169 f256_sub(P1->y, P1->y, t1); 1170 1171 /* 1172 * Compute z3 = h*z1. 1173 */ 1174 f256_montymul(P1->z, P1->z, t2); 1175 1176 /* 1177 * The "double" result, in case P1 = P2. 1178 */ 1179 1180 /* 1181 * Compute z' = 2*y2 (in t1). 1182 */ 1183 f256_add(t1, P2->y, P2->y); 1184 f256_partial_reduce(t1); 1185 1186 /* 1187 * Compute 2*(y2^2) (in t2) and s = 4*x2*(y2^2) (in t3). 1188 */ 1189 f256_montysquare(t2, P2->y); 1190 f256_add(t2, t2, t2); 1191 f256_add(t3, t2, t2); 1192 f256_montymul(t3, P2->x, t3); 1193 1194 /* 1195 * Compute m = 3*(x2^2 - 1) (in t4). 1196 */ 1197 f256_montysquare(t4, P2->x); 1198 f256_sub(t4, t4, F256_R); 1199 f256_add(t5, t4, t4); 1200 f256_add(t4, t4, t5); 1201 1202 /* 1203 * Compute x' = m^2 - 2*s (in t5). 1204 */ 1205 f256_montysquare(t5, t4); 1206 f256_sub(t5, t3); 1207 f256_sub(t5, t3); 1208 1209 /* 1210 * Compute y' = m*(s - x') - 8*y2^4 (in t6). 1211 */ 1212 f256_sub(t6, t3, t5); 1213 f256_montymul(t6, t6, t4); 1214 f256_montysquare(t7, t2); 1215 f256_sub(t6, t6, t7); 1216 f256_sub(t6, t6, t7); 1217 1218 /* 1219 * We now have the alternate (doubling) coordinates in (t5,t6,t1). 1220 * We combine them with (x3,y3,z3). 1221 */ 1222 for (i = 0; i < 5; i ++) { 1223 P1->x[i] |= tt & t5[i]; 1224 P1->y[i] |= tt & t6[i]; 1225 P1->z[i] |= tt & t1[i]; 1226 } 1227 1228 /* 1229 * If P1 = 0, then we get z3 = 0 (which is invalid); if z1 is 0, 1230 * then we want to replace the result with a copy of P2. The 1231 * test on z1 was done at the start, in the zz mask. 1232 */ 1233 for (i = 0; i < 5; i ++) { 1234 P1->x[i] ^= zz & (P1->x[i] ^ P2->x[i]); 1235 P1->y[i] ^= zz & (P1->y[i] ^ P2->y[i]); 1236 P1->z[i] ^= zz & (P1->z[i] ^ F256_R[i]); 1237 } 1238} 1239#endif 1240 1241/* 1242 * Inner function for computing a point multiplication. A window is 1243 * provided, with points 1*P to 15*P in affine coordinates. 1244 * 1245 * Assumptions: 1246 * - All provided points are valid points on the curve. 1247 * - Multiplier is non-zero, and smaller than the curve order. 1248 * - Everything is in Montgomery representation. 1249 */ 1250static void 1251point_mul_inner(p256_jacobian *R, const p256_affine *W, 1252 const unsigned char *k, size_t klen) 1253{ 1254 p256_jacobian Q; 1255 uint32_t qz; 1256 1257 memset(&Q, 0, sizeof Q); 1258 qz = 1; 1259 while (klen -- > 0) { 1260 int i; 1261 unsigned bk; 1262 1263 bk = *k ++; 1264 for (i = 0; i < 2; i ++) { 1265 uint32_t bits; 1266 uint32_t bnz; 1267 p256_affine T; 1268 p256_jacobian U; 1269 uint32_t n; 1270 int j; 1271 uint64_t m; 1272 1273 p256_double(&Q); 1274 p256_double(&Q); 1275 p256_double(&Q); 1276 p256_double(&Q); 1277 bits = (bk >> 4) & 0x0F; 1278 bnz = NEQ(bits, 0); 1279 1280 /* 1281 * Lookup point in window. If the bits are 0, 1282 * we get something invalid, which is not a 1283 * problem because we will use it only if the 1284 * bits are non-zero. 1285 */ 1286 memset(&T, 0, sizeof T); 1287 for (n = 0; n < 15; n ++) { 1288 m = -(uint64_t)EQ(bits, n + 1); 1289 T.x[0] |= m & W[n].x[0]; 1290 T.x[1] |= m & W[n].x[1]; 1291 T.x[2] |= m & W[n].x[2]; 1292 T.x[3] |= m & W[n].x[3]; 1293 T.x[4] |= m & W[n].x[4]; 1294 T.y[0] |= m & W[n].y[0]; 1295 T.y[1] |= m & W[n].y[1]; 1296 T.y[2] |= m & W[n].y[2]; 1297 T.y[3] |= m & W[n].y[3]; 1298 T.y[4] |= m & W[n].y[4]; 1299 } 1300 1301 U = Q; 1302 p256_add_mixed(&U, &T); 1303 1304 /* 1305 * If qz is still 1, then Q was all-zeros, and this 1306 * is conserved through p256_double(). 1307 */ 1308 m = -(uint64_t)(bnz & qz); 1309 for (j = 0; j < 5; j ++) { 1310 Q.x[j] ^= m & (Q.x[j] ^ T.x[j]); 1311 Q.y[j] ^= m & (Q.y[j] ^ T.y[j]); 1312 Q.z[j] ^= m & (Q.z[j] ^ F256_R[j]); 1313 } 1314 CCOPY(bnz & ~qz, &Q, &U, sizeof Q); 1315 qz &= ~bnz; 1316 bk <<= 4; 1317 } 1318 } 1319 *R = Q; 1320} 1321 1322/* 1323 * Convert a window from Jacobian to affine coordinates. A single 1324 * field inversion is used. This function works for windows up to 1325 * 32 elements. 1326 * 1327 * The destination array (aff[]) and the source array (jac[]) may 1328 * overlap, provided that the start of aff[] is not after the start of 1329 * jac[]. Even if the arrays do _not_ overlap, the source array is 1330 * modified. 1331 */ 1332static void 1333window_to_affine(p256_affine *aff, p256_jacobian *jac, int num) 1334{ 1335 /* 1336 * Convert the window points to affine coordinates. We use the 1337 * following trick to mutualize the inversion computation: if 1338 * we have z1, z2, z3, and z4, and want to invert all of them, 1339 * we compute u = 1/(z1*z2*z3*z4), and then we have: 1340 * 1/z1 = u*z2*z3*z4 1341 * 1/z2 = u*z1*z3*z4 1342 * 1/z3 = u*z1*z2*z4 1343 * 1/z4 = u*z1*z2*z3 1344 * 1345 * The partial products are computed recursively: 1346 * 1347 * - on input (z_1,z_2), return (z_2,z_1) and z_1*z_2 1348 * - on input (z_1,z_2,... z_n): 1349 * recurse on (z_1,z_2,... z_(n/2)) -> r1 and m1 1350 * recurse on (z_(n/2+1),z_(n/2+2)... z_n) -> r2 and m2 1351 * multiply elements of r1 by m2 -> s1 1352 * multiply elements of r2 by m1 -> s2 1353 * return r1||r2 and m1*m2 1354 * 1355 * In the example below, we suppose that we have 14 elements. 1356 * Let z1, z2,... zE be the 14 values to invert (index noted in 1357 * hexadecimal, starting at 1). 1358 * 1359 * - Depth 1: 1360 * swap(z1, z2); z12 = z1*z2 1361 * swap(z3, z4); z34 = z3*z4 1362 * swap(z5, z6); z56 = z5*z6 1363 * swap(z7, z8); z78 = z7*z8 1364 * swap(z9, zA); z9A = z9*zA 1365 * swap(zB, zC); zBC = zB*zC 1366 * swap(zD, zE); zDE = zD*zE 1367 * 1368 * - Depth 2: 1369 * z1 <- z1*z34, z2 <- z2*z34, z3 <- z3*z12, z4 <- z4*z12 1370 * z1234 = z12*z34 1371 * z5 <- z5*z78, z6 <- z6*z78, z7 <- z7*z56, z8 <- z8*z56 1372 * z5678 = z56*z78 1373 * z9 <- z9*zBC, zA <- zA*zBC, zB <- zB*z9A, zC <- zC*z9A 1374 * z9ABC = z9A*zBC 1375 * 1376 * - Depth 3: 1377 * z1 <- z1*z5678, z2 <- z2*z5678, z3 <- z3*z5678, z4 <- z4*z5678 1378 * z5 <- z5*z1234, z6 <- z6*z1234, z7 <- z7*z1234, z8 <- z8*z1234 1379 * z12345678 = z1234*z5678 1380 * z9 <- z9*zDE, zA <- zA*zDE, zB <- zB*zDE, zC <- zC*zDE 1381 * zD <- zD*z9ABC, zE*z9ABC 1382 * z9ABCDE = z9ABC*zDE 1383 * 1384 * - Depth 4: 1385 * multiply z1..z8 by z9ABCDE 1386 * multiply z9..zE by z12345678 1387 * final z = z12345678*z9ABCDE 1388 */ 1389 1390 uint64_t z[16][5]; 1391 int i, k, s; 1392#define zt (z[15]) 1393#define zu (z[14]) 1394#define zv (z[13]) 1395 1396 /* 1397 * First recursion step (pairwise swapping and multiplication). 1398 * If there is an odd number of elements, then we "invent" an 1399 * extra one with coordinate Z = 1 (in Montgomery representation). 1400 */ 1401 for (i = 0; (i + 1) < num; i += 2) { 1402 memcpy(zt, jac[i].z, sizeof zt); 1403 memcpy(jac[i].z, jac[i + 1].z, sizeof zt); 1404 memcpy(jac[i + 1].z, zt, sizeof zt); 1405 f256_montymul(z[i >> 1], jac[i].z, jac[i + 1].z); 1406 } 1407 if ((num & 1) != 0) { 1408 memcpy(z[num >> 1], jac[num - 1].z, sizeof zt); 1409 memcpy(jac[num - 1].z, F256_R, sizeof F256_R); 1410 } 1411 1412 /* 1413 * Perform further recursion steps. At the entry of each step, 1414 * the process has been done for groups of 's' points. The 1415 * integer k is the log2 of s. 1416 */ 1417 for (k = 1, s = 2; s < num; k ++, s <<= 1) { 1418 int n; 1419 1420 for (i = 0; i < num; i ++) { 1421 f256_montymul(jac[i].z, jac[i].z, z[(i >> k) ^ 1]); 1422 } 1423 n = (num + s - 1) >> k; 1424 for (i = 0; i < (n >> 1); i ++) { 1425 f256_montymul(z[i], z[i << 1], z[(i << 1) + 1]); 1426 } 1427 if ((n & 1) != 0) { 1428 memmove(z[n >> 1], z[n], sizeof zt); 1429 } 1430 } 1431 1432 /* 1433 * Invert the final result, and convert all points. 1434 */ 1435 f256_invert(zt, z[0]); 1436 for (i = 0; i < num; i ++) { 1437 f256_montymul(zv, jac[i].z, zt); 1438 f256_montysquare(zu, zv); 1439 f256_montymul(zv, zv, zu); 1440 f256_montymul(aff[i].x, jac[i].x, zu); 1441 f256_montymul(aff[i].y, jac[i].y, zv); 1442 } 1443} 1444 1445/* 1446 * Multiply the provided point by an integer. 1447 * Assumptions: 1448 * - Source point is a valid curve point. 1449 * - Source point is not the point-at-infinity. 1450 * - Integer is not 0, and is lower than the curve order. 1451 * If these conditions are not met, then the result is indeterminate 1452 * (but the process is still constant-time). 1453 */ 1454static void 1455p256_mul(p256_jacobian *P, const unsigned char *k, size_t klen) 1456{ 1457 union { 1458 p256_affine aff[15]; 1459 p256_jacobian jac[15]; 1460 } window; 1461 int i; 1462 1463 /* 1464 * Compute window, in Jacobian coordinates. 1465 */ 1466 window.jac[0] = *P; 1467 for (i = 2; i < 16; i ++) { 1468 window.jac[i - 1] = window.jac[(i >> 1) - 1]; 1469 if ((i & 1) == 0) { 1470 p256_double(&window.jac[i - 1]); 1471 } else { 1472 p256_add(&window.jac[i - 1], &window.jac[i >> 1]); 1473 } 1474 } 1475 1476 /* 1477 * Convert the window points to affine coordinates. Point 1478 * window[0] is the source point, already in affine coordinates. 1479 */ 1480 window_to_affine(window.aff, window.jac, 15); 1481 1482 /* 1483 * Perform point multiplication. 1484 */ 1485 point_mul_inner(P, window.aff, k, klen); 1486} 1487 1488/* 1489 * Precomputed window for the conventional generator: P256_Gwin[n] 1490 * contains (n+1)*G (affine coordinates, in Montgomery representation). 1491 */ 1492static const p256_affine P256_Gwin[] = { 1493 { 1494 { 0x30D418A9143C1, 0xC4FEDB60179E7, 0x62251075BA95F, 1495 0x5C669FB732B77, 0x08905F76B5375 }, 1496 { 0x5357CE95560A8, 0x43A19E45CDDF2, 0x21F3258B4AB8E, 1497 0xD8552E88688DD, 0x0571FF18A5885 } 1498 }, 1499 { 1500 { 0x46D410DDD64DF, 0x0B433827D8500, 0x1490D9AA6AE3C, 1501 0xA3A832205038D, 0x06BB32E52DCF3 }, 1502 { 0x48D361BEE1A57, 0xB7B236FF82F36, 0x042DBE152CD7C, 1503 0xA3AA9A8FB0E92, 0x08C577517A5B8 } 1504 }, 1505 { 1506 { 0x3F904EEBC1272, 0x9E87D81FBFFAC, 0xCBBC98B027F84, 1507 0x47E46AD77DD87, 0x06936A3FD6FF7 }, 1508 { 0x5C1FC983A7EBD, 0xC3861FE1AB04C, 0x2EE98E583E47A, 1509 0xC06A88208311A, 0x05F06A2AB587C } 1510 }, 1511 { 1512 { 0xB50D46918DCC5, 0xD7623C17374B0, 0x100AF24650A6E, 1513 0x76ABCDAACACE8, 0x077362F591B01 }, 1514 { 0xF24CE4CBABA68, 0x17AD6F4472D96, 0xDDD22E1762847, 1515 0x862EB6C36DEE5, 0x04B14C39CC5AB } 1516 }, 1517 { 1518 { 0x8AAEC45C61F5C, 0x9D4B9537DBE1B, 0x76C20C90EC649, 1519 0x3C7D41CB5AAD0, 0x0907960649052 }, 1520 { 0x9B4AE7BA4F107, 0xF75EB882BEB30, 0x7A1F6873C568E, 1521 0x915C540A9877E, 0x03A076BB9DD1E } 1522 }, 1523 { 1524 { 0x47373E77664A1, 0xF246CEE3E4039, 0x17A3AD55AE744, 1525 0x673C50A961A5B, 0x03074B5964213 }, 1526 { 0x6220D377E44BA, 0x30DFF14B593D3, 0x639F11299C2B5, 1527 0x75F5424D44CEF, 0x04C9916DEA07F } 1528 }, 1529 { 1530 { 0x354EA0173B4F1, 0x3C23C00F70746, 0x23BB082BD2021, 1531 0xE03E43EAAB50C, 0x03BA5119D3123 }, 1532 { 0xD0303F5B9D4DE, 0x17DA67BDD2847, 0xC941956742F2F, 1533 0x8670F933BDC77, 0x0AEDD9164E240 } 1534 }, 1535 { 1536 { 0x4CD19499A78FB, 0x4BF9B345527F1, 0x2CFC6B462AB5C, 1537 0x30CDF90F02AF0, 0x0763891F62652 }, 1538 { 0xA3A9532D49775, 0xD7F9EBA15F59D, 0x60BBF021E3327, 1539 0xF75C23C7B84BE, 0x06EC12F2C706D } 1540 }, 1541 { 1542 { 0x6E8F264E20E8E, 0xC79A7A84175C9, 0xC8EB00ABE6BFE, 1543 0x16A4CC09C0444, 0x005B3081D0C4E }, 1544 { 0x777AA45F33140, 0xDCE5D45E31EB7, 0xB12F1A56AF7BE, 1545 0xF9B2B6E019A88, 0x086659CDFD835 } 1546 }, 1547 { 1548 { 0xDBD19DC21EC8C, 0x94FCF81392C18, 0x250B4998F9868, 1549 0x28EB37D2CD648, 0x0C61C947E4B34 }, 1550 { 0x407880DD9E767, 0x0C83FBE080C2B, 0x9BE5D2C43A899, 1551 0xAB4EF7D2D6577, 0x08719A555B3B4 } 1552 }, 1553 { 1554 { 0x260A6245E4043, 0x53E7FDFE0EA7D, 0xAC1AB59DE4079, 1555 0x072EFF3A4158D, 0x0E7090F1949C9 }, 1556 { 0x85612B944E886, 0xE857F61C81A76, 0xAD643D250F939, 1557 0x88DAC0DAA891E, 0x089300244125B } 1558 }, 1559 { 1560 { 0x1AA7D26977684, 0x58A345A3304B7, 0x37385EABDEDEF, 1561 0x155E409D29DEE, 0x0EE1DF780B83E }, 1562 { 0x12D91CBB5B437, 0x65A8956370CAC, 0xDE6D66170ED2F, 1563 0xAC9B8228CFA8A, 0x0FF57C95C3238 } 1564 }, 1565 { 1566 { 0x25634B2ED7097, 0x9156FD30DCCC4, 0x9E98110E35676, 1567 0x7594CBCD43F55, 0x038477ACC395B }, 1568 { 0x2B90C00EE17FF, 0xF842ED2E33575, 0x1F5BC16874838, 1569 0x7968CD06422BD, 0x0BC0876AB9E7B } 1570 }, 1571 { 1572 { 0xA35BB0CF664AF, 0x68F9707E3A242, 0x832660126E48F, 1573 0x72D2717BF54C6, 0x0AAE7333ED12C }, 1574 { 0x2DB7995D586B1, 0xE732237C227B5, 0x65E7DBBE29569, 1575 0xBBBD8E4193E2A, 0x052706DC3EAA1 } 1576 }, 1577 { 1578 { 0xD8B7BC60055BE, 0xD76E27E4B72BC, 0x81937003CC23E, 1579 0xA090E337424E4, 0x02AA0E43EAD3D }, 1580 { 0x524F6383C45D2, 0x422A41B2540B8, 0x8A4797D766355, 1581 0xDF444EFA6DE77, 0x0042170A9079A } 1582 }, 1583}; 1584 1585/* 1586 * Multiply the conventional generator of the curve by the provided 1587 * integer. Return is written in *P. 1588 * 1589 * Assumptions: 1590 * - Integer is not 0, and is lower than the curve order. 1591 * If this conditions is not met, then the result is indeterminate 1592 * (but the process is still constant-time). 1593 */ 1594static void 1595p256_mulgen(p256_jacobian *P, const unsigned char *k, size_t klen) 1596{ 1597 point_mul_inner(P, P256_Gwin, k, klen); 1598} 1599 1600/* 1601 * Return 1 if all of the following hold: 1602 * - klen <= 32 1603 * - k != 0 1604 * - k is lower than the curve order 1605 * Otherwise, return 0. 1606 * 1607 * Constant-time behaviour: only klen may be observable. 1608 */ 1609static uint32_t 1610check_scalar(const unsigned char *k, size_t klen) 1611{ 1612 uint32_t z; 1613 int32_t c; 1614 size_t u; 1615 1616 if (klen > 32) { 1617 return 0; 1618 } 1619 z = 0; 1620 for (u = 0; u < klen; u ++) { 1621 z |= k[u]; 1622 } 1623 if (klen == 32) { 1624 c = 0; 1625 for (u = 0; u < klen; u ++) { 1626 c |= -(int32_t)EQ0(c) & CMP(k[u], P256_N[u]); 1627 } 1628 } else { 1629 c = -1; 1630 } 1631 return NEQ(z, 0) & LT0(c); 1632} 1633 1634static uint32_t 1635api_mul(unsigned char *G, size_t Glen, 1636 const unsigned char *k, size_t klen, int curve) 1637{ 1638 uint32_t r; 1639 p256_jacobian P; 1640 1641 (void)curve; 1642 if (Glen != 65) { 1643 return 0; 1644 } 1645 r = check_scalar(k, klen); 1646 r &= point_decode(&P, G); 1647 p256_mul(&P, k, klen); 1648 r &= point_encode(G, &P); 1649 return r; 1650} 1651 1652static size_t 1653api_mulgen(unsigned char *R, 1654 const unsigned char *k, size_t klen, int curve) 1655{ 1656 p256_jacobian P; 1657 1658 (void)curve; 1659 p256_mulgen(&P, k, klen); 1660 point_encode(R, &P); 1661 return 65; 1662} 1663 1664static uint32_t 1665api_muladd(unsigned char *A, const unsigned char *B, size_t len, 1666 const unsigned char *x, size_t xlen, 1667 const unsigned char *y, size_t ylen, int curve) 1668{ 1669 /* 1670 * We might want to use Shamir's trick here: make a composite 1671 * window of u*P+v*Q points, to merge the two doubling-ladders 1672 * into one. This, however, has some complications: 1673 * 1674 * - During the computation, we may hit the point-at-infinity. 1675 * Thus, we would need p256_add_complete_mixed() (complete 1676 * formulas for point addition), with a higher cost (17 muls 1677 * instead of 11). 1678 * 1679 * - A 4-bit window would be too large, since it would involve 1680 * 16*16-1 = 255 points. For the same window size as in the 1681 * p256_mul() case, we would need to reduce the window size 1682 * to 2 bits, and thus perform twice as many non-doubling 1683 * point additions. 1684 * 1685 * - The window may itself contain the point-at-infinity, and 1686 * thus cannot be in all generality be made of affine points. 1687 * Instead, we would need to make it a window of points in 1688 * Jacobian coordinates. Even p256_add_complete_mixed() would 1689 * be inappropriate. 1690 * 1691 * For these reasons, the code below performs two separate 1692 * point multiplications, then computes the final point addition 1693 * (which is both a "normal" addition, and a doubling, to handle 1694 * all cases). 1695 */ 1696 1697 p256_jacobian P, Q; 1698 uint32_t r, t, s; 1699 uint64_t z; 1700 1701 (void)curve; 1702 if (len != 65) { 1703 return 0; 1704 } 1705 r = point_decode(&P, A); 1706 p256_mul(&P, x, xlen); 1707 if (B == NULL) { 1708 p256_mulgen(&Q, y, ylen); 1709 } else { 1710 r &= point_decode(&Q, B); 1711 p256_mul(&Q, y, ylen); 1712 } 1713 1714 /* 1715 * The final addition may fail in case both points are equal. 1716 */ 1717 t = p256_add(&P, &Q); 1718 f256_final_reduce(P.z); 1719 z = P.z[0] | P.z[1] | P.z[2] | P.z[3] | P.z[4]; 1720 s = EQ((uint32_t)(z | (z >> 32)), 0); 1721 p256_double(&Q); 1722 1723 /* 1724 * If s is 1 then either P+Q = 0 (t = 1) or P = Q (t = 0). So we 1725 * have the following: 1726 * 1727 * s = 0, t = 0 return P (normal addition) 1728 * s = 0, t = 1 return P (normal addition) 1729 * s = 1, t = 0 return Q (a 'double' case) 1730 * s = 1, t = 1 report an error (P+Q = 0) 1731 */ 1732 CCOPY(s & ~t, &P, &Q, sizeof Q); 1733 point_encode(A, &P); 1734 r &= ~(s & t); 1735 return r; 1736} 1737 1738/* see bearssl_ec.h */ 1739const br_ec_impl br_ec_p256_m62 = { 1740 (uint32_t)0x00800000, 1741 &api_generator, 1742 &api_order, 1743 &api_xoff, 1744 &api_mul, 1745 &api_mulgen, 1746 &api_muladd 1747}; 1748 1749/* see bearssl_ec.h */ 1750const br_ec_impl * 1751br_ec_p256_m62_get(void) 1752{ 1753 return &br_ec_p256_m62; 1754} 1755 1756#else 1757 1758/* see bearssl_ec.h */ 1759const br_ec_impl * 1760br_ec_p256_m62_get(void) 1761{ 1762 return 0; 1763} 1764 1765#endif 1766