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