1 /* 2 * Copyright 2011-2019 The OpenSSL Project Authors. All Rights Reserved. 3 * 4 * Licensed under the OpenSSL license (the "License"). You may not use 5 * this file except in compliance with the License. You can obtain a copy 6 * in the file LICENSE in the source distribution or at 7 * https://www.openssl.org/source/license.html 8 */ 9 10 /* Copyright 2011 Google Inc. 11 * 12 * Licensed under the Apache License, Version 2.0 (the "License"); 13 * 14 * you may not use this file except in compliance with the License. 15 * You may obtain a copy of the License at 16 * 17 * http://www.apache.org/licenses/LICENSE-2.0 18 * 19 * Unless required by applicable law or agreed to in writing, software 20 * distributed under the License is distributed on an "AS IS" BASIS, 21 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 22 * See the License for the specific language governing permissions and 23 * limitations under the License. 24 */ 25 26 /* 27 * A 64-bit implementation of the NIST P-521 elliptic curve point multiplication 28 * 29 * OpenSSL integration was taken from Emilia Kasper's work in ecp_nistp224.c. 30 * Otherwise based on Emilia's P224 work, which was inspired by my curve25519 31 * work which got its smarts from Daniel J. Bernstein's work on the same. 32 */ 33 34 #include <openssl/e_os2.h> 35 #ifdef OPENSSL_NO_EC_NISTP_64_GCC_128 36 NON_EMPTY_TRANSLATION_UNIT 37 #else 38 39 # include <string.h> 40 # include <openssl/err.h> 41 # include "ec_lcl.h" 42 43 # if defined(__SIZEOF_INT128__) && __SIZEOF_INT128__==16 44 /* even with gcc, the typedef won't work for 32-bit platforms */ 45 typedef __uint128_t uint128_t; /* nonstandard; implemented by gcc on 64-bit 46 * platforms */ 47 # else 48 # error "Your compiler doesn't appear to support 128-bit integer types" 49 # endif 50 51 typedef uint8_t u8; 52 typedef uint64_t u64; 53 54 /* 55 * The underlying field. P521 operates over GF(2^521-1). We can serialise an 56 * element of this field into 66 bytes where the most significant byte 57 * contains only a single bit. We call this an felem_bytearray. 58 */ 59 60 typedef u8 felem_bytearray[66]; 61 62 /* 63 * These are the parameters of P521, taken from FIPS 186-3, section D.1.2.5. 64 * These values are big-endian. 65 */ 66 static const felem_bytearray nistp521_curve_params[5] = { 67 {0x01, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, /* p */ 68 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 69 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 70 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 71 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 72 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 73 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 74 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 75 0xff, 0xff}, 76 {0x01, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, /* a = -3 */ 77 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 78 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 79 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 80 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 81 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 82 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 83 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 84 0xff, 0xfc}, 85 {0x00, 0x51, 0x95, 0x3e, 0xb9, 0x61, 0x8e, 0x1c, /* b */ 86 0x9a, 0x1f, 0x92, 0x9a, 0x21, 0xa0, 0xb6, 0x85, 87 0x40, 0xee, 0xa2, 0xda, 0x72, 0x5b, 0x99, 0xb3, 88 0x15, 0xf3, 0xb8, 0xb4, 0x89, 0x91, 0x8e, 0xf1, 89 0x09, 0xe1, 0x56, 0x19, 0x39, 0x51, 0xec, 0x7e, 90 0x93, 0x7b, 0x16, 0x52, 0xc0, 0xbd, 0x3b, 0xb1, 91 0xbf, 0x07, 0x35, 0x73, 0xdf, 0x88, 0x3d, 0x2c, 92 0x34, 0xf1, 0xef, 0x45, 0x1f, 0xd4, 0x6b, 0x50, 93 0x3f, 0x00}, 94 {0x00, 0xc6, 0x85, 0x8e, 0x06, 0xb7, 0x04, 0x04, /* x */ 95 0xe9, 0xcd, 0x9e, 0x3e, 0xcb, 0x66, 0x23, 0x95, 96 0xb4, 0x42, 0x9c, 0x64, 0x81, 0x39, 0x05, 0x3f, 97 0xb5, 0x21, 0xf8, 0x28, 0xaf, 0x60, 0x6b, 0x4d, 98 0x3d, 0xba, 0xa1, 0x4b, 0x5e, 0x77, 0xef, 0xe7, 99 0x59, 0x28, 0xfe, 0x1d, 0xc1, 0x27, 0xa2, 0xff, 100 0xa8, 0xde, 0x33, 0x48, 0xb3, 0xc1, 0x85, 0x6a, 101 0x42, 0x9b, 0xf9, 0x7e, 0x7e, 0x31, 0xc2, 0xe5, 102 0xbd, 0x66}, 103 {0x01, 0x18, 0x39, 0x29, 0x6a, 0x78, 0x9a, 0x3b, /* y */ 104 0xc0, 0x04, 0x5c, 0x8a, 0x5f, 0xb4, 0x2c, 0x7d, 105 0x1b, 0xd9, 0x98, 0xf5, 0x44, 0x49, 0x57, 0x9b, 106 0x44, 0x68, 0x17, 0xaf, 0xbd, 0x17, 0x27, 0x3e, 107 0x66, 0x2c, 0x97, 0xee, 0x72, 0x99, 0x5e, 0xf4, 108 0x26, 0x40, 0xc5, 0x50, 0xb9, 0x01, 0x3f, 0xad, 109 0x07, 0x61, 0x35, 0x3c, 0x70, 0x86, 0xa2, 0x72, 110 0xc2, 0x40, 0x88, 0xbe, 0x94, 0x76, 0x9f, 0xd1, 111 0x66, 0x50} 112 }; 113 114 /*- 115 * The representation of field elements. 116 * ------------------------------------ 117 * 118 * We represent field elements with nine values. These values are either 64 or 119 * 128 bits and the field element represented is: 120 * v[0]*2^0 + v[1]*2^58 + v[2]*2^116 + ... + v[8]*2^464 (mod p) 121 * Each of the nine values is called a 'limb'. Since the limbs are spaced only 122 * 58 bits apart, but are greater than 58 bits in length, the most significant 123 * bits of each limb overlap with the least significant bits of the next. 124 * 125 * A field element with 64-bit limbs is an 'felem'. One with 128-bit limbs is a 126 * 'largefelem' */ 127 128 # define NLIMBS 9 129 130 typedef uint64_t limb; 131 typedef limb felem[NLIMBS]; 132 typedef uint128_t largefelem[NLIMBS]; 133 134 static const limb bottom57bits = 0x1ffffffffffffff; 135 static const limb bottom58bits = 0x3ffffffffffffff; 136 137 /* 138 * bin66_to_felem takes a little-endian byte array and converts it into felem 139 * form. This assumes that the CPU is little-endian. 140 */ 141 static void bin66_to_felem(felem out, const u8 in[66]) 142 { 143 out[0] = (*((limb *) & in[0])) & bottom58bits; 144 out[1] = (*((limb *) & in[7]) >> 2) & bottom58bits; 145 out[2] = (*((limb *) & in[14]) >> 4) & bottom58bits; 146 out[3] = (*((limb *) & in[21]) >> 6) & bottom58bits; 147 out[4] = (*((limb *) & in[29])) & bottom58bits; 148 out[5] = (*((limb *) & in[36]) >> 2) & bottom58bits; 149 out[6] = (*((limb *) & in[43]) >> 4) & bottom58bits; 150 out[7] = (*((limb *) & in[50]) >> 6) & bottom58bits; 151 out[8] = (*((limb *) & in[58])) & bottom57bits; 152 } 153 154 /* 155 * felem_to_bin66 takes an felem and serialises into a little endian, 66 byte 156 * array. This assumes that the CPU is little-endian. 157 */ 158 static void felem_to_bin66(u8 out[66], const felem in) 159 { 160 memset(out, 0, 66); 161 (*((limb *) & out[0])) = in[0]; 162 (*((limb *) & out[7])) |= in[1] << 2; 163 (*((limb *) & out[14])) |= in[2] << 4; 164 (*((limb *) & out[21])) |= in[3] << 6; 165 (*((limb *) & out[29])) = in[4]; 166 (*((limb *) & out[36])) |= in[5] << 2; 167 (*((limb *) & out[43])) |= in[6] << 4; 168 (*((limb *) & out[50])) |= in[7] << 6; 169 (*((limb *) & out[58])) = in[8]; 170 } 171 172 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */ 173 static void flip_endian(u8 *out, const u8 *in, unsigned len) 174 { 175 unsigned i; 176 for (i = 0; i < len; ++i) 177 out[i] = in[len - 1 - i]; 178 } 179 180 /* BN_to_felem converts an OpenSSL BIGNUM into an felem */ 181 static int BN_to_felem(felem out, const BIGNUM *bn) 182 { 183 felem_bytearray b_in; 184 felem_bytearray b_out; 185 unsigned num_bytes; 186 187 /* BN_bn2bin eats leading zeroes */ 188 memset(b_out, 0, sizeof(b_out)); 189 num_bytes = BN_num_bytes(bn); 190 if (num_bytes > sizeof(b_out)) { 191 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE); 192 return 0; 193 } 194 if (BN_is_negative(bn)) { 195 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE); 196 return 0; 197 } 198 num_bytes = BN_bn2bin(bn, b_in); 199 flip_endian(b_out, b_in, num_bytes); 200 bin66_to_felem(out, b_out); 201 return 1; 202 } 203 204 /* felem_to_BN converts an felem into an OpenSSL BIGNUM */ 205 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in) 206 { 207 felem_bytearray b_in, b_out; 208 felem_to_bin66(b_in, in); 209 flip_endian(b_out, b_in, sizeof(b_out)); 210 return BN_bin2bn(b_out, sizeof(b_out), out); 211 } 212 213 /*- 214 * Field operations 215 * ---------------- 216 */ 217 218 static void felem_one(felem out) 219 { 220 out[0] = 1; 221 out[1] = 0; 222 out[2] = 0; 223 out[3] = 0; 224 out[4] = 0; 225 out[5] = 0; 226 out[6] = 0; 227 out[7] = 0; 228 out[8] = 0; 229 } 230 231 static void felem_assign(felem out, const felem in) 232 { 233 out[0] = in[0]; 234 out[1] = in[1]; 235 out[2] = in[2]; 236 out[3] = in[3]; 237 out[4] = in[4]; 238 out[5] = in[5]; 239 out[6] = in[6]; 240 out[7] = in[7]; 241 out[8] = in[8]; 242 } 243 244 /* felem_sum64 sets out = out + in. */ 245 static void felem_sum64(felem out, const felem in) 246 { 247 out[0] += in[0]; 248 out[1] += in[1]; 249 out[2] += in[2]; 250 out[3] += in[3]; 251 out[4] += in[4]; 252 out[5] += in[5]; 253 out[6] += in[6]; 254 out[7] += in[7]; 255 out[8] += in[8]; 256 } 257 258 /* felem_scalar sets out = in * scalar */ 259 static void felem_scalar(felem out, const felem in, limb scalar) 260 { 261 out[0] = in[0] * scalar; 262 out[1] = in[1] * scalar; 263 out[2] = in[2] * scalar; 264 out[3] = in[3] * scalar; 265 out[4] = in[4] * scalar; 266 out[5] = in[5] * scalar; 267 out[6] = in[6] * scalar; 268 out[7] = in[7] * scalar; 269 out[8] = in[8] * scalar; 270 } 271 272 /* felem_scalar64 sets out = out * scalar */ 273 static void felem_scalar64(felem out, limb scalar) 274 { 275 out[0] *= scalar; 276 out[1] *= scalar; 277 out[2] *= scalar; 278 out[3] *= scalar; 279 out[4] *= scalar; 280 out[5] *= scalar; 281 out[6] *= scalar; 282 out[7] *= scalar; 283 out[8] *= scalar; 284 } 285 286 /* felem_scalar128 sets out = out * scalar */ 287 static void felem_scalar128(largefelem out, limb scalar) 288 { 289 out[0] *= scalar; 290 out[1] *= scalar; 291 out[2] *= scalar; 292 out[3] *= scalar; 293 out[4] *= scalar; 294 out[5] *= scalar; 295 out[6] *= scalar; 296 out[7] *= scalar; 297 out[8] *= scalar; 298 } 299 300 /*- 301 * felem_neg sets |out| to |-in| 302 * On entry: 303 * in[i] < 2^59 + 2^14 304 * On exit: 305 * out[i] < 2^62 306 */ 307 static void felem_neg(felem out, const felem in) 308 { 309 /* In order to prevent underflow, we subtract from 0 mod p. */ 310 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5); 311 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4); 312 313 out[0] = two62m3 - in[0]; 314 out[1] = two62m2 - in[1]; 315 out[2] = two62m2 - in[2]; 316 out[3] = two62m2 - in[3]; 317 out[4] = two62m2 - in[4]; 318 out[5] = two62m2 - in[5]; 319 out[6] = two62m2 - in[6]; 320 out[7] = two62m2 - in[7]; 321 out[8] = two62m2 - in[8]; 322 } 323 324 /*- 325 * felem_diff64 subtracts |in| from |out| 326 * On entry: 327 * in[i] < 2^59 + 2^14 328 * On exit: 329 * out[i] < out[i] + 2^62 330 */ 331 static void felem_diff64(felem out, const felem in) 332 { 333 /* 334 * In order to prevent underflow, we add 0 mod p before subtracting. 335 */ 336 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5); 337 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4); 338 339 out[0] += two62m3 - in[0]; 340 out[1] += two62m2 - in[1]; 341 out[2] += two62m2 - in[2]; 342 out[3] += two62m2 - in[3]; 343 out[4] += two62m2 - in[4]; 344 out[5] += two62m2 - in[5]; 345 out[6] += two62m2 - in[6]; 346 out[7] += two62m2 - in[7]; 347 out[8] += two62m2 - in[8]; 348 } 349 350 /*- 351 * felem_diff_128_64 subtracts |in| from |out| 352 * On entry: 353 * in[i] < 2^62 + 2^17 354 * On exit: 355 * out[i] < out[i] + 2^63 356 */ 357 static void felem_diff_128_64(largefelem out, const felem in) 358 { 359 /* 360 * In order to prevent underflow, we add 64p mod p (which is equivalent 361 * to 0 mod p) before subtracting. p is 2^521 - 1, i.e. in binary a 521 362 * digit number with all bits set to 1. See "The representation of field 363 * elements" comment above for a description of how limbs are used to 364 * represent a number. 64p is represented with 8 limbs containing a number 365 * with 58 bits set and one limb with a number with 57 bits set. 366 */ 367 static const limb two63m6 = (((limb) 1) << 63) - (((limb) 1) << 6); 368 static const limb two63m5 = (((limb) 1) << 63) - (((limb) 1) << 5); 369 370 out[0] += two63m6 - in[0]; 371 out[1] += two63m5 - in[1]; 372 out[2] += two63m5 - in[2]; 373 out[3] += two63m5 - in[3]; 374 out[4] += two63m5 - in[4]; 375 out[5] += two63m5 - in[5]; 376 out[6] += two63m5 - in[6]; 377 out[7] += two63m5 - in[7]; 378 out[8] += two63m5 - in[8]; 379 } 380 381 /*- 382 * felem_diff_128_64 subtracts |in| from |out| 383 * On entry: 384 * in[i] < 2^126 385 * On exit: 386 * out[i] < out[i] + 2^127 - 2^69 387 */ 388 static void felem_diff128(largefelem out, const largefelem in) 389 { 390 /* 391 * In order to prevent underflow, we add 0 mod p before subtracting. 392 */ 393 static const uint128_t two127m70 = 394 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 70); 395 static const uint128_t two127m69 = 396 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 69); 397 398 out[0] += (two127m70 - in[0]); 399 out[1] += (two127m69 - in[1]); 400 out[2] += (two127m69 - in[2]); 401 out[3] += (two127m69 - in[3]); 402 out[4] += (two127m69 - in[4]); 403 out[5] += (two127m69 - in[5]); 404 out[6] += (two127m69 - in[6]); 405 out[7] += (two127m69 - in[7]); 406 out[8] += (two127m69 - in[8]); 407 } 408 409 /*- 410 * felem_square sets |out| = |in|^2 411 * On entry: 412 * in[i] < 2^62 413 * On exit: 414 * out[i] < 17 * max(in[i]) * max(in[i]) 415 */ 416 static void felem_square(largefelem out, const felem in) 417 { 418 felem inx2, inx4; 419 felem_scalar(inx2, in, 2); 420 felem_scalar(inx4, in, 4); 421 422 /*- 423 * We have many cases were we want to do 424 * in[x] * in[y] + 425 * in[y] * in[x] 426 * This is obviously just 427 * 2 * in[x] * in[y] 428 * However, rather than do the doubling on the 128 bit result, we 429 * double one of the inputs to the multiplication by reading from 430 * |inx2| 431 */ 432 433 out[0] = ((uint128_t) in[0]) * in[0]; 434 out[1] = ((uint128_t) in[0]) * inx2[1]; 435 out[2] = ((uint128_t) in[0]) * inx2[2] + ((uint128_t) in[1]) * in[1]; 436 out[3] = ((uint128_t) in[0]) * inx2[3] + ((uint128_t) in[1]) * inx2[2]; 437 out[4] = ((uint128_t) in[0]) * inx2[4] + 438 ((uint128_t) in[1]) * inx2[3] + ((uint128_t) in[2]) * in[2]; 439 out[5] = ((uint128_t) in[0]) * inx2[5] + 440 ((uint128_t) in[1]) * inx2[4] + ((uint128_t) in[2]) * inx2[3]; 441 out[6] = ((uint128_t) in[0]) * inx2[6] + 442 ((uint128_t) in[1]) * inx2[5] + 443 ((uint128_t) in[2]) * inx2[4] + ((uint128_t) in[3]) * in[3]; 444 out[7] = ((uint128_t) in[0]) * inx2[7] + 445 ((uint128_t) in[1]) * inx2[6] + 446 ((uint128_t) in[2]) * inx2[5] + ((uint128_t) in[3]) * inx2[4]; 447 out[8] = ((uint128_t) in[0]) * inx2[8] + 448 ((uint128_t) in[1]) * inx2[7] + 449 ((uint128_t) in[2]) * inx2[6] + 450 ((uint128_t) in[3]) * inx2[5] + ((uint128_t) in[4]) * in[4]; 451 452 /* 453 * The remaining limbs fall above 2^521, with the first falling at 2^522. 454 * They correspond to locations one bit up from the limbs produced above 455 * so we would have to multiply by two to align them. Again, rather than 456 * operate on the 128-bit result, we double one of the inputs to the 457 * multiplication. If we want to double for both this reason, and the 458 * reason above, then we end up multiplying by four. 459 */ 460 461 /* 9 */ 462 out[0] += ((uint128_t) in[1]) * inx4[8] + 463 ((uint128_t) in[2]) * inx4[7] + 464 ((uint128_t) in[3]) * inx4[6] + ((uint128_t) in[4]) * inx4[5]; 465 466 /* 10 */ 467 out[1] += ((uint128_t) in[2]) * inx4[8] + 468 ((uint128_t) in[3]) * inx4[7] + 469 ((uint128_t) in[4]) * inx4[6] + ((uint128_t) in[5]) * inx2[5]; 470 471 /* 11 */ 472 out[2] += ((uint128_t) in[3]) * inx4[8] + 473 ((uint128_t) in[4]) * inx4[7] + ((uint128_t) in[5]) * inx4[6]; 474 475 /* 12 */ 476 out[3] += ((uint128_t) in[4]) * inx4[8] + 477 ((uint128_t) in[5]) * inx4[7] + ((uint128_t) in[6]) * inx2[6]; 478 479 /* 13 */ 480 out[4] += ((uint128_t) in[5]) * inx4[8] + ((uint128_t) in[6]) * inx4[7]; 481 482 /* 14 */ 483 out[5] += ((uint128_t) in[6]) * inx4[8] + ((uint128_t) in[7]) * inx2[7]; 484 485 /* 15 */ 486 out[6] += ((uint128_t) in[7]) * inx4[8]; 487 488 /* 16 */ 489 out[7] += ((uint128_t) in[8]) * inx2[8]; 490 } 491 492 /*- 493 * felem_mul sets |out| = |in1| * |in2| 494 * On entry: 495 * in1[i] < 2^64 496 * in2[i] < 2^63 497 * On exit: 498 * out[i] < 17 * max(in1[i]) * max(in2[i]) 499 */ 500 static void felem_mul(largefelem out, const felem in1, const felem in2) 501 { 502 felem in2x2; 503 felem_scalar(in2x2, in2, 2); 504 505 out[0] = ((uint128_t) in1[0]) * in2[0]; 506 507 out[1] = ((uint128_t) in1[0]) * in2[1] + 508 ((uint128_t) in1[1]) * in2[0]; 509 510 out[2] = ((uint128_t) in1[0]) * in2[2] + 511 ((uint128_t) in1[1]) * in2[1] + 512 ((uint128_t) in1[2]) * in2[0]; 513 514 out[3] = ((uint128_t) in1[0]) * in2[3] + 515 ((uint128_t) in1[1]) * in2[2] + 516 ((uint128_t) in1[2]) * in2[1] + 517 ((uint128_t) in1[3]) * in2[0]; 518 519 out[4] = ((uint128_t) in1[0]) * in2[4] + 520 ((uint128_t) in1[1]) * in2[3] + 521 ((uint128_t) in1[2]) * in2[2] + 522 ((uint128_t) in1[3]) * in2[1] + 523 ((uint128_t) in1[4]) * in2[0]; 524 525 out[5] = ((uint128_t) in1[0]) * in2[5] + 526 ((uint128_t) in1[1]) * in2[4] + 527 ((uint128_t) in1[2]) * in2[3] + 528 ((uint128_t) in1[3]) * in2[2] + 529 ((uint128_t) in1[4]) * in2[1] + 530 ((uint128_t) in1[5]) * in2[0]; 531 532 out[6] = ((uint128_t) in1[0]) * in2[6] + 533 ((uint128_t) in1[1]) * in2[5] + 534 ((uint128_t) in1[2]) * in2[4] + 535 ((uint128_t) in1[3]) * in2[3] + 536 ((uint128_t) in1[4]) * in2[2] + 537 ((uint128_t) in1[5]) * in2[1] + 538 ((uint128_t) in1[6]) * in2[0]; 539 540 out[7] = ((uint128_t) in1[0]) * in2[7] + 541 ((uint128_t) in1[1]) * in2[6] + 542 ((uint128_t) in1[2]) * in2[5] + 543 ((uint128_t) in1[3]) * in2[4] + 544 ((uint128_t) in1[4]) * in2[3] + 545 ((uint128_t) in1[5]) * in2[2] + 546 ((uint128_t) in1[6]) * in2[1] + 547 ((uint128_t) in1[7]) * in2[0]; 548 549 out[8] = ((uint128_t) in1[0]) * in2[8] + 550 ((uint128_t) in1[1]) * in2[7] + 551 ((uint128_t) in1[2]) * in2[6] + 552 ((uint128_t) in1[3]) * in2[5] + 553 ((uint128_t) in1[4]) * in2[4] + 554 ((uint128_t) in1[5]) * in2[3] + 555 ((uint128_t) in1[6]) * in2[2] + 556 ((uint128_t) in1[7]) * in2[1] + 557 ((uint128_t) in1[8]) * in2[0]; 558 559 /* See comment in felem_square about the use of in2x2 here */ 560 561 out[0] += ((uint128_t) in1[1]) * in2x2[8] + 562 ((uint128_t) in1[2]) * in2x2[7] + 563 ((uint128_t) in1[3]) * in2x2[6] + 564 ((uint128_t) in1[4]) * in2x2[5] + 565 ((uint128_t) in1[5]) * in2x2[4] + 566 ((uint128_t) in1[6]) * in2x2[3] + 567 ((uint128_t) in1[7]) * in2x2[2] + 568 ((uint128_t) in1[8]) * in2x2[1]; 569 570 out[1] += ((uint128_t) in1[2]) * in2x2[8] + 571 ((uint128_t) in1[3]) * in2x2[7] + 572 ((uint128_t) in1[4]) * in2x2[6] + 573 ((uint128_t) in1[5]) * in2x2[5] + 574 ((uint128_t) in1[6]) * in2x2[4] + 575 ((uint128_t) in1[7]) * in2x2[3] + 576 ((uint128_t) in1[8]) * in2x2[2]; 577 578 out[2] += ((uint128_t) in1[3]) * in2x2[8] + 579 ((uint128_t) in1[4]) * in2x2[7] + 580 ((uint128_t) in1[5]) * in2x2[6] + 581 ((uint128_t) in1[6]) * in2x2[5] + 582 ((uint128_t) in1[7]) * in2x2[4] + 583 ((uint128_t) in1[8]) * in2x2[3]; 584 585 out[3] += ((uint128_t) in1[4]) * in2x2[8] + 586 ((uint128_t) in1[5]) * in2x2[7] + 587 ((uint128_t) in1[6]) * in2x2[6] + 588 ((uint128_t) in1[7]) * in2x2[5] + 589 ((uint128_t) in1[8]) * in2x2[4]; 590 591 out[4] += ((uint128_t) in1[5]) * in2x2[8] + 592 ((uint128_t) in1[6]) * in2x2[7] + 593 ((uint128_t) in1[7]) * in2x2[6] + 594 ((uint128_t) in1[8]) * in2x2[5]; 595 596 out[5] += ((uint128_t) in1[6]) * in2x2[8] + 597 ((uint128_t) in1[7]) * in2x2[7] + 598 ((uint128_t) in1[8]) * in2x2[6]; 599 600 out[6] += ((uint128_t) in1[7]) * in2x2[8] + 601 ((uint128_t) in1[8]) * in2x2[7]; 602 603 out[7] += ((uint128_t) in1[8]) * in2x2[8]; 604 } 605 606 static const limb bottom52bits = 0xfffffffffffff; 607 608 /*- 609 * felem_reduce converts a largefelem to an felem. 610 * On entry: 611 * in[i] < 2^128 612 * On exit: 613 * out[i] < 2^59 + 2^14 614 */ 615 static void felem_reduce(felem out, const largefelem in) 616 { 617 u64 overflow1, overflow2; 618 619 out[0] = ((limb) in[0]) & bottom58bits; 620 out[1] = ((limb) in[1]) & bottom58bits; 621 out[2] = ((limb) in[2]) & bottom58bits; 622 out[3] = ((limb) in[3]) & bottom58bits; 623 out[4] = ((limb) in[4]) & bottom58bits; 624 out[5] = ((limb) in[5]) & bottom58bits; 625 out[6] = ((limb) in[6]) & bottom58bits; 626 out[7] = ((limb) in[7]) & bottom58bits; 627 out[8] = ((limb) in[8]) & bottom58bits; 628 629 /* out[i] < 2^58 */ 630 631 out[1] += ((limb) in[0]) >> 58; 632 out[1] += (((limb) (in[0] >> 64)) & bottom52bits) << 6; 633 /*- 634 * out[1] < 2^58 + 2^6 + 2^58 635 * = 2^59 + 2^6 636 */ 637 out[2] += ((limb) (in[0] >> 64)) >> 52; 638 639 out[2] += ((limb) in[1]) >> 58; 640 out[2] += (((limb) (in[1] >> 64)) & bottom52bits) << 6; 641 out[3] += ((limb) (in[1] >> 64)) >> 52; 642 643 out[3] += ((limb) in[2]) >> 58; 644 out[3] += (((limb) (in[2] >> 64)) & bottom52bits) << 6; 645 out[4] += ((limb) (in[2] >> 64)) >> 52; 646 647 out[4] += ((limb) in[3]) >> 58; 648 out[4] += (((limb) (in[3] >> 64)) & bottom52bits) << 6; 649 out[5] += ((limb) (in[3] >> 64)) >> 52; 650 651 out[5] += ((limb) in[4]) >> 58; 652 out[5] += (((limb) (in[4] >> 64)) & bottom52bits) << 6; 653 out[6] += ((limb) (in[4] >> 64)) >> 52; 654 655 out[6] += ((limb) in[5]) >> 58; 656 out[6] += (((limb) (in[5] >> 64)) & bottom52bits) << 6; 657 out[7] += ((limb) (in[5] >> 64)) >> 52; 658 659 out[7] += ((limb) in[6]) >> 58; 660 out[7] += (((limb) (in[6] >> 64)) & bottom52bits) << 6; 661 out[8] += ((limb) (in[6] >> 64)) >> 52; 662 663 out[8] += ((limb) in[7]) >> 58; 664 out[8] += (((limb) (in[7] >> 64)) & bottom52bits) << 6; 665 /*- 666 * out[x > 1] < 2^58 + 2^6 + 2^58 + 2^12 667 * < 2^59 + 2^13 668 */ 669 overflow1 = ((limb) (in[7] >> 64)) >> 52; 670 671 overflow1 += ((limb) in[8]) >> 58; 672 overflow1 += (((limb) (in[8] >> 64)) & bottom52bits) << 6; 673 overflow2 = ((limb) (in[8] >> 64)) >> 52; 674 675 overflow1 <<= 1; /* overflow1 < 2^13 + 2^7 + 2^59 */ 676 overflow2 <<= 1; /* overflow2 < 2^13 */ 677 678 out[0] += overflow1; /* out[0] < 2^60 */ 679 out[1] += overflow2; /* out[1] < 2^59 + 2^6 + 2^13 */ 680 681 out[1] += out[0] >> 58; 682 out[0] &= bottom58bits; 683 /*- 684 * out[0] < 2^58 685 * out[1] < 2^59 + 2^6 + 2^13 + 2^2 686 * < 2^59 + 2^14 687 */ 688 } 689 690 static void felem_square_reduce(felem out, const felem in) 691 { 692 largefelem tmp; 693 felem_square(tmp, in); 694 felem_reduce(out, tmp); 695 } 696 697 static void felem_mul_reduce(felem out, const felem in1, const felem in2) 698 { 699 largefelem tmp; 700 felem_mul(tmp, in1, in2); 701 felem_reduce(out, tmp); 702 } 703 704 /*- 705 * felem_inv calculates |out| = |in|^{-1} 706 * 707 * Based on Fermat's Little Theorem: 708 * a^p = a (mod p) 709 * a^{p-1} = 1 (mod p) 710 * a^{p-2} = a^{-1} (mod p) 711 */ 712 static void felem_inv(felem out, const felem in) 713 { 714 felem ftmp, ftmp2, ftmp3, ftmp4; 715 largefelem tmp; 716 unsigned i; 717 718 felem_square(tmp, in); 719 felem_reduce(ftmp, tmp); /* 2^1 */ 720 felem_mul(tmp, in, ftmp); 721 felem_reduce(ftmp, tmp); /* 2^2 - 2^0 */ 722 felem_assign(ftmp2, ftmp); 723 felem_square(tmp, ftmp); 724 felem_reduce(ftmp, tmp); /* 2^3 - 2^1 */ 725 felem_mul(tmp, in, ftmp); 726 felem_reduce(ftmp, tmp); /* 2^3 - 2^0 */ 727 felem_square(tmp, ftmp); 728 felem_reduce(ftmp, tmp); /* 2^4 - 2^1 */ 729 730 felem_square(tmp, ftmp2); 731 felem_reduce(ftmp3, tmp); /* 2^3 - 2^1 */ 732 felem_square(tmp, ftmp3); 733 felem_reduce(ftmp3, tmp); /* 2^4 - 2^2 */ 734 felem_mul(tmp, ftmp3, ftmp2); 735 felem_reduce(ftmp3, tmp); /* 2^4 - 2^0 */ 736 737 felem_assign(ftmp2, ftmp3); 738 felem_square(tmp, ftmp3); 739 felem_reduce(ftmp3, tmp); /* 2^5 - 2^1 */ 740 felem_square(tmp, ftmp3); 741 felem_reduce(ftmp3, tmp); /* 2^6 - 2^2 */ 742 felem_square(tmp, ftmp3); 743 felem_reduce(ftmp3, tmp); /* 2^7 - 2^3 */ 744 felem_square(tmp, ftmp3); 745 felem_reduce(ftmp3, tmp); /* 2^8 - 2^4 */ 746 felem_assign(ftmp4, ftmp3); 747 felem_mul(tmp, ftmp3, ftmp); 748 felem_reduce(ftmp4, tmp); /* 2^8 - 2^1 */ 749 felem_square(tmp, ftmp4); 750 felem_reduce(ftmp4, tmp); /* 2^9 - 2^2 */ 751 felem_mul(tmp, ftmp3, ftmp2); 752 felem_reduce(ftmp3, tmp); /* 2^8 - 2^0 */ 753 felem_assign(ftmp2, ftmp3); 754 755 for (i = 0; i < 8; i++) { 756 felem_square(tmp, ftmp3); 757 felem_reduce(ftmp3, tmp); /* 2^16 - 2^8 */ 758 } 759 felem_mul(tmp, ftmp3, ftmp2); 760 felem_reduce(ftmp3, tmp); /* 2^16 - 2^0 */ 761 felem_assign(ftmp2, ftmp3); 762 763 for (i = 0; i < 16; i++) { 764 felem_square(tmp, ftmp3); 765 felem_reduce(ftmp3, tmp); /* 2^32 - 2^16 */ 766 } 767 felem_mul(tmp, ftmp3, ftmp2); 768 felem_reduce(ftmp3, tmp); /* 2^32 - 2^0 */ 769 felem_assign(ftmp2, ftmp3); 770 771 for (i = 0; i < 32; i++) { 772 felem_square(tmp, ftmp3); 773 felem_reduce(ftmp3, tmp); /* 2^64 - 2^32 */ 774 } 775 felem_mul(tmp, ftmp3, ftmp2); 776 felem_reduce(ftmp3, tmp); /* 2^64 - 2^0 */ 777 felem_assign(ftmp2, ftmp3); 778 779 for (i = 0; i < 64; i++) { 780 felem_square(tmp, ftmp3); 781 felem_reduce(ftmp3, tmp); /* 2^128 - 2^64 */ 782 } 783 felem_mul(tmp, ftmp3, ftmp2); 784 felem_reduce(ftmp3, tmp); /* 2^128 - 2^0 */ 785 felem_assign(ftmp2, ftmp3); 786 787 for (i = 0; i < 128; i++) { 788 felem_square(tmp, ftmp3); 789 felem_reduce(ftmp3, tmp); /* 2^256 - 2^128 */ 790 } 791 felem_mul(tmp, ftmp3, ftmp2); 792 felem_reduce(ftmp3, tmp); /* 2^256 - 2^0 */ 793 felem_assign(ftmp2, ftmp3); 794 795 for (i = 0; i < 256; i++) { 796 felem_square(tmp, ftmp3); 797 felem_reduce(ftmp3, tmp); /* 2^512 - 2^256 */ 798 } 799 felem_mul(tmp, ftmp3, ftmp2); 800 felem_reduce(ftmp3, tmp); /* 2^512 - 2^0 */ 801 802 for (i = 0; i < 9; i++) { 803 felem_square(tmp, ftmp3); 804 felem_reduce(ftmp3, tmp); /* 2^521 - 2^9 */ 805 } 806 felem_mul(tmp, ftmp3, ftmp4); 807 felem_reduce(ftmp3, tmp); /* 2^512 - 2^2 */ 808 felem_mul(tmp, ftmp3, in); 809 felem_reduce(out, tmp); /* 2^512 - 3 */ 810 } 811 812 /* This is 2^521-1, expressed as an felem */ 813 static const felem kPrime = { 814 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff, 815 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff, 816 0x03ffffffffffffff, 0x03ffffffffffffff, 0x01ffffffffffffff 817 }; 818 819 /*- 820 * felem_is_zero returns a limb with all bits set if |in| == 0 (mod p) and 0 821 * otherwise. 822 * On entry: 823 * in[i] < 2^59 + 2^14 824 */ 825 static limb felem_is_zero(const felem in) 826 { 827 felem ftmp; 828 limb is_zero, is_p; 829 felem_assign(ftmp, in); 830 831 ftmp[0] += ftmp[8] >> 57; 832 ftmp[8] &= bottom57bits; 833 /* ftmp[8] < 2^57 */ 834 ftmp[1] += ftmp[0] >> 58; 835 ftmp[0] &= bottom58bits; 836 ftmp[2] += ftmp[1] >> 58; 837 ftmp[1] &= bottom58bits; 838 ftmp[3] += ftmp[2] >> 58; 839 ftmp[2] &= bottom58bits; 840 ftmp[4] += ftmp[3] >> 58; 841 ftmp[3] &= bottom58bits; 842 ftmp[5] += ftmp[4] >> 58; 843 ftmp[4] &= bottom58bits; 844 ftmp[6] += ftmp[5] >> 58; 845 ftmp[5] &= bottom58bits; 846 ftmp[7] += ftmp[6] >> 58; 847 ftmp[6] &= bottom58bits; 848 ftmp[8] += ftmp[7] >> 58; 849 ftmp[7] &= bottom58bits; 850 /* ftmp[8] < 2^57 + 4 */ 851 852 /* 853 * The ninth limb of 2*(2^521-1) is 0x03ffffffffffffff, which is greater 854 * than our bound for ftmp[8]. Therefore we only have to check if the 855 * zero is zero or 2^521-1. 856 */ 857 858 is_zero = 0; 859 is_zero |= ftmp[0]; 860 is_zero |= ftmp[1]; 861 is_zero |= ftmp[2]; 862 is_zero |= ftmp[3]; 863 is_zero |= ftmp[4]; 864 is_zero |= ftmp[5]; 865 is_zero |= ftmp[6]; 866 is_zero |= ftmp[7]; 867 is_zero |= ftmp[8]; 868 869 is_zero--; 870 /* 871 * We know that ftmp[i] < 2^63, therefore the only way that the top bit 872 * can be set is if is_zero was 0 before the decrement. 873 */ 874 is_zero = 0 - (is_zero >> 63); 875 876 is_p = ftmp[0] ^ kPrime[0]; 877 is_p |= ftmp[1] ^ kPrime[1]; 878 is_p |= ftmp[2] ^ kPrime[2]; 879 is_p |= ftmp[3] ^ kPrime[3]; 880 is_p |= ftmp[4] ^ kPrime[4]; 881 is_p |= ftmp[5] ^ kPrime[5]; 882 is_p |= ftmp[6] ^ kPrime[6]; 883 is_p |= ftmp[7] ^ kPrime[7]; 884 is_p |= ftmp[8] ^ kPrime[8]; 885 886 is_p--; 887 is_p = 0 - (is_p >> 63); 888 889 is_zero |= is_p; 890 return is_zero; 891 } 892 893 static int felem_is_zero_int(const void *in) 894 { 895 return (int)(felem_is_zero(in) & ((limb) 1)); 896 } 897 898 /*- 899 * felem_contract converts |in| to its unique, minimal representation. 900 * On entry: 901 * in[i] < 2^59 + 2^14 902 */ 903 static void felem_contract(felem out, const felem in) 904 { 905 limb is_p, is_greater, sign; 906 static const limb two58 = ((limb) 1) << 58; 907 908 felem_assign(out, in); 909 910 out[0] += out[8] >> 57; 911 out[8] &= bottom57bits; 912 /* out[8] < 2^57 */ 913 out[1] += out[0] >> 58; 914 out[0] &= bottom58bits; 915 out[2] += out[1] >> 58; 916 out[1] &= bottom58bits; 917 out[3] += out[2] >> 58; 918 out[2] &= bottom58bits; 919 out[4] += out[3] >> 58; 920 out[3] &= bottom58bits; 921 out[5] += out[4] >> 58; 922 out[4] &= bottom58bits; 923 out[6] += out[5] >> 58; 924 out[5] &= bottom58bits; 925 out[7] += out[6] >> 58; 926 out[6] &= bottom58bits; 927 out[8] += out[7] >> 58; 928 out[7] &= bottom58bits; 929 /* out[8] < 2^57 + 4 */ 930 931 /* 932 * If the value is greater than 2^521-1 then we have to subtract 2^521-1 933 * out. See the comments in felem_is_zero regarding why we don't test for 934 * other multiples of the prime. 935 */ 936 937 /* 938 * First, if |out| is equal to 2^521-1, we subtract it out to get zero. 939 */ 940 941 is_p = out[0] ^ kPrime[0]; 942 is_p |= out[1] ^ kPrime[1]; 943 is_p |= out[2] ^ kPrime[2]; 944 is_p |= out[3] ^ kPrime[3]; 945 is_p |= out[4] ^ kPrime[4]; 946 is_p |= out[5] ^ kPrime[5]; 947 is_p |= out[6] ^ kPrime[6]; 948 is_p |= out[7] ^ kPrime[7]; 949 is_p |= out[8] ^ kPrime[8]; 950 951 is_p--; 952 is_p &= is_p << 32; 953 is_p &= is_p << 16; 954 is_p &= is_p << 8; 955 is_p &= is_p << 4; 956 is_p &= is_p << 2; 957 is_p &= is_p << 1; 958 is_p = 0 - (is_p >> 63); 959 is_p = ~is_p; 960 961 /* is_p is 0 iff |out| == 2^521-1 and all ones otherwise */ 962 963 out[0] &= is_p; 964 out[1] &= is_p; 965 out[2] &= is_p; 966 out[3] &= is_p; 967 out[4] &= is_p; 968 out[5] &= is_p; 969 out[6] &= is_p; 970 out[7] &= is_p; 971 out[8] &= is_p; 972 973 /* 974 * In order to test that |out| >= 2^521-1 we need only test if out[8] >> 975 * 57 is greater than zero as (2^521-1) + x >= 2^522 976 */ 977 is_greater = out[8] >> 57; 978 is_greater |= is_greater << 32; 979 is_greater |= is_greater << 16; 980 is_greater |= is_greater << 8; 981 is_greater |= is_greater << 4; 982 is_greater |= is_greater << 2; 983 is_greater |= is_greater << 1; 984 is_greater = 0 - (is_greater >> 63); 985 986 out[0] -= kPrime[0] & is_greater; 987 out[1] -= kPrime[1] & is_greater; 988 out[2] -= kPrime[2] & is_greater; 989 out[3] -= kPrime[3] & is_greater; 990 out[4] -= kPrime[4] & is_greater; 991 out[5] -= kPrime[5] & is_greater; 992 out[6] -= kPrime[6] & is_greater; 993 out[7] -= kPrime[7] & is_greater; 994 out[8] -= kPrime[8] & is_greater; 995 996 /* Eliminate negative coefficients */ 997 sign = -(out[0] >> 63); 998 out[0] += (two58 & sign); 999 out[1] -= (1 & sign); 1000 sign = -(out[1] >> 63); 1001 out[1] += (two58 & sign); 1002 out[2] -= (1 & sign); 1003 sign = -(out[2] >> 63); 1004 out[2] += (two58 & sign); 1005 out[3] -= (1 & sign); 1006 sign = -(out[3] >> 63); 1007 out[3] += (two58 & sign); 1008 out[4] -= (1 & sign); 1009 sign = -(out[4] >> 63); 1010 out[4] += (two58 & sign); 1011 out[5] -= (1 & sign); 1012 sign = -(out[0] >> 63); 1013 out[5] += (two58 & sign); 1014 out[6] -= (1 & sign); 1015 sign = -(out[6] >> 63); 1016 out[6] += (two58 & sign); 1017 out[7] -= (1 & sign); 1018 sign = -(out[7] >> 63); 1019 out[7] += (two58 & sign); 1020 out[8] -= (1 & sign); 1021 sign = -(out[5] >> 63); 1022 out[5] += (two58 & sign); 1023 out[6] -= (1 & sign); 1024 sign = -(out[6] >> 63); 1025 out[6] += (two58 & sign); 1026 out[7] -= (1 & sign); 1027 sign = -(out[7] >> 63); 1028 out[7] += (two58 & sign); 1029 out[8] -= (1 & sign); 1030 } 1031 1032 /*- 1033 * Group operations 1034 * ---------------- 1035 * 1036 * Building on top of the field operations we have the operations on the 1037 * elliptic curve group itself. Points on the curve are represented in Jacobian 1038 * coordinates */ 1039 1040 /*- 1041 * point_double calculates 2*(x_in, y_in, z_in) 1042 * 1043 * The method is taken from: 1044 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b 1045 * 1046 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed. 1047 * while x_out == y_in is not (maybe this works, but it's not tested). */ 1048 static void 1049 point_double(felem x_out, felem y_out, felem z_out, 1050 const felem x_in, const felem y_in, const felem z_in) 1051 { 1052 largefelem tmp, tmp2; 1053 felem delta, gamma, beta, alpha, ftmp, ftmp2; 1054 1055 felem_assign(ftmp, x_in); 1056 felem_assign(ftmp2, x_in); 1057 1058 /* delta = z^2 */ 1059 felem_square(tmp, z_in); 1060 felem_reduce(delta, tmp); /* delta[i] < 2^59 + 2^14 */ 1061 1062 /* gamma = y^2 */ 1063 felem_square(tmp, y_in); 1064 felem_reduce(gamma, tmp); /* gamma[i] < 2^59 + 2^14 */ 1065 1066 /* beta = x*gamma */ 1067 felem_mul(tmp, x_in, gamma); 1068 felem_reduce(beta, tmp); /* beta[i] < 2^59 + 2^14 */ 1069 1070 /* alpha = 3*(x-delta)*(x+delta) */ 1071 felem_diff64(ftmp, delta); 1072 /* ftmp[i] < 2^61 */ 1073 felem_sum64(ftmp2, delta); 1074 /* ftmp2[i] < 2^60 + 2^15 */ 1075 felem_scalar64(ftmp2, 3); 1076 /* ftmp2[i] < 3*2^60 + 3*2^15 */ 1077 felem_mul(tmp, ftmp, ftmp2); 1078 /*- 1079 * tmp[i] < 17(3*2^121 + 3*2^76) 1080 * = 61*2^121 + 61*2^76 1081 * < 64*2^121 + 64*2^76 1082 * = 2^127 + 2^82 1083 * < 2^128 1084 */ 1085 felem_reduce(alpha, tmp); 1086 1087 /* x' = alpha^2 - 8*beta */ 1088 felem_square(tmp, alpha); 1089 /* 1090 * tmp[i] < 17*2^120 < 2^125 1091 */ 1092 felem_assign(ftmp, beta); 1093 felem_scalar64(ftmp, 8); 1094 /* ftmp[i] < 2^62 + 2^17 */ 1095 felem_diff_128_64(tmp, ftmp); 1096 /* tmp[i] < 2^125 + 2^63 + 2^62 + 2^17 */ 1097 felem_reduce(x_out, tmp); 1098 1099 /* z' = (y + z)^2 - gamma - delta */ 1100 felem_sum64(delta, gamma); 1101 /* delta[i] < 2^60 + 2^15 */ 1102 felem_assign(ftmp, y_in); 1103 felem_sum64(ftmp, z_in); 1104 /* ftmp[i] < 2^60 + 2^15 */ 1105 felem_square(tmp, ftmp); 1106 /* 1107 * tmp[i] < 17(2^122) < 2^127 1108 */ 1109 felem_diff_128_64(tmp, delta); 1110 /* tmp[i] < 2^127 + 2^63 */ 1111 felem_reduce(z_out, tmp); 1112 1113 /* y' = alpha*(4*beta - x') - 8*gamma^2 */ 1114 felem_scalar64(beta, 4); 1115 /* beta[i] < 2^61 + 2^16 */ 1116 felem_diff64(beta, x_out); 1117 /* beta[i] < 2^61 + 2^60 + 2^16 */ 1118 felem_mul(tmp, alpha, beta); 1119 /*- 1120 * tmp[i] < 17*((2^59 + 2^14)(2^61 + 2^60 + 2^16)) 1121 * = 17*(2^120 + 2^75 + 2^119 + 2^74 + 2^75 + 2^30) 1122 * = 17*(2^120 + 2^119 + 2^76 + 2^74 + 2^30) 1123 * < 2^128 1124 */ 1125 felem_square(tmp2, gamma); 1126 /*- 1127 * tmp2[i] < 17*(2^59 + 2^14)^2 1128 * = 17*(2^118 + 2^74 + 2^28) 1129 */ 1130 felem_scalar128(tmp2, 8); 1131 /*- 1132 * tmp2[i] < 8*17*(2^118 + 2^74 + 2^28) 1133 * = 2^125 + 2^121 + 2^81 + 2^77 + 2^35 + 2^31 1134 * < 2^126 1135 */ 1136 felem_diff128(tmp, tmp2); 1137 /*- 1138 * tmp[i] < 2^127 - 2^69 + 17(2^120 + 2^119 + 2^76 + 2^74 + 2^30) 1139 * = 2^127 + 2^124 + 2^122 + 2^120 + 2^118 + 2^80 + 2^78 + 2^76 + 1140 * 2^74 + 2^69 + 2^34 + 2^30 1141 * < 2^128 1142 */ 1143 felem_reduce(y_out, tmp); 1144 } 1145 1146 /* copy_conditional copies in to out iff mask is all ones. */ 1147 static void copy_conditional(felem out, const felem in, limb mask) 1148 { 1149 unsigned i; 1150 for (i = 0; i < NLIMBS; ++i) { 1151 const limb tmp = mask & (in[i] ^ out[i]); 1152 out[i] ^= tmp; 1153 } 1154 } 1155 1156 /*- 1157 * point_add calculates (x1, y1, z1) + (x2, y2, z2) 1158 * 1159 * The method is taken from 1160 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl, 1161 * adapted for mixed addition (z2 = 1, or z2 = 0 for the point at infinity). 1162 * 1163 * This function includes a branch for checking whether the two input points 1164 * are equal (while not equal to the point at infinity). See comment below 1165 * on constant-time. 1166 */ 1167 static void point_add(felem x3, felem y3, felem z3, 1168 const felem x1, const felem y1, const felem z1, 1169 const int mixed, const felem x2, const felem y2, 1170 const felem z2) 1171 { 1172 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, ftmp6, x_out, y_out, z_out; 1173 largefelem tmp, tmp2; 1174 limb x_equal, y_equal, z1_is_zero, z2_is_zero; 1175 1176 z1_is_zero = felem_is_zero(z1); 1177 z2_is_zero = felem_is_zero(z2); 1178 1179 /* ftmp = z1z1 = z1**2 */ 1180 felem_square(tmp, z1); 1181 felem_reduce(ftmp, tmp); 1182 1183 if (!mixed) { 1184 /* ftmp2 = z2z2 = z2**2 */ 1185 felem_square(tmp, z2); 1186 felem_reduce(ftmp2, tmp); 1187 1188 /* u1 = ftmp3 = x1*z2z2 */ 1189 felem_mul(tmp, x1, ftmp2); 1190 felem_reduce(ftmp3, tmp); 1191 1192 /* ftmp5 = z1 + z2 */ 1193 felem_assign(ftmp5, z1); 1194 felem_sum64(ftmp5, z2); 1195 /* ftmp5[i] < 2^61 */ 1196 1197 /* ftmp5 = (z1 + z2)**2 - z1z1 - z2z2 = 2*z1z2 */ 1198 felem_square(tmp, ftmp5); 1199 /* tmp[i] < 17*2^122 */ 1200 felem_diff_128_64(tmp, ftmp); 1201 /* tmp[i] < 17*2^122 + 2^63 */ 1202 felem_diff_128_64(tmp, ftmp2); 1203 /* tmp[i] < 17*2^122 + 2^64 */ 1204 felem_reduce(ftmp5, tmp); 1205 1206 /* ftmp2 = z2 * z2z2 */ 1207 felem_mul(tmp, ftmp2, z2); 1208 felem_reduce(ftmp2, tmp); 1209 1210 /* s1 = ftmp6 = y1 * z2**3 */ 1211 felem_mul(tmp, y1, ftmp2); 1212 felem_reduce(ftmp6, tmp); 1213 } else { 1214 /* 1215 * We'll assume z2 = 1 (special case z2 = 0 is handled later) 1216 */ 1217 1218 /* u1 = ftmp3 = x1*z2z2 */ 1219 felem_assign(ftmp3, x1); 1220 1221 /* ftmp5 = 2*z1z2 */ 1222 felem_scalar(ftmp5, z1, 2); 1223 1224 /* s1 = ftmp6 = y1 * z2**3 */ 1225 felem_assign(ftmp6, y1); 1226 } 1227 1228 /* u2 = x2*z1z1 */ 1229 felem_mul(tmp, x2, ftmp); 1230 /* tmp[i] < 17*2^120 */ 1231 1232 /* h = ftmp4 = u2 - u1 */ 1233 felem_diff_128_64(tmp, ftmp3); 1234 /* tmp[i] < 17*2^120 + 2^63 */ 1235 felem_reduce(ftmp4, tmp); 1236 1237 x_equal = felem_is_zero(ftmp4); 1238 1239 /* z_out = ftmp5 * h */ 1240 felem_mul(tmp, ftmp5, ftmp4); 1241 felem_reduce(z_out, tmp); 1242 1243 /* ftmp = z1 * z1z1 */ 1244 felem_mul(tmp, ftmp, z1); 1245 felem_reduce(ftmp, tmp); 1246 1247 /* s2 = tmp = y2 * z1**3 */ 1248 felem_mul(tmp, y2, ftmp); 1249 /* tmp[i] < 17*2^120 */ 1250 1251 /* r = ftmp5 = (s2 - s1)*2 */ 1252 felem_diff_128_64(tmp, ftmp6); 1253 /* tmp[i] < 17*2^120 + 2^63 */ 1254 felem_reduce(ftmp5, tmp); 1255 y_equal = felem_is_zero(ftmp5); 1256 felem_scalar64(ftmp5, 2); 1257 /* ftmp5[i] < 2^61 */ 1258 1259 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) { 1260 /* 1261 * This is obviously not constant-time but it will almost-never happen 1262 * for ECDH / ECDSA. The case where it can happen is during scalar-mult 1263 * where the intermediate value gets very close to the group order. 1264 * Since |ec_GFp_nistp_recode_scalar_bits| produces signed digits for 1265 * the scalar, it's possible for the intermediate value to be a small 1266 * negative multiple of the base point, and for the final signed digit 1267 * to be the same value. We believe that this only occurs for the scalar 1268 * 1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffff 1269 * ffffffa51868783bf2f966b7fcc0148f709a5d03bb5c9b8899c47aebb6fb 1270 * 71e913863f7, in that case the penultimate intermediate is -9G and 1271 * the final digit is also -9G. Since this only happens for a single 1272 * scalar, the timing leak is irrelevent. (Any attacker who wanted to 1273 * check whether a secret scalar was that exact value, can already do 1274 * so.) 1275 */ 1276 point_double(x3, y3, z3, x1, y1, z1); 1277 return; 1278 } 1279 1280 /* I = ftmp = (2h)**2 */ 1281 felem_assign(ftmp, ftmp4); 1282 felem_scalar64(ftmp, 2); 1283 /* ftmp[i] < 2^61 */ 1284 felem_square(tmp, ftmp); 1285 /* tmp[i] < 17*2^122 */ 1286 felem_reduce(ftmp, tmp); 1287 1288 /* J = ftmp2 = h * I */ 1289 felem_mul(tmp, ftmp4, ftmp); 1290 felem_reduce(ftmp2, tmp); 1291 1292 /* V = ftmp4 = U1 * I */ 1293 felem_mul(tmp, ftmp3, ftmp); 1294 felem_reduce(ftmp4, tmp); 1295 1296 /* x_out = r**2 - J - 2V */ 1297 felem_square(tmp, ftmp5); 1298 /* tmp[i] < 17*2^122 */ 1299 felem_diff_128_64(tmp, ftmp2); 1300 /* tmp[i] < 17*2^122 + 2^63 */ 1301 felem_assign(ftmp3, ftmp4); 1302 felem_scalar64(ftmp4, 2); 1303 /* ftmp4[i] < 2^61 */ 1304 felem_diff_128_64(tmp, ftmp4); 1305 /* tmp[i] < 17*2^122 + 2^64 */ 1306 felem_reduce(x_out, tmp); 1307 1308 /* y_out = r(V-x_out) - 2 * s1 * J */ 1309 felem_diff64(ftmp3, x_out); 1310 /* 1311 * ftmp3[i] < 2^60 + 2^60 = 2^61 1312 */ 1313 felem_mul(tmp, ftmp5, ftmp3); 1314 /* tmp[i] < 17*2^122 */ 1315 felem_mul(tmp2, ftmp6, ftmp2); 1316 /* tmp2[i] < 17*2^120 */ 1317 felem_scalar128(tmp2, 2); 1318 /* tmp2[i] < 17*2^121 */ 1319 felem_diff128(tmp, tmp2); 1320 /*- 1321 * tmp[i] < 2^127 - 2^69 + 17*2^122 1322 * = 2^126 - 2^122 - 2^6 - 2^2 - 1 1323 * < 2^127 1324 */ 1325 felem_reduce(y_out, tmp); 1326 1327 copy_conditional(x_out, x2, z1_is_zero); 1328 copy_conditional(x_out, x1, z2_is_zero); 1329 copy_conditional(y_out, y2, z1_is_zero); 1330 copy_conditional(y_out, y1, z2_is_zero); 1331 copy_conditional(z_out, z2, z1_is_zero); 1332 copy_conditional(z_out, z1, z2_is_zero); 1333 felem_assign(x3, x_out); 1334 felem_assign(y3, y_out); 1335 felem_assign(z3, z_out); 1336 } 1337 1338 /*- 1339 * Base point pre computation 1340 * -------------------------- 1341 * 1342 * Two different sorts of precomputed tables are used in the following code. 1343 * Each contain various points on the curve, where each point is three field 1344 * elements (x, y, z). 1345 * 1346 * For the base point table, z is usually 1 (0 for the point at infinity). 1347 * This table has 16 elements: 1348 * index | bits | point 1349 * ------+---------+------------------------------ 1350 * 0 | 0 0 0 0 | 0G 1351 * 1 | 0 0 0 1 | 1G 1352 * 2 | 0 0 1 0 | 2^130G 1353 * 3 | 0 0 1 1 | (2^130 + 1)G 1354 * 4 | 0 1 0 0 | 2^260G 1355 * 5 | 0 1 0 1 | (2^260 + 1)G 1356 * 6 | 0 1 1 0 | (2^260 + 2^130)G 1357 * 7 | 0 1 1 1 | (2^260 + 2^130 + 1)G 1358 * 8 | 1 0 0 0 | 2^390G 1359 * 9 | 1 0 0 1 | (2^390 + 1)G 1360 * 10 | 1 0 1 0 | (2^390 + 2^130)G 1361 * 11 | 1 0 1 1 | (2^390 + 2^130 + 1)G 1362 * 12 | 1 1 0 0 | (2^390 + 2^260)G 1363 * 13 | 1 1 0 1 | (2^390 + 2^260 + 1)G 1364 * 14 | 1 1 1 0 | (2^390 + 2^260 + 2^130)G 1365 * 15 | 1 1 1 1 | (2^390 + 2^260 + 2^130 + 1)G 1366 * 1367 * The reason for this is so that we can clock bits into four different 1368 * locations when doing simple scalar multiplies against the base point. 1369 * 1370 * Tables for other points have table[i] = iG for i in 0 .. 16. */ 1371 1372 /* gmul is the table of precomputed base points */ 1373 static const felem gmul[16][3] = { 1374 {{0, 0, 0, 0, 0, 0, 0, 0, 0}, 1375 {0, 0, 0, 0, 0, 0, 0, 0, 0}, 1376 {0, 0, 0, 0, 0, 0, 0, 0, 0}}, 1377 {{0x017e7e31c2e5bd66, 0x022cf0615a90a6fe, 0x00127a2ffa8de334, 1378 0x01dfbf9d64a3f877, 0x006b4d3dbaa14b5e, 0x014fed487e0a2bd8, 1379 0x015b4429c6481390, 0x03a73678fb2d988e, 0x00c6858e06b70404}, 1380 {0x00be94769fd16650, 0x031c21a89cb09022, 0x039013fad0761353, 1381 0x02657bd099031542, 0x03273e662c97ee72, 0x01e6d11a05ebef45, 1382 0x03d1bd998f544495, 0x03001172297ed0b1, 0x011839296a789a3b}, 1383 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1384 {{0x0373faacbc875bae, 0x00f325023721c671, 0x00f666fd3dbde5ad, 1385 0x01a6932363f88ea7, 0x01fc6d9e13f9c47b, 0x03bcbffc2bbf734e, 1386 0x013ee3c3647f3a92, 0x029409fefe75d07d, 0x00ef9199963d85e5}, 1387 {0x011173743ad5b178, 0x02499c7c21bf7d46, 0x035beaeabb8b1a58, 1388 0x00f989c4752ea0a3, 0x0101e1de48a9c1a3, 0x01a20076be28ba6c, 1389 0x02f8052e5eb2de95, 0x01bfe8f82dea117c, 0x0160074d3c36ddb7}, 1390 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1391 {{0x012f3fc373393b3b, 0x03d3d6172f1419fa, 0x02adc943c0b86873, 1392 0x00d475584177952b, 0x012a4d1673750ee2, 0x00512517a0f13b0c, 1393 0x02b184671a7b1734, 0x0315b84236f1a50a, 0x00a4afc472edbdb9}, 1394 {0x00152a7077f385c4, 0x03044007d8d1c2ee, 0x0065829d61d52b52, 1395 0x00494ff6b6631d0d, 0x00a11d94d5f06bcf, 0x02d2f89474d9282e, 1396 0x0241c5727c06eeb9, 0x0386928710fbdb9d, 0x01f883f727b0dfbe}, 1397 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1398 {{0x019b0c3c9185544d, 0x006243a37c9d97db, 0x02ee3cbe030a2ad2, 1399 0x00cfdd946bb51e0d, 0x0271c00932606b91, 0x03f817d1ec68c561, 1400 0x03f37009806a369c, 0x03c1f30baf184fd5, 0x01091022d6d2f065}, 1401 {0x0292c583514c45ed, 0x0316fca51f9a286c, 0x00300af507c1489a, 1402 0x0295f69008298cf1, 0x02c0ed8274943d7b, 0x016509b9b47a431e, 1403 0x02bc9de9634868ce, 0x005b34929bffcb09, 0x000c1a0121681524}, 1404 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1405 {{0x0286abc0292fb9f2, 0x02665eee9805b3f7, 0x01ed7455f17f26d6, 1406 0x0346355b83175d13, 0x006284944cd0a097, 0x0191895bcdec5e51, 1407 0x02e288370afda7d9, 0x03b22312bfefa67a, 0x01d104d3fc0613fe}, 1408 {0x0092421a12f7e47f, 0x0077a83fa373c501, 0x03bd25c5f696bd0d, 1409 0x035c41e4d5459761, 0x01ca0d1742b24f53, 0x00aaab27863a509c, 1410 0x018b6de47df73917, 0x025c0b771705cd01, 0x01fd51d566d760a7}, 1411 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1412 {{0x01dd92ff6b0d1dbd, 0x039c5e2e8f8afa69, 0x0261ed13242c3b27, 1413 0x0382c6e67026e6a0, 0x01d60b10be2089f9, 0x03c15f3dce86723f, 1414 0x03c764a32d2a062d, 0x017307eac0fad056, 0x018207c0b96c5256}, 1415 {0x0196a16d60e13154, 0x03e6ce74c0267030, 0x00ddbf2b4e52a5aa, 1416 0x012738241bbf31c8, 0x00ebe8dc04685a28, 0x024c2ad6d380d4a2, 1417 0x035ee062a6e62d0e, 0x0029ed74af7d3a0f, 0x00eef32aec142ebd}, 1418 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1419 {{0x00c31ec398993b39, 0x03a9f45bcda68253, 0x00ac733c24c70890, 1420 0x00872b111401ff01, 0x01d178c23195eafb, 0x03bca2c816b87f74, 1421 0x0261a9af46fbad7a, 0x0324b2a8dd3d28f9, 0x00918121d8f24e23}, 1422 {0x032bc8c1ca983cd7, 0x00d869dfb08fc8c6, 0x01693cb61fce1516, 1423 0x012a5ea68f4e88a8, 0x010869cab88d7ae3, 0x009081ad277ceee1, 1424 0x033a77166d064cdc, 0x03955235a1fb3a95, 0x01251a4a9b25b65e}, 1425 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1426 {{0x00148a3a1b27f40b, 0x0123186df1b31fdc, 0x00026e7beaad34ce, 1427 0x01db446ac1d3dbba, 0x0299c1a33437eaec, 0x024540610183cbb7, 1428 0x0173bb0e9ce92e46, 0x02b937e43921214b, 0x01ab0436a9bf01b5}, 1429 {0x0383381640d46948, 0x008dacbf0e7f330f, 0x03602122bcc3f318, 1430 0x01ee596b200620d6, 0x03bd0585fda430b3, 0x014aed77fd123a83, 1431 0x005ace749e52f742, 0x0390fe041da2b842, 0x0189a8ceb3299242}, 1432 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1433 {{0x012a19d6b3282473, 0x00c0915918b423ce, 0x023a954eb94405ae, 1434 0x00529f692be26158, 0x0289fa1b6fa4b2aa, 0x0198ae4ceea346ef, 1435 0x0047d8cdfbdedd49, 0x00cc8c8953f0f6b8, 0x001424abbff49203}, 1436 {0x0256732a1115a03a, 0x0351bc38665c6733, 0x03f7b950fb4a6447, 1437 0x000afffa94c22155, 0x025763d0a4dab540, 0x000511e92d4fc283, 1438 0x030a7e9eda0ee96c, 0x004c3cd93a28bf0a, 0x017edb3a8719217f}, 1439 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1440 {{0x011de5675a88e673, 0x031d7d0f5e567fbe, 0x0016b2062c970ae5, 1441 0x03f4a2be49d90aa7, 0x03cef0bd13822866, 0x03f0923dcf774a6c, 1442 0x0284bebc4f322f72, 0x016ab2645302bb2c, 0x01793f95dace0e2a}, 1443 {0x010646e13527a28f, 0x01ca1babd59dc5e7, 0x01afedfd9a5595df, 1444 0x01f15785212ea6b1, 0x0324e5d64f6ae3f4, 0x02d680f526d00645, 1445 0x0127920fadf627a7, 0x03b383f75df4f684, 0x0089e0057e783b0a}, 1446 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1447 {{0x00f334b9eb3c26c6, 0x0298fdaa98568dce, 0x01c2d24843a82292, 1448 0x020bcb24fa1b0711, 0x02cbdb3d2b1875e6, 0x0014907598f89422, 1449 0x03abe3aa43b26664, 0x02cbf47f720bc168, 0x0133b5e73014b79b}, 1450 {0x034aab5dab05779d, 0x00cdc5d71fee9abb, 0x0399f16bd4bd9d30, 1451 0x03582fa592d82647, 0x02be1cdfb775b0e9, 0x0034f7cea32e94cb, 1452 0x0335a7f08f56f286, 0x03b707e9565d1c8b, 0x0015c946ea5b614f}, 1453 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1454 {{0x024676f6cff72255, 0x00d14625cac96378, 0x00532b6008bc3767, 1455 0x01fc16721b985322, 0x023355ea1b091668, 0x029de7afdc0317c3, 1456 0x02fc8a7ca2da037c, 0x02de1217d74a6f30, 0x013f7173175b73bf}, 1457 {0x0344913f441490b5, 0x0200f9e272b61eca, 0x0258a246b1dd55d2, 1458 0x03753db9ea496f36, 0x025e02937a09c5ef, 0x030cbd3d14012692, 1459 0x01793a67e70dc72a, 0x03ec1d37048a662e, 0x006550f700c32a8d}, 1460 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1461 {{0x00d3f48a347eba27, 0x008e636649b61bd8, 0x00d3b93716778fb3, 1462 0x004d1915757bd209, 0x019d5311a3da44e0, 0x016d1afcbbe6aade, 1463 0x0241bf5f73265616, 0x0384672e5d50d39b, 0x005009fee522b684}, 1464 {0x029b4fab064435fe, 0x018868ee095bbb07, 0x01ea3d6936cc92b8, 1465 0x000608b00f78a2f3, 0x02db911073d1c20f, 0x018205938470100a, 1466 0x01f1e4964cbe6ff2, 0x021a19a29eed4663, 0x01414485f42afa81}, 1467 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1468 {{0x01612b3a17f63e34, 0x03813992885428e6, 0x022b3c215b5a9608, 1469 0x029b4057e19f2fcb, 0x0384059a587af7e6, 0x02d6400ace6fe610, 1470 0x029354d896e8e331, 0x00c047ee6dfba65e, 0x0037720542e9d49d}, 1471 {0x02ce9eed7c5e9278, 0x0374ed703e79643b, 0x01316c54c4072006, 1472 0x005aaa09054b2ee8, 0x002824000c840d57, 0x03d4eba24771ed86, 1473 0x0189c50aabc3bdae, 0x0338c01541e15510, 0x00466d56e38eed42}, 1474 {1, 0, 0, 0, 0, 0, 0, 0, 0}}, 1475 {{0x007efd8330ad8bd6, 0x02465ed48047710b, 0x0034c6606b215e0c, 1476 0x016ae30c53cbf839, 0x01fa17bd37161216, 0x018ead4e61ce8ab9, 1477 0x005482ed5f5dee46, 0x037543755bba1d7f, 0x005e5ac7e70a9d0f}, 1478 {0x0117e1bb2fdcb2a2, 0x03deea36249f40c4, 0x028d09b4a6246cb7, 1479 0x03524b8855bcf756, 0x023d7d109d5ceb58, 0x0178e43e3223ef9c, 1480 0x0154536a0c6e966a, 0x037964d1286ee9fe, 0x0199bcd90e125055}, 1481 {1, 0, 0, 0, 0, 0, 0, 0, 0}} 1482 }; 1483 1484 /* 1485 * select_point selects the |idx|th point from a precomputation table and 1486 * copies it to out. 1487 */ 1488 /* pre_comp below is of the size provided in |size| */ 1489 static void select_point(const limb idx, unsigned int size, 1490 const felem pre_comp[][3], felem out[3]) 1491 { 1492 unsigned i, j; 1493 limb *outlimbs = &out[0][0]; 1494 1495 memset(out, 0, sizeof(*out) * 3); 1496 1497 for (i = 0; i < size; i++) { 1498 const limb *inlimbs = &pre_comp[i][0][0]; 1499 limb mask = i ^ idx; 1500 mask |= mask >> 4; 1501 mask |= mask >> 2; 1502 mask |= mask >> 1; 1503 mask &= 1; 1504 mask--; 1505 for (j = 0; j < NLIMBS * 3; j++) 1506 outlimbs[j] |= inlimbs[j] & mask; 1507 } 1508 } 1509 1510 /* get_bit returns the |i|th bit in |in| */ 1511 static char get_bit(const felem_bytearray in, int i) 1512 { 1513 if (i < 0) 1514 return 0; 1515 return (in[i >> 3] >> (i & 7)) & 1; 1516 } 1517 1518 /* 1519 * Interleaved point multiplication using precomputed point multiples: The 1520 * small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[], the scalars 1521 * in scalars[]. If g_scalar is non-NULL, we also add this multiple of the 1522 * generator, using certain (large) precomputed multiples in g_pre_comp. 1523 * Output point (X, Y, Z) is stored in x_out, y_out, z_out 1524 */ 1525 static void batch_mul(felem x_out, felem y_out, felem z_out, 1526 const felem_bytearray scalars[], 1527 const unsigned num_points, const u8 *g_scalar, 1528 const int mixed, const felem pre_comp[][17][3], 1529 const felem g_pre_comp[16][3]) 1530 { 1531 int i, skip; 1532 unsigned num, gen_mul = (g_scalar != NULL); 1533 felem nq[3], tmp[4]; 1534 limb bits; 1535 u8 sign, digit; 1536 1537 /* set nq to the point at infinity */ 1538 memset(nq, 0, sizeof(nq)); 1539 1540 /* 1541 * Loop over all scalars msb-to-lsb, interleaving additions of multiples 1542 * of the generator (last quarter of rounds) and additions of other 1543 * points multiples (every 5th round). 1544 */ 1545 skip = 1; /* save two point operations in the first 1546 * round */ 1547 for (i = (num_points ? 520 : 130); i >= 0; --i) { 1548 /* double */ 1549 if (!skip) 1550 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]); 1551 1552 /* add multiples of the generator */ 1553 if (gen_mul && (i <= 130)) { 1554 bits = get_bit(g_scalar, i + 390) << 3; 1555 if (i < 130) { 1556 bits |= get_bit(g_scalar, i + 260) << 2; 1557 bits |= get_bit(g_scalar, i + 130) << 1; 1558 bits |= get_bit(g_scalar, i); 1559 } 1560 /* select the point to add, in constant time */ 1561 select_point(bits, 16, g_pre_comp, tmp); 1562 if (!skip) { 1563 /* The 1 argument below is for "mixed" */ 1564 point_add(nq[0], nq[1], nq[2], 1565 nq[0], nq[1], nq[2], 1, tmp[0], tmp[1], tmp[2]); 1566 } else { 1567 memcpy(nq, tmp, 3 * sizeof(felem)); 1568 skip = 0; 1569 } 1570 } 1571 1572 /* do other additions every 5 doublings */ 1573 if (num_points && (i % 5 == 0)) { 1574 /* loop over all scalars */ 1575 for (num = 0; num < num_points; ++num) { 1576 bits = get_bit(scalars[num], i + 4) << 5; 1577 bits |= get_bit(scalars[num], i + 3) << 4; 1578 bits |= get_bit(scalars[num], i + 2) << 3; 1579 bits |= get_bit(scalars[num], i + 1) << 2; 1580 bits |= get_bit(scalars[num], i) << 1; 1581 bits |= get_bit(scalars[num], i - 1); 1582 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits); 1583 1584 /* 1585 * select the point to add or subtract, in constant time 1586 */ 1587 select_point(digit, 17, pre_comp[num], tmp); 1588 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative 1589 * point */ 1590 copy_conditional(tmp[1], tmp[3], (-(limb) sign)); 1591 1592 if (!skip) { 1593 point_add(nq[0], nq[1], nq[2], 1594 nq[0], nq[1], nq[2], 1595 mixed, tmp[0], tmp[1], tmp[2]); 1596 } else { 1597 memcpy(nq, tmp, 3 * sizeof(felem)); 1598 skip = 0; 1599 } 1600 } 1601 } 1602 } 1603 felem_assign(x_out, nq[0]); 1604 felem_assign(y_out, nq[1]); 1605 felem_assign(z_out, nq[2]); 1606 } 1607 1608 /* Precomputation for the group generator. */ 1609 struct nistp521_pre_comp_st { 1610 felem g_pre_comp[16][3]; 1611 CRYPTO_REF_COUNT references; 1612 CRYPTO_RWLOCK *lock; 1613 }; 1614 1615 const EC_METHOD *EC_GFp_nistp521_method(void) 1616 { 1617 static const EC_METHOD ret = { 1618 EC_FLAGS_DEFAULT_OCT, 1619 NID_X9_62_prime_field, 1620 ec_GFp_nistp521_group_init, 1621 ec_GFp_simple_group_finish, 1622 ec_GFp_simple_group_clear_finish, 1623 ec_GFp_nist_group_copy, 1624 ec_GFp_nistp521_group_set_curve, 1625 ec_GFp_simple_group_get_curve, 1626 ec_GFp_simple_group_get_degree, 1627 ec_group_simple_order_bits, 1628 ec_GFp_simple_group_check_discriminant, 1629 ec_GFp_simple_point_init, 1630 ec_GFp_simple_point_finish, 1631 ec_GFp_simple_point_clear_finish, 1632 ec_GFp_simple_point_copy, 1633 ec_GFp_simple_point_set_to_infinity, 1634 ec_GFp_simple_set_Jprojective_coordinates_GFp, 1635 ec_GFp_simple_get_Jprojective_coordinates_GFp, 1636 ec_GFp_simple_point_set_affine_coordinates, 1637 ec_GFp_nistp521_point_get_affine_coordinates, 1638 0 /* point_set_compressed_coordinates */ , 1639 0 /* point2oct */ , 1640 0 /* oct2point */ , 1641 ec_GFp_simple_add, 1642 ec_GFp_simple_dbl, 1643 ec_GFp_simple_invert, 1644 ec_GFp_simple_is_at_infinity, 1645 ec_GFp_simple_is_on_curve, 1646 ec_GFp_simple_cmp, 1647 ec_GFp_simple_make_affine, 1648 ec_GFp_simple_points_make_affine, 1649 ec_GFp_nistp521_points_mul, 1650 ec_GFp_nistp521_precompute_mult, 1651 ec_GFp_nistp521_have_precompute_mult, 1652 ec_GFp_nist_field_mul, 1653 ec_GFp_nist_field_sqr, 1654 0 /* field_div */ , 1655 ec_GFp_simple_field_inv, 1656 0 /* field_encode */ , 1657 0 /* field_decode */ , 1658 0, /* field_set_to_one */ 1659 ec_key_simple_priv2oct, 1660 ec_key_simple_oct2priv, 1661 0, /* set private */ 1662 ec_key_simple_generate_key, 1663 ec_key_simple_check_key, 1664 ec_key_simple_generate_public_key, 1665 0, /* keycopy */ 1666 0, /* keyfinish */ 1667 ecdh_simple_compute_key, 1668 0, /* field_inverse_mod_ord */ 1669 0, /* blind_coordinates */ 1670 0, /* ladder_pre */ 1671 0, /* ladder_step */ 1672 0 /* ladder_post */ 1673 }; 1674 1675 return &ret; 1676 } 1677 1678 /******************************************************************************/ 1679 /* 1680 * FUNCTIONS TO MANAGE PRECOMPUTATION 1681 */ 1682 1683 static NISTP521_PRE_COMP *nistp521_pre_comp_new(void) 1684 { 1685 NISTP521_PRE_COMP *ret = OPENSSL_zalloc(sizeof(*ret)); 1686 1687 if (ret == NULL) { 1688 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE); 1689 return ret; 1690 } 1691 1692 ret->references = 1; 1693 1694 ret->lock = CRYPTO_THREAD_lock_new(); 1695 if (ret->lock == NULL) { 1696 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE); 1697 OPENSSL_free(ret); 1698 return NULL; 1699 } 1700 return ret; 1701 } 1702 1703 NISTP521_PRE_COMP *EC_nistp521_pre_comp_dup(NISTP521_PRE_COMP *p) 1704 { 1705 int i; 1706 if (p != NULL) 1707 CRYPTO_UP_REF(&p->references, &i, p->lock); 1708 return p; 1709 } 1710 1711 void EC_nistp521_pre_comp_free(NISTP521_PRE_COMP *p) 1712 { 1713 int i; 1714 1715 if (p == NULL) 1716 return; 1717 1718 CRYPTO_DOWN_REF(&p->references, &i, p->lock); 1719 REF_PRINT_COUNT("EC_nistp521", x); 1720 if (i > 0) 1721 return; 1722 REF_ASSERT_ISNT(i < 0); 1723 1724 CRYPTO_THREAD_lock_free(p->lock); 1725 OPENSSL_free(p); 1726 } 1727 1728 /******************************************************************************/ 1729 /* 1730 * OPENSSL EC_METHOD FUNCTIONS 1731 */ 1732 1733 int ec_GFp_nistp521_group_init(EC_GROUP *group) 1734 { 1735 int ret; 1736 ret = ec_GFp_simple_group_init(group); 1737 group->a_is_minus3 = 1; 1738 return ret; 1739 } 1740 1741 int ec_GFp_nistp521_group_set_curve(EC_GROUP *group, const BIGNUM *p, 1742 const BIGNUM *a, const BIGNUM *b, 1743 BN_CTX *ctx) 1744 { 1745 int ret = 0; 1746 BN_CTX *new_ctx = NULL; 1747 BIGNUM *curve_p, *curve_a, *curve_b; 1748 1749 if (ctx == NULL) 1750 if ((ctx = new_ctx = BN_CTX_new()) == NULL) 1751 return 0; 1752 BN_CTX_start(ctx); 1753 curve_p = BN_CTX_get(ctx); 1754 curve_a = BN_CTX_get(ctx); 1755 curve_b = BN_CTX_get(ctx); 1756 if (curve_b == NULL) 1757 goto err; 1758 BN_bin2bn(nistp521_curve_params[0], sizeof(felem_bytearray), curve_p); 1759 BN_bin2bn(nistp521_curve_params[1], sizeof(felem_bytearray), curve_a); 1760 BN_bin2bn(nistp521_curve_params[2], sizeof(felem_bytearray), curve_b); 1761 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) || (BN_cmp(curve_b, b))) { 1762 ECerr(EC_F_EC_GFP_NISTP521_GROUP_SET_CURVE, 1763 EC_R_WRONG_CURVE_PARAMETERS); 1764 goto err; 1765 } 1766 group->field_mod_func = BN_nist_mod_521; 1767 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx); 1768 err: 1769 BN_CTX_end(ctx); 1770 BN_CTX_free(new_ctx); 1771 return ret; 1772 } 1773 1774 /* 1775 * Takes the Jacobian coordinates (X, Y, Z) of a point and returns (X', Y') = 1776 * (X/Z^2, Y/Z^3) 1777 */ 1778 int ec_GFp_nistp521_point_get_affine_coordinates(const EC_GROUP *group, 1779 const EC_POINT *point, 1780 BIGNUM *x, BIGNUM *y, 1781 BN_CTX *ctx) 1782 { 1783 felem z1, z2, x_in, y_in, x_out, y_out; 1784 largefelem tmp; 1785 1786 if (EC_POINT_is_at_infinity(group, point)) { 1787 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES, 1788 EC_R_POINT_AT_INFINITY); 1789 return 0; 1790 } 1791 if ((!BN_to_felem(x_in, point->X)) || (!BN_to_felem(y_in, point->Y)) || 1792 (!BN_to_felem(z1, point->Z))) 1793 return 0; 1794 felem_inv(z2, z1); 1795 felem_square(tmp, z2); 1796 felem_reduce(z1, tmp); 1797 felem_mul(tmp, x_in, z1); 1798 felem_reduce(x_in, tmp); 1799 felem_contract(x_out, x_in); 1800 if (x != NULL) { 1801 if (!felem_to_BN(x, x_out)) { 1802 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES, 1803 ERR_R_BN_LIB); 1804 return 0; 1805 } 1806 } 1807 felem_mul(tmp, z1, z2); 1808 felem_reduce(z1, tmp); 1809 felem_mul(tmp, y_in, z1); 1810 felem_reduce(y_in, tmp); 1811 felem_contract(y_out, y_in); 1812 if (y != NULL) { 1813 if (!felem_to_BN(y, y_out)) { 1814 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES, 1815 ERR_R_BN_LIB); 1816 return 0; 1817 } 1818 } 1819 return 1; 1820 } 1821 1822 /* points below is of size |num|, and tmp_felems is of size |num+1/ */ 1823 static void make_points_affine(size_t num, felem points[][3], 1824 felem tmp_felems[]) 1825 { 1826 /* 1827 * Runs in constant time, unless an input is the point at infinity (which 1828 * normally shouldn't happen). 1829 */ 1830 ec_GFp_nistp_points_make_affine_internal(num, 1831 points, 1832 sizeof(felem), 1833 tmp_felems, 1834 (void (*)(void *))felem_one, 1835 felem_is_zero_int, 1836 (void (*)(void *, const void *)) 1837 felem_assign, 1838 (void (*)(void *, const void *)) 1839 felem_square_reduce, (void (*) 1840 (void *, 1841 const void 1842 *, 1843 const void 1844 *)) 1845 felem_mul_reduce, 1846 (void (*)(void *, const void *)) 1847 felem_inv, 1848 (void (*)(void *, const void *)) 1849 felem_contract); 1850 } 1851 1852 /* 1853 * Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL 1854 * values Result is stored in r (r can equal one of the inputs). 1855 */ 1856 int ec_GFp_nistp521_points_mul(const EC_GROUP *group, EC_POINT *r, 1857 const BIGNUM *scalar, size_t num, 1858 const EC_POINT *points[], 1859 const BIGNUM *scalars[], BN_CTX *ctx) 1860 { 1861 int ret = 0; 1862 int j; 1863 int mixed = 0; 1864 BIGNUM *x, *y, *z, *tmp_scalar; 1865 felem_bytearray g_secret; 1866 felem_bytearray *secrets = NULL; 1867 felem (*pre_comp)[17][3] = NULL; 1868 felem *tmp_felems = NULL; 1869 felem_bytearray tmp; 1870 unsigned i, num_bytes; 1871 int have_pre_comp = 0; 1872 size_t num_points = num; 1873 felem x_in, y_in, z_in, x_out, y_out, z_out; 1874 NISTP521_PRE_COMP *pre = NULL; 1875 felem(*g_pre_comp)[3] = NULL; 1876 EC_POINT *generator = NULL; 1877 const EC_POINT *p = NULL; 1878 const BIGNUM *p_scalar = NULL; 1879 1880 BN_CTX_start(ctx); 1881 x = BN_CTX_get(ctx); 1882 y = BN_CTX_get(ctx); 1883 z = BN_CTX_get(ctx); 1884 tmp_scalar = BN_CTX_get(ctx); 1885 if (tmp_scalar == NULL) 1886 goto err; 1887 1888 if (scalar != NULL) { 1889 pre = group->pre_comp.nistp521; 1890 if (pre) 1891 /* we have precomputation, try to use it */ 1892 g_pre_comp = &pre->g_pre_comp[0]; 1893 else 1894 /* try to use the standard precomputation */ 1895 g_pre_comp = (felem(*)[3]) gmul; 1896 generator = EC_POINT_new(group); 1897 if (generator == NULL) 1898 goto err; 1899 /* get the generator from precomputation */ 1900 if (!felem_to_BN(x, g_pre_comp[1][0]) || 1901 !felem_to_BN(y, g_pre_comp[1][1]) || 1902 !felem_to_BN(z, g_pre_comp[1][2])) { 1903 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB); 1904 goto err; 1905 } 1906 if (!EC_POINT_set_Jprojective_coordinates_GFp(group, 1907 generator, x, y, z, 1908 ctx)) 1909 goto err; 1910 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) 1911 /* precomputation matches generator */ 1912 have_pre_comp = 1; 1913 else 1914 /* 1915 * we don't have valid precomputation: treat the generator as a 1916 * random point 1917 */ 1918 num_points++; 1919 } 1920 1921 if (num_points > 0) { 1922 if (num_points >= 2) { 1923 /* 1924 * unless we precompute multiples for just one point, converting 1925 * those into affine form is time well spent 1926 */ 1927 mixed = 1; 1928 } 1929 secrets = OPENSSL_zalloc(sizeof(*secrets) * num_points); 1930 pre_comp = OPENSSL_zalloc(sizeof(*pre_comp) * num_points); 1931 if (mixed) 1932 tmp_felems = 1933 OPENSSL_malloc(sizeof(*tmp_felems) * (num_points * 17 + 1)); 1934 if ((secrets == NULL) || (pre_comp == NULL) 1935 || (mixed && (tmp_felems == NULL))) { 1936 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_MALLOC_FAILURE); 1937 goto err; 1938 } 1939 1940 /* 1941 * we treat NULL scalars as 0, and NULL points as points at infinity, 1942 * i.e., they contribute nothing to the linear combination 1943 */ 1944 for (i = 0; i < num_points; ++i) { 1945 if (i == num) 1946 /* 1947 * we didn't have a valid precomputation, so we pick the 1948 * generator 1949 */ 1950 { 1951 p = EC_GROUP_get0_generator(group); 1952 p_scalar = scalar; 1953 } else 1954 /* the i^th point */ 1955 { 1956 p = points[i]; 1957 p_scalar = scalars[i]; 1958 } 1959 if ((p_scalar != NULL) && (p != NULL)) { 1960 /* reduce scalar to 0 <= scalar < 2^521 */ 1961 if ((BN_num_bits(p_scalar) > 521) 1962 || (BN_is_negative(p_scalar))) { 1963 /* 1964 * this is an unusual input, and we don't guarantee 1965 * constant-timeness 1966 */ 1967 if (!BN_nnmod(tmp_scalar, p_scalar, group->order, ctx)) { 1968 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB); 1969 goto err; 1970 } 1971 num_bytes = BN_bn2bin(tmp_scalar, tmp); 1972 } else 1973 num_bytes = BN_bn2bin(p_scalar, tmp); 1974 flip_endian(secrets[i], tmp, num_bytes); 1975 /* precompute multiples */ 1976 if ((!BN_to_felem(x_out, p->X)) || 1977 (!BN_to_felem(y_out, p->Y)) || 1978 (!BN_to_felem(z_out, p->Z))) 1979 goto err; 1980 memcpy(pre_comp[i][1][0], x_out, sizeof(felem)); 1981 memcpy(pre_comp[i][1][1], y_out, sizeof(felem)); 1982 memcpy(pre_comp[i][1][2], z_out, sizeof(felem)); 1983 for (j = 2; j <= 16; ++j) { 1984 if (j & 1) { 1985 point_add(pre_comp[i][j][0], pre_comp[i][j][1], 1986 pre_comp[i][j][2], pre_comp[i][1][0], 1987 pre_comp[i][1][1], pre_comp[i][1][2], 0, 1988 pre_comp[i][j - 1][0], 1989 pre_comp[i][j - 1][1], 1990 pre_comp[i][j - 1][2]); 1991 } else { 1992 point_double(pre_comp[i][j][0], pre_comp[i][j][1], 1993 pre_comp[i][j][2], pre_comp[i][j / 2][0], 1994 pre_comp[i][j / 2][1], 1995 pre_comp[i][j / 2][2]); 1996 } 1997 } 1998 } 1999 } 2000 if (mixed) 2001 make_points_affine(num_points * 17, pre_comp[0], tmp_felems); 2002 } 2003 2004 /* the scalar for the generator */ 2005 if ((scalar != NULL) && (have_pre_comp)) { 2006 memset(g_secret, 0, sizeof(g_secret)); 2007 /* reduce scalar to 0 <= scalar < 2^521 */ 2008 if ((BN_num_bits(scalar) > 521) || (BN_is_negative(scalar))) { 2009 /* 2010 * this is an unusual input, and we don't guarantee 2011 * constant-timeness 2012 */ 2013 if (!BN_nnmod(tmp_scalar, scalar, group->order, ctx)) { 2014 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB); 2015 goto err; 2016 } 2017 num_bytes = BN_bn2bin(tmp_scalar, tmp); 2018 } else 2019 num_bytes = BN_bn2bin(scalar, tmp); 2020 flip_endian(g_secret, tmp, num_bytes); 2021 /* do the multiplication with generator precomputation */ 2022 batch_mul(x_out, y_out, z_out, 2023 (const felem_bytearray(*))secrets, num_points, 2024 g_secret, 2025 mixed, (const felem(*)[17][3])pre_comp, 2026 (const felem(*)[3])g_pre_comp); 2027 } else 2028 /* do the multiplication without generator precomputation */ 2029 batch_mul(x_out, y_out, z_out, 2030 (const felem_bytearray(*))secrets, num_points, 2031 NULL, mixed, (const felem(*)[17][3])pre_comp, NULL); 2032 /* reduce the output to its unique minimal representation */ 2033 felem_contract(x_in, x_out); 2034 felem_contract(y_in, y_out); 2035 felem_contract(z_in, z_out); 2036 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) || 2037 (!felem_to_BN(z, z_in))) { 2038 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB); 2039 goto err; 2040 } 2041 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx); 2042 2043 err: 2044 BN_CTX_end(ctx); 2045 EC_POINT_free(generator); 2046 OPENSSL_free(secrets); 2047 OPENSSL_free(pre_comp); 2048 OPENSSL_free(tmp_felems); 2049 return ret; 2050 } 2051 2052 int ec_GFp_nistp521_precompute_mult(EC_GROUP *group, BN_CTX *ctx) 2053 { 2054 int ret = 0; 2055 NISTP521_PRE_COMP *pre = NULL; 2056 int i, j; 2057 BN_CTX *new_ctx = NULL; 2058 BIGNUM *x, *y; 2059 EC_POINT *generator = NULL; 2060 felem tmp_felems[16]; 2061 2062 /* throw away old precomputation */ 2063 EC_pre_comp_free(group); 2064 if (ctx == NULL) 2065 if ((ctx = new_ctx = BN_CTX_new()) == NULL) 2066 return 0; 2067 BN_CTX_start(ctx); 2068 x = BN_CTX_get(ctx); 2069 y = BN_CTX_get(ctx); 2070 if (y == NULL) 2071 goto err; 2072 /* get the generator */ 2073 if (group->generator == NULL) 2074 goto err; 2075 generator = EC_POINT_new(group); 2076 if (generator == NULL) 2077 goto err; 2078 BN_bin2bn(nistp521_curve_params[3], sizeof(felem_bytearray), x); 2079 BN_bin2bn(nistp521_curve_params[4], sizeof(felem_bytearray), y); 2080 if (!EC_POINT_set_affine_coordinates(group, generator, x, y, ctx)) 2081 goto err; 2082 if ((pre = nistp521_pre_comp_new()) == NULL) 2083 goto err; 2084 /* 2085 * if the generator is the standard one, use built-in precomputation 2086 */ 2087 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) { 2088 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp)); 2089 goto done; 2090 } 2091 if ((!BN_to_felem(pre->g_pre_comp[1][0], group->generator->X)) || 2092 (!BN_to_felem(pre->g_pre_comp[1][1], group->generator->Y)) || 2093 (!BN_to_felem(pre->g_pre_comp[1][2], group->generator->Z))) 2094 goto err; 2095 /* compute 2^130*G, 2^260*G, 2^390*G */ 2096 for (i = 1; i <= 4; i <<= 1) { 2097 point_double(pre->g_pre_comp[2 * i][0], pre->g_pre_comp[2 * i][1], 2098 pre->g_pre_comp[2 * i][2], pre->g_pre_comp[i][0], 2099 pre->g_pre_comp[i][1], pre->g_pre_comp[i][2]); 2100 for (j = 0; j < 129; ++j) { 2101 point_double(pre->g_pre_comp[2 * i][0], 2102 pre->g_pre_comp[2 * i][1], 2103 pre->g_pre_comp[2 * i][2], 2104 pre->g_pre_comp[2 * i][0], 2105 pre->g_pre_comp[2 * i][1], 2106 pre->g_pre_comp[2 * i][2]); 2107 } 2108 } 2109 /* g_pre_comp[0] is the point at infinity */ 2110 memset(pre->g_pre_comp[0], 0, sizeof(pre->g_pre_comp[0])); 2111 /* the remaining multiples */ 2112 /* 2^130*G + 2^260*G */ 2113 point_add(pre->g_pre_comp[6][0], pre->g_pre_comp[6][1], 2114 pre->g_pre_comp[6][2], pre->g_pre_comp[4][0], 2115 pre->g_pre_comp[4][1], pre->g_pre_comp[4][2], 2116 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1], 2117 pre->g_pre_comp[2][2]); 2118 /* 2^130*G + 2^390*G */ 2119 point_add(pre->g_pre_comp[10][0], pre->g_pre_comp[10][1], 2120 pre->g_pre_comp[10][2], pre->g_pre_comp[8][0], 2121 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2], 2122 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1], 2123 pre->g_pre_comp[2][2]); 2124 /* 2^260*G + 2^390*G */ 2125 point_add(pre->g_pre_comp[12][0], pre->g_pre_comp[12][1], 2126 pre->g_pre_comp[12][2], pre->g_pre_comp[8][0], 2127 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2], 2128 0, pre->g_pre_comp[4][0], pre->g_pre_comp[4][1], 2129 pre->g_pre_comp[4][2]); 2130 /* 2^130*G + 2^260*G + 2^390*G */ 2131 point_add(pre->g_pre_comp[14][0], pre->g_pre_comp[14][1], 2132 pre->g_pre_comp[14][2], pre->g_pre_comp[12][0], 2133 pre->g_pre_comp[12][1], pre->g_pre_comp[12][2], 2134 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1], 2135 pre->g_pre_comp[2][2]); 2136 for (i = 1; i < 8; ++i) { 2137 /* odd multiples: add G */ 2138 point_add(pre->g_pre_comp[2 * i + 1][0], 2139 pre->g_pre_comp[2 * i + 1][1], 2140 pre->g_pre_comp[2 * i + 1][2], pre->g_pre_comp[2 * i][0], 2141 pre->g_pre_comp[2 * i][1], pre->g_pre_comp[2 * i][2], 0, 2142 pre->g_pre_comp[1][0], pre->g_pre_comp[1][1], 2143 pre->g_pre_comp[1][2]); 2144 } 2145 make_points_affine(15, &(pre->g_pre_comp[1]), tmp_felems); 2146 2147 done: 2148 SETPRECOMP(group, nistp521, pre); 2149 ret = 1; 2150 pre = NULL; 2151 err: 2152 BN_CTX_end(ctx); 2153 EC_POINT_free(generator); 2154 BN_CTX_free(new_ctx); 2155 EC_nistp521_pre_comp_free(pre); 2156 return ret; 2157 } 2158 2159 int ec_GFp_nistp521_have_precompute_mult(const EC_GROUP *group) 2160 { 2161 return HAVEPRECOMP(group, nistp521); 2162 } 2163 2164 #endif 2165