2 * Copyright 2011-2016 The OpenSSL Project Authors. All Rights Reserved.
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
10 /* Copyright 2011 Google Inc.
12 * Licensed under the Apache License, Version 2.0 (the "License");
14 * you may not use this file except in compliance with the License.
15 * You may obtain a copy of the License at
17 * http://www.apache.org/licenses/LICENSE-2.0
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.
27 * A 64-bit implementation of the NIST P-521 elliptic curve point multiplication
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.
34 #include <openssl/e_os2.h>
35 #ifdef OPENSSL_NO_EC_NISTP_64_GCC_128
36 NON_EMPTY_TRANSLATION_UNIT
40 # include <openssl/err.h>
43 # if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 1))
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
48 # error "Need GCC 3.1 or later to define type uint128_t"
56 * The underlying field. P521 operates over GF(2^521-1). We can serialise an
57 * element of this field into 66 bytes where the most significant byte
58 * contains only a single bit. We call this an felem_bytearray.
61 typedef u8 felem_bytearray[66];
64 * These are the parameters of P521, taken from FIPS 186-3, section D.1.2.5.
65 * These values are big-endian.
67 static const felem_bytearray nistp521_curve_params[5] = {
68 {0x01, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, /* p */
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, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
77 {0x01, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, /* a = -3 */
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, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
86 {0x00, 0x51, 0x95, 0x3e, 0xb9, 0x61, 0x8e, 0x1c, /* b */
87 0x9a, 0x1f, 0x92, 0x9a, 0x21, 0xa0, 0xb6, 0x85,
88 0x40, 0xee, 0xa2, 0xda, 0x72, 0x5b, 0x99, 0xb3,
89 0x15, 0xf3, 0xb8, 0xb4, 0x89, 0x91, 0x8e, 0xf1,
90 0x09, 0xe1, 0x56, 0x19, 0x39, 0x51, 0xec, 0x7e,
91 0x93, 0x7b, 0x16, 0x52, 0xc0, 0xbd, 0x3b, 0xb1,
92 0xbf, 0x07, 0x35, 0x73, 0xdf, 0x88, 0x3d, 0x2c,
93 0x34, 0xf1, 0xef, 0x45, 0x1f, 0xd4, 0x6b, 0x50,
95 {0x00, 0xc6, 0x85, 0x8e, 0x06, 0xb7, 0x04, 0x04, /* x */
96 0xe9, 0xcd, 0x9e, 0x3e, 0xcb, 0x66, 0x23, 0x95,
97 0xb4, 0x42, 0x9c, 0x64, 0x81, 0x39, 0x05, 0x3f,
98 0xb5, 0x21, 0xf8, 0x28, 0xaf, 0x60, 0x6b, 0x4d,
99 0x3d, 0xba, 0xa1, 0x4b, 0x5e, 0x77, 0xef, 0xe7,
100 0x59, 0x28, 0xfe, 0x1d, 0xc1, 0x27, 0xa2, 0xff,
101 0xa8, 0xde, 0x33, 0x48, 0xb3, 0xc1, 0x85, 0x6a,
102 0x42, 0x9b, 0xf9, 0x7e, 0x7e, 0x31, 0xc2, 0xe5,
104 {0x01, 0x18, 0x39, 0x29, 0x6a, 0x78, 0x9a, 0x3b, /* y */
105 0xc0, 0x04, 0x5c, 0x8a, 0x5f, 0xb4, 0x2c, 0x7d,
106 0x1b, 0xd9, 0x98, 0xf5, 0x44, 0x49, 0x57, 0x9b,
107 0x44, 0x68, 0x17, 0xaf, 0xbd, 0x17, 0x27, 0x3e,
108 0x66, 0x2c, 0x97, 0xee, 0x72, 0x99, 0x5e, 0xf4,
109 0x26, 0x40, 0xc5, 0x50, 0xb9, 0x01, 0x3f, 0xad,
110 0x07, 0x61, 0x35, 0x3c, 0x70, 0x86, 0xa2, 0x72,
111 0xc2, 0x40, 0x88, 0xbe, 0x94, 0x76, 0x9f, 0xd1,
116 * The representation of field elements.
117 * ------------------------------------
119 * We represent field elements with nine values. These values are either 64 or
120 * 128 bits and the field element represented is:
121 * v[0]*2^0 + v[1]*2^58 + v[2]*2^116 + ... + v[8]*2^464 (mod p)
122 * Each of the nine values is called a 'limb'. Since the limbs are spaced only
123 * 58 bits apart, but are greater than 58 bits in length, the most significant
124 * bits of each limb overlap with the least significant bits of the next.
126 * A field element with 64-bit limbs is an 'felem'. One with 128-bit limbs is a
131 typedef uint64_t limb;
132 typedef limb felem[NLIMBS];
133 typedef uint128_t largefelem[NLIMBS];
135 static const limb bottom57bits = 0x1ffffffffffffff;
136 static const limb bottom58bits = 0x3ffffffffffffff;
139 * bin66_to_felem takes a little-endian byte array and converts it into felem
140 * form. This assumes that the CPU is little-endian.
142 static void bin66_to_felem(felem out, const u8 in[66])
144 out[0] = (*((limb *) & in[0])) & bottom58bits;
145 out[1] = (*((limb *) & in[7]) >> 2) & bottom58bits;
146 out[2] = (*((limb *) & in[14]) >> 4) & bottom58bits;
147 out[3] = (*((limb *) & in[21]) >> 6) & bottom58bits;
148 out[4] = (*((limb *) & in[29])) & bottom58bits;
149 out[5] = (*((limb *) & in[36]) >> 2) & bottom58bits;
150 out[6] = (*((limb *) & in[43]) >> 4) & bottom58bits;
151 out[7] = (*((limb *) & in[50]) >> 6) & bottom58bits;
152 out[8] = (*((limb *) & in[58])) & bottom57bits;
156 * felem_to_bin66 takes an felem and serialises into a little endian, 66 byte
157 * array. This assumes that the CPU is little-endian.
159 static void felem_to_bin66(u8 out[66], const felem in)
162 (*((limb *) & out[0])) = in[0];
163 (*((limb *) & out[7])) |= in[1] << 2;
164 (*((limb *) & out[14])) |= in[2] << 4;
165 (*((limb *) & out[21])) |= in[3] << 6;
166 (*((limb *) & out[29])) = in[4];
167 (*((limb *) & out[36])) |= in[5] << 2;
168 (*((limb *) & out[43])) |= in[6] << 4;
169 (*((limb *) & out[50])) |= in[7] << 6;
170 (*((limb *) & out[58])) = in[8];
173 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */
174 static void flip_endian(u8 *out, const u8 *in, unsigned len)
177 for (i = 0; i < len; ++i)
178 out[i] = in[len - 1 - i];
181 /* BN_to_felem converts an OpenSSL BIGNUM into an felem */
182 static int BN_to_felem(felem out, const BIGNUM *bn)
184 felem_bytearray b_in;
185 felem_bytearray b_out;
188 /* BN_bn2bin eats leading zeroes */
189 memset(b_out, 0, sizeof(b_out));
190 num_bytes = BN_num_bytes(bn);
191 if (num_bytes > sizeof(b_out)) {
192 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
195 if (BN_is_negative(bn)) {
196 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
199 num_bytes = BN_bn2bin(bn, b_in);
200 flip_endian(b_out, b_in, num_bytes);
201 bin66_to_felem(out, b_out);
205 /* felem_to_BN converts an felem into an OpenSSL BIGNUM */
206 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
208 felem_bytearray b_in, b_out;
209 felem_to_bin66(b_in, in);
210 flip_endian(b_out, b_in, sizeof(b_out));
211 return BN_bin2bn(b_out, sizeof(b_out), out);
219 static void felem_one(felem out)
232 static void felem_assign(felem out, const felem in)
245 /* felem_sum64 sets out = out + in. */
246 static void felem_sum64(felem out, const felem in)
259 /* felem_scalar sets out = in * scalar */
260 static void felem_scalar(felem out, const felem in, limb scalar)
262 out[0] = in[0] * scalar;
263 out[1] = in[1] * scalar;
264 out[2] = in[2] * scalar;
265 out[3] = in[3] * scalar;
266 out[4] = in[4] * scalar;
267 out[5] = in[5] * scalar;
268 out[6] = in[6] * scalar;
269 out[7] = in[7] * scalar;
270 out[8] = in[8] * scalar;
273 /* felem_scalar64 sets out = out * scalar */
274 static void felem_scalar64(felem out, limb scalar)
287 /* felem_scalar128 sets out = out * scalar */
288 static void felem_scalar128(largefelem out, limb scalar)
302 * felem_neg sets |out| to |-in|
304 * in[i] < 2^59 + 2^14
308 static void felem_neg(felem out, const felem in)
310 /* In order to prevent underflow, we subtract from 0 mod p. */
311 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5);
312 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4);
314 out[0] = two62m3 - in[0];
315 out[1] = two62m2 - in[1];
316 out[2] = two62m2 - in[2];
317 out[3] = two62m2 - in[3];
318 out[4] = two62m2 - in[4];
319 out[5] = two62m2 - in[5];
320 out[6] = two62m2 - in[6];
321 out[7] = two62m2 - in[7];
322 out[8] = two62m2 - in[8];
326 * felem_diff64 subtracts |in| from |out|
328 * in[i] < 2^59 + 2^14
330 * out[i] < out[i] + 2^62
332 static void felem_diff64(felem out, const felem in)
335 * In order to prevent underflow, we add 0 mod p before subtracting.
337 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5);
338 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4);
340 out[0] += two62m3 - in[0];
341 out[1] += two62m2 - in[1];
342 out[2] += two62m2 - in[2];
343 out[3] += two62m2 - in[3];
344 out[4] += two62m2 - in[4];
345 out[5] += two62m2 - in[5];
346 out[6] += two62m2 - in[6];
347 out[7] += two62m2 - in[7];
348 out[8] += two62m2 - in[8];
352 * felem_diff_128_64 subtracts |in| from |out|
354 * in[i] < 2^62 + 2^17
356 * out[i] < out[i] + 2^63
358 static void felem_diff_128_64(largefelem out, const felem in)
361 * In order to prevent underflow, we add 0 mod p before subtracting.
363 static const limb two63m6 = (((limb) 1) << 62) - (((limb) 1) << 5);
364 static const limb two63m5 = (((limb) 1) << 62) - (((limb) 1) << 4);
366 out[0] += two63m6 - in[0];
367 out[1] += two63m5 - in[1];
368 out[2] += two63m5 - in[2];
369 out[3] += two63m5 - in[3];
370 out[4] += two63m5 - in[4];
371 out[5] += two63m5 - in[5];
372 out[6] += two63m5 - in[6];
373 out[7] += two63m5 - in[7];
374 out[8] += two63m5 - in[8];
378 * felem_diff_128_64 subtracts |in| from |out|
382 * out[i] < out[i] + 2^127 - 2^69
384 static void felem_diff128(largefelem out, const largefelem in)
387 * In order to prevent underflow, we add 0 mod p before subtracting.
389 static const uint128_t two127m70 =
390 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 70);
391 static const uint128_t two127m69 =
392 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 69);
394 out[0] += (two127m70 - in[0]);
395 out[1] += (two127m69 - in[1]);
396 out[2] += (two127m69 - in[2]);
397 out[3] += (two127m69 - in[3]);
398 out[4] += (two127m69 - in[4]);
399 out[5] += (two127m69 - in[5]);
400 out[6] += (two127m69 - in[6]);
401 out[7] += (two127m69 - in[7]);
402 out[8] += (two127m69 - in[8]);
406 * felem_square sets |out| = |in|^2
410 * out[i] < 17 * max(in[i]) * max(in[i])
412 static void felem_square(largefelem out, const felem in)
415 felem_scalar(inx2, in, 2);
416 felem_scalar(inx4, in, 4);
419 * We have many cases were we want to do
422 * This is obviously just
424 * However, rather than do the doubling on the 128 bit result, we
425 * double one of the inputs to the multiplication by reading from
429 out[0] = ((uint128_t) in[0]) * in[0];
430 out[1] = ((uint128_t) in[0]) * inx2[1];
431 out[2] = ((uint128_t) in[0]) * inx2[2] + ((uint128_t) in[1]) * in[1];
432 out[3] = ((uint128_t) in[0]) * inx2[3] + ((uint128_t) in[1]) * inx2[2];
433 out[4] = ((uint128_t) in[0]) * inx2[4] +
434 ((uint128_t) in[1]) * inx2[3] + ((uint128_t) in[2]) * in[2];
435 out[5] = ((uint128_t) in[0]) * inx2[5] +
436 ((uint128_t) in[1]) * inx2[4] + ((uint128_t) in[2]) * inx2[3];
437 out[6] = ((uint128_t) in[0]) * inx2[6] +
438 ((uint128_t) in[1]) * inx2[5] +
439 ((uint128_t) in[2]) * inx2[4] + ((uint128_t) in[3]) * in[3];
440 out[7] = ((uint128_t) in[0]) * inx2[7] +
441 ((uint128_t) in[1]) * inx2[6] +
442 ((uint128_t) in[2]) * inx2[5] + ((uint128_t) in[3]) * inx2[4];
443 out[8] = ((uint128_t) in[0]) * inx2[8] +
444 ((uint128_t) in[1]) * inx2[7] +
445 ((uint128_t) in[2]) * inx2[6] +
446 ((uint128_t) in[3]) * inx2[5] + ((uint128_t) in[4]) * in[4];
449 * The remaining limbs fall above 2^521, with the first falling at 2^522.
450 * They correspond to locations one bit up from the limbs produced above
451 * so we would have to multiply by two to align them. Again, rather than
452 * operate on the 128-bit result, we double one of the inputs to the
453 * multiplication. If we want to double for both this reason, and the
454 * reason above, then we end up multiplying by four.
458 out[0] += ((uint128_t) in[1]) * inx4[8] +
459 ((uint128_t) in[2]) * inx4[7] +
460 ((uint128_t) in[3]) * inx4[6] + ((uint128_t) in[4]) * inx4[5];
463 out[1] += ((uint128_t) in[2]) * inx4[8] +
464 ((uint128_t) in[3]) * inx4[7] +
465 ((uint128_t) in[4]) * inx4[6] + ((uint128_t) in[5]) * inx2[5];
468 out[2] += ((uint128_t) in[3]) * inx4[8] +
469 ((uint128_t) in[4]) * inx4[7] + ((uint128_t) in[5]) * inx4[6];
472 out[3] += ((uint128_t) in[4]) * inx4[8] +
473 ((uint128_t) in[5]) * inx4[7] + ((uint128_t) in[6]) * inx2[6];
476 out[4] += ((uint128_t) in[5]) * inx4[8] + ((uint128_t) in[6]) * inx4[7];
479 out[5] += ((uint128_t) in[6]) * inx4[8] + ((uint128_t) in[7]) * inx2[7];
482 out[6] += ((uint128_t) in[7]) * inx4[8];
485 out[7] += ((uint128_t) in[8]) * inx2[8];
489 * felem_mul sets |out| = |in1| * |in2|
494 * out[i] < 17 * max(in1[i]) * max(in2[i])
496 static void felem_mul(largefelem out, const felem in1, const felem in2)
499 felem_scalar(in2x2, in2, 2);
501 out[0] = ((uint128_t) in1[0]) * in2[0];
503 out[1] = ((uint128_t) in1[0]) * in2[1] +
504 ((uint128_t) in1[1]) * in2[0];
506 out[2] = ((uint128_t) in1[0]) * in2[2] +
507 ((uint128_t) in1[1]) * in2[1] +
508 ((uint128_t) in1[2]) * in2[0];
510 out[3] = ((uint128_t) in1[0]) * in2[3] +
511 ((uint128_t) in1[1]) * in2[2] +
512 ((uint128_t) in1[2]) * in2[1] +
513 ((uint128_t) in1[3]) * in2[0];
515 out[4] = ((uint128_t) in1[0]) * in2[4] +
516 ((uint128_t) in1[1]) * in2[3] +
517 ((uint128_t) in1[2]) * in2[2] +
518 ((uint128_t) in1[3]) * in2[1] +
519 ((uint128_t) in1[4]) * in2[0];
521 out[5] = ((uint128_t) in1[0]) * in2[5] +
522 ((uint128_t) in1[1]) * in2[4] +
523 ((uint128_t) in1[2]) * in2[3] +
524 ((uint128_t) in1[3]) * in2[2] +
525 ((uint128_t) in1[4]) * in2[1] +
526 ((uint128_t) in1[5]) * in2[0];
528 out[6] = ((uint128_t) in1[0]) * in2[6] +
529 ((uint128_t) in1[1]) * in2[5] +
530 ((uint128_t) in1[2]) * in2[4] +
531 ((uint128_t) in1[3]) * in2[3] +
532 ((uint128_t) in1[4]) * in2[2] +
533 ((uint128_t) in1[5]) * in2[1] +
534 ((uint128_t) in1[6]) * in2[0];
536 out[7] = ((uint128_t) in1[0]) * in2[7] +
537 ((uint128_t) in1[1]) * in2[6] +
538 ((uint128_t) in1[2]) * in2[5] +
539 ((uint128_t) in1[3]) * in2[4] +
540 ((uint128_t) in1[4]) * in2[3] +
541 ((uint128_t) in1[5]) * in2[2] +
542 ((uint128_t) in1[6]) * in2[1] +
543 ((uint128_t) in1[7]) * in2[0];
545 out[8] = ((uint128_t) in1[0]) * in2[8] +
546 ((uint128_t) in1[1]) * in2[7] +
547 ((uint128_t) in1[2]) * in2[6] +
548 ((uint128_t) in1[3]) * in2[5] +
549 ((uint128_t) in1[4]) * in2[4] +
550 ((uint128_t) in1[5]) * in2[3] +
551 ((uint128_t) in1[6]) * in2[2] +
552 ((uint128_t) in1[7]) * in2[1] +
553 ((uint128_t) in1[8]) * in2[0];
555 /* See comment in felem_square about the use of in2x2 here */
557 out[0] += ((uint128_t) in1[1]) * in2x2[8] +
558 ((uint128_t) in1[2]) * in2x2[7] +
559 ((uint128_t) in1[3]) * in2x2[6] +
560 ((uint128_t) in1[4]) * in2x2[5] +
561 ((uint128_t) in1[5]) * in2x2[4] +
562 ((uint128_t) in1[6]) * in2x2[3] +
563 ((uint128_t) in1[7]) * in2x2[2] +
564 ((uint128_t) in1[8]) * in2x2[1];
566 out[1] += ((uint128_t) in1[2]) * in2x2[8] +
567 ((uint128_t) in1[3]) * in2x2[7] +
568 ((uint128_t) in1[4]) * in2x2[6] +
569 ((uint128_t) in1[5]) * in2x2[5] +
570 ((uint128_t) in1[6]) * in2x2[4] +
571 ((uint128_t) in1[7]) * in2x2[3] +
572 ((uint128_t) in1[8]) * in2x2[2];
574 out[2] += ((uint128_t) in1[3]) * in2x2[8] +
575 ((uint128_t) in1[4]) * in2x2[7] +
576 ((uint128_t) in1[5]) * in2x2[6] +
577 ((uint128_t) in1[6]) * in2x2[5] +
578 ((uint128_t) in1[7]) * in2x2[4] +
579 ((uint128_t) in1[8]) * in2x2[3];
581 out[3] += ((uint128_t) in1[4]) * in2x2[8] +
582 ((uint128_t) in1[5]) * in2x2[7] +
583 ((uint128_t) in1[6]) * in2x2[6] +
584 ((uint128_t) in1[7]) * in2x2[5] +
585 ((uint128_t) in1[8]) * in2x2[4];
587 out[4] += ((uint128_t) in1[5]) * in2x2[8] +
588 ((uint128_t) in1[6]) * in2x2[7] +
589 ((uint128_t) in1[7]) * in2x2[6] +
590 ((uint128_t) in1[8]) * in2x2[5];
592 out[5] += ((uint128_t) in1[6]) * in2x2[8] +
593 ((uint128_t) in1[7]) * in2x2[7] +
594 ((uint128_t) in1[8]) * in2x2[6];
596 out[6] += ((uint128_t) in1[7]) * in2x2[8] +
597 ((uint128_t) in1[8]) * in2x2[7];
599 out[7] += ((uint128_t) in1[8]) * in2x2[8];
602 static const limb bottom52bits = 0xfffffffffffff;
605 * felem_reduce converts a largefelem to an felem.
609 * out[i] < 2^59 + 2^14
611 static void felem_reduce(felem out, const largefelem in)
613 u64 overflow1, overflow2;
615 out[0] = ((limb) in[0]) & bottom58bits;
616 out[1] = ((limb) in[1]) & bottom58bits;
617 out[2] = ((limb) in[2]) & bottom58bits;
618 out[3] = ((limb) in[3]) & bottom58bits;
619 out[4] = ((limb) in[4]) & bottom58bits;
620 out[5] = ((limb) in[5]) & bottom58bits;
621 out[6] = ((limb) in[6]) & bottom58bits;
622 out[7] = ((limb) in[7]) & bottom58bits;
623 out[8] = ((limb) in[8]) & bottom58bits;
627 out[1] += ((limb) in[0]) >> 58;
628 out[1] += (((limb) (in[0] >> 64)) & bottom52bits) << 6;
630 * out[1] < 2^58 + 2^6 + 2^58
633 out[2] += ((limb) (in[0] >> 64)) >> 52;
635 out[2] += ((limb) in[1]) >> 58;
636 out[2] += (((limb) (in[1] >> 64)) & bottom52bits) << 6;
637 out[3] += ((limb) (in[1] >> 64)) >> 52;
639 out[3] += ((limb) in[2]) >> 58;
640 out[3] += (((limb) (in[2] >> 64)) & bottom52bits) << 6;
641 out[4] += ((limb) (in[2] >> 64)) >> 52;
643 out[4] += ((limb) in[3]) >> 58;
644 out[4] += (((limb) (in[3] >> 64)) & bottom52bits) << 6;
645 out[5] += ((limb) (in[3] >> 64)) >> 52;
647 out[5] += ((limb) in[4]) >> 58;
648 out[5] += (((limb) (in[4] >> 64)) & bottom52bits) << 6;
649 out[6] += ((limb) (in[4] >> 64)) >> 52;
651 out[6] += ((limb) in[5]) >> 58;
652 out[6] += (((limb) (in[5] >> 64)) & bottom52bits) << 6;
653 out[7] += ((limb) (in[5] >> 64)) >> 52;
655 out[7] += ((limb) in[6]) >> 58;
656 out[7] += (((limb) (in[6] >> 64)) & bottom52bits) << 6;
657 out[8] += ((limb) (in[6] >> 64)) >> 52;
659 out[8] += ((limb) in[7]) >> 58;
660 out[8] += (((limb) (in[7] >> 64)) & bottom52bits) << 6;
662 * out[x > 1] < 2^58 + 2^6 + 2^58 + 2^12
665 overflow1 = ((limb) (in[7] >> 64)) >> 52;
667 overflow1 += ((limb) in[8]) >> 58;
668 overflow1 += (((limb) (in[8] >> 64)) & bottom52bits) << 6;
669 overflow2 = ((limb) (in[8] >> 64)) >> 52;
671 overflow1 <<= 1; /* overflow1 < 2^13 + 2^7 + 2^59 */
672 overflow2 <<= 1; /* overflow2 < 2^13 */
674 out[0] += overflow1; /* out[0] < 2^60 */
675 out[1] += overflow2; /* out[1] < 2^59 + 2^6 + 2^13 */
677 out[1] += out[0] >> 58;
678 out[0] &= bottom58bits;
681 * out[1] < 2^59 + 2^6 + 2^13 + 2^2
686 static void felem_square_reduce(felem out, const felem in)
689 felem_square(tmp, in);
690 felem_reduce(out, tmp);
693 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
696 felem_mul(tmp, in1, in2);
697 felem_reduce(out, tmp);
701 * felem_inv calculates |out| = |in|^{-1}
703 * Based on Fermat's Little Theorem:
705 * a^{p-1} = 1 (mod p)
706 * a^{p-2} = a^{-1} (mod p)
708 static void felem_inv(felem out, const felem in)
710 felem ftmp, ftmp2, ftmp3, ftmp4;
714 felem_square(tmp, in);
715 felem_reduce(ftmp, tmp); /* 2^1 */
716 felem_mul(tmp, in, ftmp);
717 felem_reduce(ftmp, tmp); /* 2^2 - 2^0 */
718 felem_assign(ftmp2, ftmp);
719 felem_square(tmp, ftmp);
720 felem_reduce(ftmp, tmp); /* 2^3 - 2^1 */
721 felem_mul(tmp, in, ftmp);
722 felem_reduce(ftmp, tmp); /* 2^3 - 2^0 */
723 felem_square(tmp, ftmp);
724 felem_reduce(ftmp, tmp); /* 2^4 - 2^1 */
726 felem_square(tmp, ftmp2);
727 felem_reduce(ftmp3, tmp); /* 2^3 - 2^1 */
728 felem_square(tmp, ftmp3);
729 felem_reduce(ftmp3, tmp); /* 2^4 - 2^2 */
730 felem_mul(tmp, ftmp3, ftmp2);
731 felem_reduce(ftmp3, tmp); /* 2^4 - 2^0 */
733 felem_assign(ftmp2, ftmp3);
734 felem_square(tmp, ftmp3);
735 felem_reduce(ftmp3, tmp); /* 2^5 - 2^1 */
736 felem_square(tmp, ftmp3);
737 felem_reduce(ftmp3, tmp); /* 2^6 - 2^2 */
738 felem_square(tmp, ftmp3);
739 felem_reduce(ftmp3, tmp); /* 2^7 - 2^3 */
740 felem_square(tmp, ftmp3);
741 felem_reduce(ftmp3, tmp); /* 2^8 - 2^4 */
742 felem_assign(ftmp4, ftmp3);
743 felem_mul(tmp, ftmp3, ftmp);
744 felem_reduce(ftmp4, tmp); /* 2^8 - 2^1 */
745 felem_square(tmp, ftmp4);
746 felem_reduce(ftmp4, tmp); /* 2^9 - 2^2 */
747 felem_mul(tmp, ftmp3, ftmp2);
748 felem_reduce(ftmp3, tmp); /* 2^8 - 2^0 */
749 felem_assign(ftmp2, ftmp3);
751 for (i = 0; i < 8; i++) {
752 felem_square(tmp, ftmp3);
753 felem_reduce(ftmp3, tmp); /* 2^16 - 2^8 */
755 felem_mul(tmp, ftmp3, ftmp2);
756 felem_reduce(ftmp3, tmp); /* 2^16 - 2^0 */
757 felem_assign(ftmp2, ftmp3);
759 for (i = 0; i < 16; i++) {
760 felem_square(tmp, ftmp3);
761 felem_reduce(ftmp3, tmp); /* 2^32 - 2^16 */
763 felem_mul(tmp, ftmp3, ftmp2);
764 felem_reduce(ftmp3, tmp); /* 2^32 - 2^0 */
765 felem_assign(ftmp2, ftmp3);
767 for (i = 0; i < 32; i++) {
768 felem_square(tmp, ftmp3);
769 felem_reduce(ftmp3, tmp); /* 2^64 - 2^32 */
771 felem_mul(tmp, ftmp3, ftmp2);
772 felem_reduce(ftmp3, tmp); /* 2^64 - 2^0 */
773 felem_assign(ftmp2, ftmp3);
775 for (i = 0; i < 64; i++) {
776 felem_square(tmp, ftmp3);
777 felem_reduce(ftmp3, tmp); /* 2^128 - 2^64 */
779 felem_mul(tmp, ftmp3, ftmp2);
780 felem_reduce(ftmp3, tmp); /* 2^128 - 2^0 */
781 felem_assign(ftmp2, ftmp3);
783 for (i = 0; i < 128; i++) {
784 felem_square(tmp, ftmp3);
785 felem_reduce(ftmp3, tmp); /* 2^256 - 2^128 */
787 felem_mul(tmp, ftmp3, ftmp2);
788 felem_reduce(ftmp3, tmp); /* 2^256 - 2^0 */
789 felem_assign(ftmp2, ftmp3);
791 for (i = 0; i < 256; i++) {
792 felem_square(tmp, ftmp3);
793 felem_reduce(ftmp3, tmp); /* 2^512 - 2^256 */
795 felem_mul(tmp, ftmp3, ftmp2);
796 felem_reduce(ftmp3, tmp); /* 2^512 - 2^0 */
798 for (i = 0; i < 9; i++) {
799 felem_square(tmp, ftmp3);
800 felem_reduce(ftmp3, tmp); /* 2^521 - 2^9 */
802 felem_mul(tmp, ftmp3, ftmp4);
803 felem_reduce(ftmp3, tmp); /* 2^512 - 2^2 */
804 felem_mul(tmp, ftmp3, in);
805 felem_reduce(out, tmp); /* 2^512 - 3 */
808 /* This is 2^521-1, expressed as an felem */
809 static const felem kPrime = {
810 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff,
811 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff,
812 0x03ffffffffffffff, 0x03ffffffffffffff, 0x01ffffffffffffff
816 * felem_is_zero returns a limb with all bits set if |in| == 0 (mod p) and 0
819 * in[i] < 2^59 + 2^14
821 static limb felem_is_zero(const felem in)
825 felem_assign(ftmp, in);
827 ftmp[0] += ftmp[8] >> 57;
828 ftmp[8] &= bottom57bits;
830 ftmp[1] += ftmp[0] >> 58;
831 ftmp[0] &= bottom58bits;
832 ftmp[2] += ftmp[1] >> 58;
833 ftmp[1] &= bottom58bits;
834 ftmp[3] += ftmp[2] >> 58;
835 ftmp[2] &= bottom58bits;
836 ftmp[4] += ftmp[3] >> 58;
837 ftmp[3] &= bottom58bits;
838 ftmp[5] += ftmp[4] >> 58;
839 ftmp[4] &= bottom58bits;
840 ftmp[6] += ftmp[5] >> 58;
841 ftmp[5] &= bottom58bits;
842 ftmp[7] += ftmp[6] >> 58;
843 ftmp[6] &= bottom58bits;
844 ftmp[8] += ftmp[7] >> 58;
845 ftmp[7] &= bottom58bits;
846 /* ftmp[8] < 2^57 + 4 */
849 * The ninth limb of 2*(2^521-1) is 0x03ffffffffffffff, which is greater
850 * than our bound for ftmp[8]. Therefore we only have to check if the
851 * zero is zero or 2^521-1.
867 * We know that ftmp[i] < 2^63, therefore the only way that the top bit
868 * can be set is if is_zero was 0 before the decrement.
870 is_zero = ((s64) is_zero) >> 63;
872 is_p = ftmp[0] ^ kPrime[0];
873 is_p |= ftmp[1] ^ kPrime[1];
874 is_p |= ftmp[2] ^ kPrime[2];
875 is_p |= ftmp[3] ^ kPrime[3];
876 is_p |= ftmp[4] ^ kPrime[4];
877 is_p |= ftmp[5] ^ kPrime[5];
878 is_p |= ftmp[6] ^ kPrime[6];
879 is_p |= ftmp[7] ^ kPrime[7];
880 is_p |= ftmp[8] ^ kPrime[8];
883 is_p = ((s64) is_p) >> 63;
889 static int felem_is_zero_int(const void *in)
891 return (int)(felem_is_zero(in) & ((limb) 1));
895 * felem_contract converts |in| to its unique, minimal representation.
897 * in[i] < 2^59 + 2^14
899 static void felem_contract(felem out, const felem in)
901 limb is_p, is_greater, sign;
902 static const limb two58 = ((limb) 1) << 58;
904 felem_assign(out, in);
906 out[0] += out[8] >> 57;
907 out[8] &= bottom57bits;
909 out[1] += out[0] >> 58;
910 out[0] &= bottom58bits;
911 out[2] += out[1] >> 58;
912 out[1] &= bottom58bits;
913 out[3] += out[2] >> 58;
914 out[2] &= bottom58bits;
915 out[4] += out[3] >> 58;
916 out[3] &= bottom58bits;
917 out[5] += out[4] >> 58;
918 out[4] &= bottom58bits;
919 out[6] += out[5] >> 58;
920 out[5] &= bottom58bits;
921 out[7] += out[6] >> 58;
922 out[6] &= bottom58bits;
923 out[8] += out[7] >> 58;
924 out[7] &= bottom58bits;
925 /* out[8] < 2^57 + 4 */
928 * If the value is greater than 2^521-1 then we have to subtract 2^521-1
929 * out. See the comments in felem_is_zero regarding why we don't test for
930 * other multiples of the prime.
934 * First, if |out| is equal to 2^521-1, we subtract it out to get zero.
937 is_p = out[0] ^ kPrime[0];
938 is_p |= out[1] ^ kPrime[1];
939 is_p |= out[2] ^ kPrime[2];
940 is_p |= out[3] ^ kPrime[3];
941 is_p |= out[4] ^ kPrime[4];
942 is_p |= out[5] ^ kPrime[5];
943 is_p |= out[6] ^ kPrime[6];
944 is_p |= out[7] ^ kPrime[7];
945 is_p |= out[8] ^ kPrime[8];
954 is_p = ((s64) is_p) >> 63;
957 /* is_p is 0 iff |out| == 2^521-1 and all ones otherwise */
970 * In order to test that |out| >= 2^521-1 we need only test if out[8] >>
971 * 57 is greater than zero as (2^521-1) + x >= 2^522
973 is_greater = out[8] >> 57;
974 is_greater |= is_greater << 32;
975 is_greater |= is_greater << 16;
976 is_greater |= is_greater << 8;
977 is_greater |= is_greater << 4;
978 is_greater |= is_greater << 2;
979 is_greater |= is_greater << 1;
980 is_greater = ((s64) is_greater) >> 63;
982 out[0] -= kPrime[0] & is_greater;
983 out[1] -= kPrime[1] & is_greater;
984 out[2] -= kPrime[2] & is_greater;
985 out[3] -= kPrime[3] & is_greater;
986 out[4] -= kPrime[4] & is_greater;
987 out[5] -= kPrime[5] & is_greater;
988 out[6] -= kPrime[6] & is_greater;
989 out[7] -= kPrime[7] & is_greater;
990 out[8] -= kPrime[8] & is_greater;
992 /* Eliminate negative coefficients */
993 sign = -(out[0] >> 63);
994 out[0] += (two58 & sign);
995 out[1] -= (1 & sign);
996 sign = -(out[1] >> 63);
997 out[1] += (two58 & sign);
998 out[2] -= (1 & sign);
999 sign = -(out[2] >> 63);
1000 out[2] += (two58 & sign);
1001 out[3] -= (1 & sign);
1002 sign = -(out[3] >> 63);
1003 out[3] += (two58 & sign);
1004 out[4] -= (1 & sign);
1005 sign = -(out[4] >> 63);
1006 out[4] += (two58 & sign);
1007 out[5] -= (1 & sign);
1008 sign = -(out[0] >> 63);
1009 out[5] += (two58 & sign);
1010 out[6] -= (1 & sign);
1011 sign = -(out[6] >> 63);
1012 out[6] += (two58 & sign);
1013 out[7] -= (1 & sign);
1014 sign = -(out[7] >> 63);
1015 out[7] += (two58 & sign);
1016 out[8] -= (1 & sign);
1017 sign = -(out[5] >> 63);
1018 out[5] += (two58 & sign);
1019 out[6] -= (1 & sign);
1020 sign = -(out[6] >> 63);
1021 out[6] += (two58 & sign);
1022 out[7] -= (1 & sign);
1023 sign = -(out[7] >> 63);
1024 out[7] += (two58 & sign);
1025 out[8] -= (1 & sign);
1032 * Building on top of the field operations we have the operations on the
1033 * elliptic curve group itself. Points on the curve are represented in Jacobian
1037 * point_double calculates 2*(x_in, y_in, z_in)
1039 * The method is taken from:
1040 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
1042 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed.
1043 * while x_out == y_in is not (maybe this works, but it's not tested). */
1045 point_double(felem x_out, felem y_out, felem z_out,
1046 const felem x_in, const felem y_in, const felem z_in)
1048 largefelem tmp, tmp2;
1049 felem delta, gamma, beta, alpha, ftmp, ftmp2;
1051 felem_assign(ftmp, x_in);
1052 felem_assign(ftmp2, x_in);
1055 felem_square(tmp, z_in);
1056 felem_reduce(delta, tmp); /* delta[i] < 2^59 + 2^14 */
1059 felem_square(tmp, y_in);
1060 felem_reduce(gamma, tmp); /* gamma[i] < 2^59 + 2^14 */
1062 /* beta = x*gamma */
1063 felem_mul(tmp, x_in, gamma);
1064 felem_reduce(beta, tmp); /* beta[i] < 2^59 + 2^14 */
1066 /* alpha = 3*(x-delta)*(x+delta) */
1067 felem_diff64(ftmp, delta);
1068 /* ftmp[i] < 2^61 */
1069 felem_sum64(ftmp2, delta);
1070 /* ftmp2[i] < 2^60 + 2^15 */
1071 felem_scalar64(ftmp2, 3);
1072 /* ftmp2[i] < 3*2^60 + 3*2^15 */
1073 felem_mul(tmp, ftmp, ftmp2);
1075 * tmp[i] < 17(3*2^121 + 3*2^76)
1076 * = 61*2^121 + 61*2^76
1077 * < 64*2^121 + 64*2^76
1081 felem_reduce(alpha, tmp);
1083 /* x' = alpha^2 - 8*beta */
1084 felem_square(tmp, alpha);
1086 * tmp[i] < 17*2^120 < 2^125
1088 felem_assign(ftmp, beta);
1089 felem_scalar64(ftmp, 8);
1090 /* ftmp[i] < 2^62 + 2^17 */
1091 felem_diff_128_64(tmp, ftmp);
1092 /* tmp[i] < 2^125 + 2^63 + 2^62 + 2^17 */
1093 felem_reduce(x_out, tmp);
1095 /* z' = (y + z)^2 - gamma - delta */
1096 felem_sum64(delta, gamma);
1097 /* delta[i] < 2^60 + 2^15 */
1098 felem_assign(ftmp, y_in);
1099 felem_sum64(ftmp, z_in);
1100 /* ftmp[i] < 2^60 + 2^15 */
1101 felem_square(tmp, ftmp);
1103 * tmp[i] < 17(2^122) < 2^127
1105 felem_diff_128_64(tmp, delta);
1106 /* tmp[i] < 2^127 + 2^63 */
1107 felem_reduce(z_out, tmp);
1109 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
1110 felem_scalar64(beta, 4);
1111 /* beta[i] < 2^61 + 2^16 */
1112 felem_diff64(beta, x_out);
1113 /* beta[i] < 2^61 + 2^60 + 2^16 */
1114 felem_mul(tmp, alpha, beta);
1116 * tmp[i] < 17*((2^59 + 2^14)(2^61 + 2^60 + 2^16))
1117 * = 17*(2^120 + 2^75 + 2^119 + 2^74 + 2^75 + 2^30)
1118 * = 17*(2^120 + 2^119 + 2^76 + 2^74 + 2^30)
1121 felem_square(tmp2, gamma);
1123 * tmp2[i] < 17*(2^59 + 2^14)^2
1124 * = 17*(2^118 + 2^74 + 2^28)
1126 felem_scalar128(tmp2, 8);
1128 * tmp2[i] < 8*17*(2^118 + 2^74 + 2^28)
1129 * = 2^125 + 2^121 + 2^81 + 2^77 + 2^35 + 2^31
1132 felem_diff128(tmp, tmp2);
1134 * tmp[i] < 2^127 - 2^69 + 17(2^120 + 2^119 + 2^76 + 2^74 + 2^30)
1135 * = 2^127 + 2^124 + 2^122 + 2^120 + 2^118 + 2^80 + 2^78 + 2^76 +
1136 * 2^74 + 2^69 + 2^34 + 2^30
1139 felem_reduce(y_out, tmp);
1142 /* copy_conditional copies in to out iff mask is all ones. */
1143 static void copy_conditional(felem out, const felem in, limb mask)
1146 for (i = 0; i < NLIMBS; ++i) {
1147 const limb tmp = mask & (in[i] ^ out[i]);
1153 * point_add calculates (x1, y1, z1) + (x2, y2, z2)
1155 * The method is taken from
1156 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl,
1157 * adapted for mixed addition (z2 = 1, or z2 = 0 for the point at infinity).
1159 * This function includes a branch for checking whether the two input points
1160 * are equal (while not equal to the point at infinity). See comment below
1163 static void point_add(felem x3, felem y3, felem z3,
1164 const felem x1, const felem y1, const felem z1,
1165 const int mixed, const felem x2, const felem y2,
1168 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, ftmp6, x_out, y_out, z_out;
1169 largefelem tmp, tmp2;
1170 limb x_equal, y_equal, z1_is_zero, z2_is_zero;
1172 z1_is_zero = felem_is_zero(z1);
1173 z2_is_zero = felem_is_zero(z2);
1175 /* ftmp = z1z1 = z1**2 */
1176 felem_square(tmp, z1);
1177 felem_reduce(ftmp, tmp);
1180 /* ftmp2 = z2z2 = z2**2 */
1181 felem_square(tmp, z2);
1182 felem_reduce(ftmp2, tmp);
1184 /* u1 = ftmp3 = x1*z2z2 */
1185 felem_mul(tmp, x1, ftmp2);
1186 felem_reduce(ftmp3, tmp);
1188 /* ftmp5 = z1 + z2 */
1189 felem_assign(ftmp5, z1);
1190 felem_sum64(ftmp5, z2);
1191 /* ftmp5[i] < 2^61 */
1193 /* ftmp5 = (z1 + z2)**2 - z1z1 - z2z2 = 2*z1z2 */
1194 felem_square(tmp, ftmp5);
1195 /* tmp[i] < 17*2^122 */
1196 felem_diff_128_64(tmp, ftmp);
1197 /* tmp[i] < 17*2^122 + 2^63 */
1198 felem_diff_128_64(tmp, ftmp2);
1199 /* tmp[i] < 17*2^122 + 2^64 */
1200 felem_reduce(ftmp5, tmp);
1202 /* ftmp2 = z2 * z2z2 */
1203 felem_mul(tmp, ftmp2, z2);
1204 felem_reduce(ftmp2, tmp);
1206 /* s1 = ftmp6 = y1 * z2**3 */
1207 felem_mul(tmp, y1, ftmp2);
1208 felem_reduce(ftmp6, tmp);
1211 * We'll assume z2 = 1 (special case z2 = 0 is handled later)
1214 /* u1 = ftmp3 = x1*z2z2 */
1215 felem_assign(ftmp3, x1);
1217 /* ftmp5 = 2*z1z2 */
1218 felem_scalar(ftmp5, z1, 2);
1220 /* s1 = ftmp6 = y1 * z2**3 */
1221 felem_assign(ftmp6, y1);
1225 felem_mul(tmp, x2, ftmp);
1226 /* tmp[i] < 17*2^120 */
1228 /* h = ftmp4 = u2 - u1 */
1229 felem_diff_128_64(tmp, ftmp3);
1230 /* tmp[i] < 17*2^120 + 2^63 */
1231 felem_reduce(ftmp4, tmp);
1233 x_equal = felem_is_zero(ftmp4);
1235 /* z_out = ftmp5 * h */
1236 felem_mul(tmp, ftmp5, ftmp4);
1237 felem_reduce(z_out, tmp);
1239 /* ftmp = z1 * z1z1 */
1240 felem_mul(tmp, ftmp, z1);
1241 felem_reduce(ftmp, tmp);
1243 /* s2 = tmp = y2 * z1**3 */
1244 felem_mul(tmp, y2, ftmp);
1245 /* tmp[i] < 17*2^120 */
1247 /* r = ftmp5 = (s2 - s1)*2 */
1248 felem_diff_128_64(tmp, ftmp6);
1249 /* tmp[i] < 17*2^120 + 2^63 */
1250 felem_reduce(ftmp5, tmp);
1251 y_equal = felem_is_zero(ftmp5);
1252 felem_scalar64(ftmp5, 2);
1253 /* ftmp5[i] < 2^61 */
1255 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
1257 * This is obviously not constant-time but it will almost-never happen
1258 * for ECDH / ECDSA. The case where it can happen is during scalar-mult
1259 * where the intermediate value gets very close to the group order.
1260 * Since |ec_GFp_nistp_recode_scalar_bits| produces signed digits for
1261 * the scalar, it's possible for the intermediate value to be a small
1262 * negative multiple of the base point, and for the final signed digit
1263 * to be the same value. We believe that this only occurs for the scalar
1264 * 1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
1265 * ffffffa51868783bf2f966b7fcc0148f709a5d03bb5c9b8899c47aebb6fb
1266 * 71e913863f7, in that case the penultimate intermediate is -9G and
1267 * the final digit is also -9G. Since this only happens for a single
1268 * scalar, the timing leak is irrelevent. (Any attacker who wanted to
1269 * check whether a secret scalar was that exact value, can already do
1272 point_double(x3, y3, z3, x1, y1, z1);
1276 /* I = ftmp = (2h)**2 */
1277 felem_assign(ftmp, ftmp4);
1278 felem_scalar64(ftmp, 2);
1279 /* ftmp[i] < 2^61 */
1280 felem_square(tmp, ftmp);
1281 /* tmp[i] < 17*2^122 */
1282 felem_reduce(ftmp, tmp);
1284 /* J = ftmp2 = h * I */
1285 felem_mul(tmp, ftmp4, ftmp);
1286 felem_reduce(ftmp2, tmp);
1288 /* V = ftmp4 = U1 * I */
1289 felem_mul(tmp, ftmp3, ftmp);
1290 felem_reduce(ftmp4, tmp);
1292 /* x_out = r**2 - J - 2V */
1293 felem_square(tmp, ftmp5);
1294 /* tmp[i] < 17*2^122 */
1295 felem_diff_128_64(tmp, ftmp2);
1296 /* tmp[i] < 17*2^122 + 2^63 */
1297 felem_assign(ftmp3, ftmp4);
1298 felem_scalar64(ftmp4, 2);
1299 /* ftmp4[i] < 2^61 */
1300 felem_diff_128_64(tmp, ftmp4);
1301 /* tmp[i] < 17*2^122 + 2^64 */
1302 felem_reduce(x_out, tmp);
1304 /* y_out = r(V-x_out) - 2 * s1 * J */
1305 felem_diff64(ftmp3, x_out);
1307 * ftmp3[i] < 2^60 + 2^60 = 2^61
1309 felem_mul(tmp, ftmp5, ftmp3);
1310 /* tmp[i] < 17*2^122 */
1311 felem_mul(tmp2, ftmp6, ftmp2);
1312 /* tmp2[i] < 17*2^120 */
1313 felem_scalar128(tmp2, 2);
1314 /* tmp2[i] < 17*2^121 */
1315 felem_diff128(tmp, tmp2);
1317 * tmp[i] < 2^127 - 2^69 + 17*2^122
1318 * = 2^126 - 2^122 - 2^6 - 2^2 - 1
1321 felem_reduce(y_out, tmp);
1323 copy_conditional(x_out, x2, z1_is_zero);
1324 copy_conditional(x_out, x1, z2_is_zero);
1325 copy_conditional(y_out, y2, z1_is_zero);
1326 copy_conditional(y_out, y1, z2_is_zero);
1327 copy_conditional(z_out, z2, z1_is_zero);
1328 copy_conditional(z_out, z1, z2_is_zero);
1329 felem_assign(x3, x_out);
1330 felem_assign(y3, y_out);
1331 felem_assign(z3, z_out);
1335 * Base point pre computation
1336 * --------------------------
1338 * Two different sorts of precomputed tables are used in the following code.
1339 * Each contain various points on the curve, where each point is three field
1340 * elements (x, y, z).
1342 * For the base point table, z is usually 1 (0 for the point at infinity).
1343 * This table has 16 elements:
1344 * index | bits | point
1345 * ------+---------+------------------------------
1348 * 2 | 0 0 1 0 | 2^130G
1349 * 3 | 0 0 1 1 | (2^130 + 1)G
1350 * 4 | 0 1 0 0 | 2^260G
1351 * 5 | 0 1 0 1 | (2^260 + 1)G
1352 * 6 | 0 1 1 0 | (2^260 + 2^130)G
1353 * 7 | 0 1 1 1 | (2^260 + 2^130 + 1)G
1354 * 8 | 1 0 0 0 | 2^390G
1355 * 9 | 1 0 0 1 | (2^390 + 1)G
1356 * 10 | 1 0 1 0 | (2^390 + 2^130)G
1357 * 11 | 1 0 1 1 | (2^390 + 2^130 + 1)G
1358 * 12 | 1 1 0 0 | (2^390 + 2^260)G
1359 * 13 | 1 1 0 1 | (2^390 + 2^260 + 1)G
1360 * 14 | 1 1 1 0 | (2^390 + 2^260 + 2^130)G
1361 * 15 | 1 1 1 1 | (2^390 + 2^260 + 2^130 + 1)G
1363 * The reason for this is so that we can clock bits into four different
1364 * locations when doing simple scalar multiplies against the base point.
1366 * Tables for other points have table[i] = iG for i in 0 .. 16. */
1368 /* gmul is the table of precomputed base points */
1369 static const felem gmul[16][3] = {
1370 {{0, 0, 0, 0, 0, 0, 0, 0, 0},
1371 {0, 0, 0, 0, 0, 0, 0, 0, 0},
1372 {0, 0, 0, 0, 0, 0, 0, 0, 0}},
1373 {{0x017e7e31c2e5bd66, 0x022cf0615a90a6fe, 0x00127a2ffa8de334,
1374 0x01dfbf9d64a3f877, 0x006b4d3dbaa14b5e, 0x014fed487e0a2bd8,
1375 0x015b4429c6481390, 0x03a73678fb2d988e, 0x00c6858e06b70404},
1376 {0x00be94769fd16650, 0x031c21a89cb09022, 0x039013fad0761353,
1377 0x02657bd099031542, 0x03273e662c97ee72, 0x01e6d11a05ebef45,
1378 0x03d1bd998f544495, 0x03001172297ed0b1, 0x011839296a789a3b},
1379 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1380 {{0x0373faacbc875bae, 0x00f325023721c671, 0x00f666fd3dbde5ad,
1381 0x01a6932363f88ea7, 0x01fc6d9e13f9c47b, 0x03bcbffc2bbf734e,
1382 0x013ee3c3647f3a92, 0x029409fefe75d07d, 0x00ef9199963d85e5},
1383 {0x011173743ad5b178, 0x02499c7c21bf7d46, 0x035beaeabb8b1a58,
1384 0x00f989c4752ea0a3, 0x0101e1de48a9c1a3, 0x01a20076be28ba6c,
1385 0x02f8052e5eb2de95, 0x01bfe8f82dea117c, 0x0160074d3c36ddb7},
1386 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1387 {{0x012f3fc373393b3b, 0x03d3d6172f1419fa, 0x02adc943c0b86873,
1388 0x00d475584177952b, 0x012a4d1673750ee2, 0x00512517a0f13b0c,
1389 0x02b184671a7b1734, 0x0315b84236f1a50a, 0x00a4afc472edbdb9},
1390 {0x00152a7077f385c4, 0x03044007d8d1c2ee, 0x0065829d61d52b52,
1391 0x00494ff6b6631d0d, 0x00a11d94d5f06bcf, 0x02d2f89474d9282e,
1392 0x0241c5727c06eeb9, 0x0386928710fbdb9d, 0x01f883f727b0dfbe},
1393 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1394 {{0x019b0c3c9185544d, 0x006243a37c9d97db, 0x02ee3cbe030a2ad2,
1395 0x00cfdd946bb51e0d, 0x0271c00932606b91, 0x03f817d1ec68c561,
1396 0x03f37009806a369c, 0x03c1f30baf184fd5, 0x01091022d6d2f065},
1397 {0x0292c583514c45ed, 0x0316fca51f9a286c, 0x00300af507c1489a,
1398 0x0295f69008298cf1, 0x02c0ed8274943d7b, 0x016509b9b47a431e,
1399 0x02bc9de9634868ce, 0x005b34929bffcb09, 0x000c1a0121681524},
1400 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1401 {{0x0286abc0292fb9f2, 0x02665eee9805b3f7, 0x01ed7455f17f26d6,
1402 0x0346355b83175d13, 0x006284944cd0a097, 0x0191895bcdec5e51,
1403 0x02e288370afda7d9, 0x03b22312bfefa67a, 0x01d104d3fc0613fe},
1404 {0x0092421a12f7e47f, 0x0077a83fa373c501, 0x03bd25c5f696bd0d,
1405 0x035c41e4d5459761, 0x01ca0d1742b24f53, 0x00aaab27863a509c,
1406 0x018b6de47df73917, 0x025c0b771705cd01, 0x01fd51d566d760a7},
1407 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1408 {{0x01dd92ff6b0d1dbd, 0x039c5e2e8f8afa69, 0x0261ed13242c3b27,
1409 0x0382c6e67026e6a0, 0x01d60b10be2089f9, 0x03c15f3dce86723f,
1410 0x03c764a32d2a062d, 0x017307eac0fad056, 0x018207c0b96c5256},
1411 {0x0196a16d60e13154, 0x03e6ce74c0267030, 0x00ddbf2b4e52a5aa,
1412 0x012738241bbf31c8, 0x00ebe8dc04685a28, 0x024c2ad6d380d4a2,
1413 0x035ee062a6e62d0e, 0x0029ed74af7d3a0f, 0x00eef32aec142ebd},
1414 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1415 {{0x00c31ec398993b39, 0x03a9f45bcda68253, 0x00ac733c24c70890,
1416 0x00872b111401ff01, 0x01d178c23195eafb, 0x03bca2c816b87f74,
1417 0x0261a9af46fbad7a, 0x0324b2a8dd3d28f9, 0x00918121d8f24e23},
1418 {0x032bc8c1ca983cd7, 0x00d869dfb08fc8c6, 0x01693cb61fce1516,
1419 0x012a5ea68f4e88a8, 0x010869cab88d7ae3, 0x009081ad277ceee1,
1420 0x033a77166d064cdc, 0x03955235a1fb3a95, 0x01251a4a9b25b65e},
1421 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1422 {{0x00148a3a1b27f40b, 0x0123186df1b31fdc, 0x00026e7beaad34ce,
1423 0x01db446ac1d3dbba, 0x0299c1a33437eaec, 0x024540610183cbb7,
1424 0x0173bb0e9ce92e46, 0x02b937e43921214b, 0x01ab0436a9bf01b5},
1425 {0x0383381640d46948, 0x008dacbf0e7f330f, 0x03602122bcc3f318,
1426 0x01ee596b200620d6, 0x03bd0585fda430b3, 0x014aed77fd123a83,
1427 0x005ace749e52f742, 0x0390fe041da2b842, 0x0189a8ceb3299242},
1428 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1429 {{0x012a19d6b3282473, 0x00c0915918b423ce, 0x023a954eb94405ae,
1430 0x00529f692be26158, 0x0289fa1b6fa4b2aa, 0x0198ae4ceea346ef,
1431 0x0047d8cdfbdedd49, 0x00cc8c8953f0f6b8, 0x001424abbff49203},
1432 {0x0256732a1115a03a, 0x0351bc38665c6733, 0x03f7b950fb4a6447,
1433 0x000afffa94c22155, 0x025763d0a4dab540, 0x000511e92d4fc283,
1434 0x030a7e9eda0ee96c, 0x004c3cd93a28bf0a, 0x017edb3a8719217f},
1435 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1436 {{0x011de5675a88e673, 0x031d7d0f5e567fbe, 0x0016b2062c970ae5,
1437 0x03f4a2be49d90aa7, 0x03cef0bd13822866, 0x03f0923dcf774a6c,
1438 0x0284bebc4f322f72, 0x016ab2645302bb2c, 0x01793f95dace0e2a},
1439 {0x010646e13527a28f, 0x01ca1babd59dc5e7, 0x01afedfd9a5595df,
1440 0x01f15785212ea6b1, 0x0324e5d64f6ae3f4, 0x02d680f526d00645,
1441 0x0127920fadf627a7, 0x03b383f75df4f684, 0x0089e0057e783b0a},
1442 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1443 {{0x00f334b9eb3c26c6, 0x0298fdaa98568dce, 0x01c2d24843a82292,
1444 0x020bcb24fa1b0711, 0x02cbdb3d2b1875e6, 0x0014907598f89422,
1445 0x03abe3aa43b26664, 0x02cbf47f720bc168, 0x0133b5e73014b79b},
1446 {0x034aab5dab05779d, 0x00cdc5d71fee9abb, 0x0399f16bd4bd9d30,
1447 0x03582fa592d82647, 0x02be1cdfb775b0e9, 0x0034f7cea32e94cb,
1448 0x0335a7f08f56f286, 0x03b707e9565d1c8b, 0x0015c946ea5b614f},
1449 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1450 {{0x024676f6cff72255, 0x00d14625cac96378, 0x00532b6008bc3767,
1451 0x01fc16721b985322, 0x023355ea1b091668, 0x029de7afdc0317c3,
1452 0x02fc8a7ca2da037c, 0x02de1217d74a6f30, 0x013f7173175b73bf},
1453 {0x0344913f441490b5, 0x0200f9e272b61eca, 0x0258a246b1dd55d2,
1454 0x03753db9ea496f36, 0x025e02937a09c5ef, 0x030cbd3d14012692,
1455 0x01793a67e70dc72a, 0x03ec1d37048a662e, 0x006550f700c32a8d},
1456 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1457 {{0x00d3f48a347eba27, 0x008e636649b61bd8, 0x00d3b93716778fb3,
1458 0x004d1915757bd209, 0x019d5311a3da44e0, 0x016d1afcbbe6aade,
1459 0x0241bf5f73265616, 0x0384672e5d50d39b, 0x005009fee522b684},
1460 {0x029b4fab064435fe, 0x018868ee095bbb07, 0x01ea3d6936cc92b8,
1461 0x000608b00f78a2f3, 0x02db911073d1c20f, 0x018205938470100a,
1462 0x01f1e4964cbe6ff2, 0x021a19a29eed4663, 0x01414485f42afa81},
1463 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1464 {{0x01612b3a17f63e34, 0x03813992885428e6, 0x022b3c215b5a9608,
1465 0x029b4057e19f2fcb, 0x0384059a587af7e6, 0x02d6400ace6fe610,
1466 0x029354d896e8e331, 0x00c047ee6dfba65e, 0x0037720542e9d49d},
1467 {0x02ce9eed7c5e9278, 0x0374ed703e79643b, 0x01316c54c4072006,
1468 0x005aaa09054b2ee8, 0x002824000c840d57, 0x03d4eba24771ed86,
1469 0x0189c50aabc3bdae, 0x0338c01541e15510, 0x00466d56e38eed42},
1470 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1471 {{0x007efd8330ad8bd6, 0x02465ed48047710b, 0x0034c6606b215e0c,
1472 0x016ae30c53cbf839, 0x01fa17bd37161216, 0x018ead4e61ce8ab9,
1473 0x005482ed5f5dee46, 0x037543755bba1d7f, 0x005e5ac7e70a9d0f},
1474 {0x0117e1bb2fdcb2a2, 0x03deea36249f40c4, 0x028d09b4a6246cb7,
1475 0x03524b8855bcf756, 0x023d7d109d5ceb58, 0x0178e43e3223ef9c,
1476 0x0154536a0c6e966a, 0x037964d1286ee9fe, 0x0199bcd90e125055},
1477 {1, 0, 0, 0, 0, 0, 0, 0, 0}}
1481 * select_point selects the |idx|th point from a precomputation table and
1484 /* pre_comp below is of the size provided in |size| */
1485 static void select_point(const limb idx, unsigned int size,
1486 const felem pre_comp[][3], felem out[3])
1489 limb *outlimbs = &out[0][0];
1491 memset(out, 0, sizeof(*out) * 3);
1493 for (i = 0; i < size; i++) {
1494 const limb *inlimbs = &pre_comp[i][0][0];
1495 limb mask = i ^ idx;
1501 for (j = 0; j < NLIMBS * 3; j++)
1502 outlimbs[j] |= inlimbs[j] & mask;
1506 /* get_bit returns the |i|th bit in |in| */
1507 static char get_bit(const felem_bytearray in, int i)
1511 return (in[i >> 3] >> (i & 7)) & 1;
1515 * Interleaved point multiplication using precomputed point multiples: The
1516 * small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[], the scalars
1517 * in scalars[]. If g_scalar is non-NULL, we also add this multiple of the
1518 * generator, using certain (large) precomputed multiples in g_pre_comp.
1519 * Output point (X, Y, Z) is stored in x_out, y_out, z_out
1521 static void batch_mul(felem x_out, felem y_out, felem z_out,
1522 const felem_bytearray scalars[],
1523 const unsigned num_points, const u8 *g_scalar,
1524 const int mixed, const felem pre_comp[][17][3],
1525 const felem g_pre_comp[16][3])
1528 unsigned num, gen_mul = (g_scalar != NULL);
1529 felem nq[3], tmp[4];
1533 /* set nq to the point at infinity */
1534 memset(nq, 0, sizeof(nq));
1537 * Loop over all scalars msb-to-lsb, interleaving additions of multiples
1538 * of the generator (last quarter of rounds) and additions of other
1539 * points multiples (every 5th round).
1541 skip = 1; /* save two point operations in the first
1543 for (i = (num_points ? 520 : 130); i >= 0; --i) {
1546 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1548 /* add multiples of the generator */
1549 if (gen_mul && (i <= 130)) {
1550 bits = get_bit(g_scalar, i + 390) << 3;
1552 bits |= get_bit(g_scalar, i + 260) << 2;
1553 bits |= get_bit(g_scalar, i + 130) << 1;
1554 bits |= get_bit(g_scalar, i);
1556 /* select the point to add, in constant time */
1557 select_point(bits, 16, g_pre_comp, tmp);
1559 /* The 1 argument below is for "mixed" */
1560 point_add(nq[0], nq[1], nq[2],
1561 nq[0], nq[1], nq[2], 1, tmp[0], tmp[1], tmp[2]);
1563 memcpy(nq, tmp, 3 * sizeof(felem));
1568 /* do other additions every 5 doublings */
1569 if (num_points && (i % 5 == 0)) {
1570 /* loop over all scalars */
1571 for (num = 0; num < num_points; ++num) {
1572 bits = get_bit(scalars[num], i + 4) << 5;
1573 bits |= get_bit(scalars[num], i + 3) << 4;
1574 bits |= get_bit(scalars[num], i + 2) << 3;
1575 bits |= get_bit(scalars[num], i + 1) << 2;
1576 bits |= get_bit(scalars[num], i) << 1;
1577 bits |= get_bit(scalars[num], i - 1);
1578 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1581 * select the point to add or subtract, in constant time
1583 select_point(digit, 17, pre_comp[num], tmp);
1584 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative
1586 copy_conditional(tmp[1], tmp[3], (-(limb) sign));
1589 point_add(nq[0], nq[1], nq[2],
1590 nq[0], nq[1], nq[2],
1591 mixed, tmp[0], tmp[1], tmp[2]);
1593 memcpy(nq, tmp, 3 * sizeof(felem));
1599 felem_assign(x_out, nq[0]);
1600 felem_assign(y_out, nq[1]);
1601 felem_assign(z_out, nq[2]);
1604 /* Precomputation for the group generator. */
1605 struct nistp521_pre_comp_st {
1606 felem g_pre_comp[16][3];
1607 CRYPTO_REF_COUNT references;
1608 CRYPTO_RWLOCK *lock;
1611 const EC_METHOD *EC_GFp_nistp521_method(void)
1613 static const EC_METHOD ret = {
1614 EC_FLAGS_DEFAULT_OCT,
1615 NID_X9_62_prime_field,
1616 ec_GFp_nistp521_group_init,
1617 ec_GFp_simple_group_finish,
1618 ec_GFp_simple_group_clear_finish,
1619 ec_GFp_nist_group_copy,
1620 ec_GFp_nistp521_group_set_curve,
1621 ec_GFp_simple_group_get_curve,
1622 ec_GFp_simple_group_get_degree,
1623 ec_group_simple_order_bits,
1624 ec_GFp_simple_group_check_discriminant,
1625 ec_GFp_simple_point_init,
1626 ec_GFp_simple_point_finish,
1627 ec_GFp_simple_point_clear_finish,
1628 ec_GFp_simple_point_copy,
1629 ec_GFp_simple_point_set_to_infinity,
1630 ec_GFp_simple_set_Jprojective_coordinates_GFp,
1631 ec_GFp_simple_get_Jprojective_coordinates_GFp,
1632 ec_GFp_simple_point_set_affine_coordinates,
1633 ec_GFp_nistp521_point_get_affine_coordinates,
1634 0 /* point_set_compressed_coordinates */ ,
1639 ec_GFp_simple_invert,
1640 ec_GFp_simple_is_at_infinity,
1641 ec_GFp_simple_is_on_curve,
1643 ec_GFp_simple_make_affine,
1644 ec_GFp_simple_points_make_affine,
1645 ec_GFp_nistp521_points_mul,
1646 ec_GFp_nistp521_precompute_mult,
1647 ec_GFp_nistp521_have_precompute_mult,
1648 ec_GFp_nist_field_mul,
1649 ec_GFp_nist_field_sqr,
1651 0 /* field_encode */ ,
1652 0 /* field_decode */ ,
1653 0, /* field_set_to_one */
1654 ec_key_simple_priv2oct,
1655 ec_key_simple_oct2priv,
1656 0, /* set private */
1657 ec_key_simple_generate_key,
1658 ec_key_simple_check_key,
1659 ec_key_simple_generate_public_key,
1662 ecdh_simple_compute_key
1668 /******************************************************************************/
1670 * FUNCTIONS TO MANAGE PRECOMPUTATION
1673 static NISTP521_PRE_COMP *nistp521_pre_comp_new()
1675 NISTP521_PRE_COMP *ret = OPENSSL_zalloc(sizeof(*ret));
1678 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1682 ret->references = 1;
1684 ret->lock = CRYPTO_THREAD_lock_new();
1685 if (ret->lock == NULL) {
1686 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1693 NISTP521_PRE_COMP *EC_nistp521_pre_comp_dup(NISTP521_PRE_COMP *p)
1697 CRYPTO_UP_REF(&p->references, &i, p->lock);
1701 void EC_nistp521_pre_comp_free(NISTP521_PRE_COMP *p)
1708 CRYPTO_DOWN_REF(&p->references, &i, p->lock);
1709 REF_PRINT_COUNT("EC_nistp521", x);
1712 REF_ASSERT_ISNT(i < 0);
1714 CRYPTO_THREAD_lock_free(p->lock);
1718 /******************************************************************************/
1720 * OPENSSL EC_METHOD FUNCTIONS
1723 int ec_GFp_nistp521_group_init(EC_GROUP *group)
1726 ret = ec_GFp_simple_group_init(group);
1727 group->a_is_minus3 = 1;
1731 int ec_GFp_nistp521_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1732 const BIGNUM *a, const BIGNUM *b,
1736 BN_CTX *new_ctx = NULL;
1737 BIGNUM *curve_p, *curve_a, *curve_b;
1740 if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1743 curve_p = BN_CTX_get(ctx);
1744 curve_a = BN_CTX_get(ctx);
1745 curve_b = BN_CTX_get(ctx);
1746 if (curve_b == NULL)
1748 BN_bin2bn(nistp521_curve_params[0], sizeof(felem_bytearray), curve_p);
1749 BN_bin2bn(nistp521_curve_params[1], sizeof(felem_bytearray), curve_a);
1750 BN_bin2bn(nistp521_curve_params[2], sizeof(felem_bytearray), curve_b);
1751 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) || (BN_cmp(curve_b, b))) {
1752 ECerr(EC_F_EC_GFP_NISTP521_GROUP_SET_CURVE,
1753 EC_R_WRONG_CURVE_PARAMETERS);
1756 group->field_mod_func = BN_nist_mod_521;
1757 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1760 BN_CTX_free(new_ctx);
1765 * Takes the Jacobian coordinates (X, Y, Z) of a point and returns (X', Y') =
1768 int ec_GFp_nistp521_point_get_affine_coordinates(const EC_GROUP *group,
1769 const EC_POINT *point,
1770 BIGNUM *x, BIGNUM *y,
1773 felem z1, z2, x_in, y_in, x_out, y_out;
1776 if (EC_POINT_is_at_infinity(group, point)) {
1777 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1778 EC_R_POINT_AT_INFINITY);
1781 if ((!BN_to_felem(x_in, point->X)) || (!BN_to_felem(y_in, point->Y)) ||
1782 (!BN_to_felem(z1, point->Z)))
1785 felem_square(tmp, z2);
1786 felem_reduce(z1, tmp);
1787 felem_mul(tmp, x_in, z1);
1788 felem_reduce(x_in, tmp);
1789 felem_contract(x_out, x_in);
1791 if (!felem_to_BN(x, x_out)) {
1792 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1797 felem_mul(tmp, z1, z2);
1798 felem_reduce(z1, tmp);
1799 felem_mul(tmp, y_in, z1);
1800 felem_reduce(y_in, tmp);
1801 felem_contract(y_out, y_in);
1803 if (!felem_to_BN(y, y_out)) {
1804 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1812 /* points below is of size |num|, and tmp_felems is of size |num+1/ */
1813 static void make_points_affine(size_t num, felem points[][3],
1817 * Runs in constant time, unless an input is the point at infinity (which
1818 * normally shouldn't happen).
1820 ec_GFp_nistp_points_make_affine_internal(num,
1824 (void (*)(void *))felem_one,
1826 (void (*)(void *, const void *))
1828 (void (*)(void *, const void *))
1829 felem_square_reduce, (void (*)
1836 (void (*)(void *, const void *))
1838 (void (*)(void *, const void *))
1843 * Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL
1844 * values Result is stored in r (r can equal one of the inputs).
1846 int ec_GFp_nistp521_points_mul(const EC_GROUP *group, EC_POINT *r,
1847 const BIGNUM *scalar, size_t num,
1848 const EC_POINT *points[],
1849 const BIGNUM *scalars[], BN_CTX *ctx)
1854 BN_CTX *new_ctx = NULL;
1855 BIGNUM *x, *y, *z, *tmp_scalar;
1856 felem_bytearray g_secret;
1857 felem_bytearray *secrets = NULL;
1858 felem (*pre_comp)[17][3] = NULL;
1859 felem *tmp_felems = NULL;
1860 felem_bytearray tmp;
1861 unsigned i, num_bytes;
1862 int have_pre_comp = 0;
1863 size_t num_points = num;
1864 felem x_in, y_in, z_in, x_out, y_out, z_out;
1865 NISTP521_PRE_COMP *pre = NULL;
1866 felem(*g_pre_comp)[3] = NULL;
1867 EC_POINT *generator = NULL;
1868 const EC_POINT *p = NULL;
1869 const BIGNUM *p_scalar = NULL;
1872 if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1875 x = BN_CTX_get(ctx);
1876 y = BN_CTX_get(ctx);
1877 z = BN_CTX_get(ctx);
1878 tmp_scalar = BN_CTX_get(ctx);
1879 if (tmp_scalar == NULL)
1882 if (scalar != NULL) {
1883 pre = group->pre_comp.nistp521;
1885 /* we have precomputation, try to use it */
1886 g_pre_comp = &pre->g_pre_comp[0];
1888 /* try to use the standard precomputation */
1889 g_pre_comp = (felem(*)[3]) gmul;
1890 generator = EC_POINT_new(group);
1891 if (generator == NULL)
1893 /* get the generator from precomputation */
1894 if (!felem_to_BN(x, g_pre_comp[1][0]) ||
1895 !felem_to_BN(y, g_pre_comp[1][1]) ||
1896 !felem_to_BN(z, g_pre_comp[1][2])) {
1897 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1900 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1904 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1905 /* precomputation matches generator */
1909 * we don't have valid precomputation: treat the generator as a
1915 if (num_points > 0) {
1916 if (num_points >= 2) {
1918 * unless we precompute multiples for just one point, converting
1919 * those into affine form is time well spent
1923 secrets = OPENSSL_zalloc(sizeof(*secrets) * num_points);
1924 pre_comp = OPENSSL_zalloc(sizeof(*pre_comp) * num_points);
1927 OPENSSL_malloc(sizeof(*tmp_felems) * (num_points * 17 + 1));
1928 if ((secrets == NULL) || (pre_comp == NULL)
1929 || (mixed && (tmp_felems == NULL))) {
1930 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1935 * we treat NULL scalars as 0, and NULL points as points at infinity,
1936 * i.e., they contribute nothing to the linear combination
1938 for (i = 0; i < num_points; ++i) {
1941 * we didn't have a valid precomputation, so we pick the
1945 p = EC_GROUP_get0_generator(group);
1948 /* the i^th point */
1951 p_scalar = scalars[i];
1953 if ((p_scalar != NULL) && (p != NULL)) {
1954 /* reduce scalar to 0 <= scalar < 2^521 */
1955 if ((BN_num_bits(p_scalar) > 521)
1956 || (BN_is_negative(p_scalar))) {
1958 * this is an unusual input, and we don't guarantee
1961 if (!BN_nnmod(tmp_scalar, p_scalar, group->order, ctx)) {
1962 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1965 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1967 num_bytes = BN_bn2bin(p_scalar, tmp);
1968 flip_endian(secrets[i], tmp, num_bytes);
1969 /* precompute multiples */
1970 if ((!BN_to_felem(x_out, p->X)) ||
1971 (!BN_to_felem(y_out, p->Y)) ||
1972 (!BN_to_felem(z_out, p->Z)))
1974 memcpy(pre_comp[i][1][0], x_out, sizeof(felem));
1975 memcpy(pre_comp[i][1][1], y_out, sizeof(felem));
1976 memcpy(pre_comp[i][1][2], z_out, sizeof(felem));
1977 for (j = 2; j <= 16; ++j) {
1979 point_add(pre_comp[i][j][0], pre_comp[i][j][1],
1980 pre_comp[i][j][2], pre_comp[i][1][0],
1981 pre_comp[i][1][1], pre_comp[i][1][2], 0,
1982 pre_comp[i][j - 1][0],
1983 pre_comp[i][j - 1][1],
1984 pre_comp[i][j - 1][2]);
1986 point_double(pre_comp[i][j][0], pre_comp[i][j][1],
1987 pre_comp[i][j][2], pre_comp[i][j / 2][0],
1988 pre_comp[i][j / 2][1],
1989 pre_comp[i][j / 2][2]);
1995 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
1998 /* the scalar for the generator */
1999 if ((scalar != NULL) && (have_pre_comp)) {
2000 memset(g_secret, 0, sizeof(g_secret));
2001 /* reduce scalar to 0 <= scalar < 2^521 */
2002 if ((BN_num_bits(scalar) > 521) || (BN_is_negative(scalar))) {
2004 * this is an unusual input, and we don't guarantee
2007 if (!BN_nnmod(tmp_scalar, scalar, group->order, ctx)) {
2008 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
2011 num_bytes = BN_bn2bin(tmp_scalar, tmp);
2013 num_bytes = BN_bn2bin(scalar, tmp);
2014 flip_endian(g_secret, tmp, num_bytes);
2015 /* do the multiplication with generator precomputation */
2016 batch_mul(x_out, y_out, z_out,
2017 (const felem_bytearray(*))secrets, num_points,
2019 mixed, (const felem(*)[17][3])pre_comp,
2020 (const felem(*)[3])g_pre_comp);
2022 /* do the multiplication without generator precomputation */
2023 batch_mul(x_out, y_out, z_out,
2024 (const felem_bytearray(*))secrets, num_points,
2025 NULL, mixed, (const felem(*)[17][3])pre_comp, NULL);
2026 /* reduce the output to its unique minimal representation */
2027 felem_contract(x_in, x_out);
2028 felem_contract(y_in, y_out);
2029 felem_contract(z_in, z_out);
2030 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
2031 (!felem_to_BN(z, z_in))) {
2032 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
2035 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
2039 EC_POINT_free(generator);
2040 BN_CTX_free(new_ctx);
2041 OPENSSL_free(secrets);
2042 OPENSSL_free(pre_comp);
2043 OPENSSL_free(tmp_felems);
2047 int ec_GFp_nistp521_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
2050 NISTP521_PRE_COMP *pre = NULL;
2052 BN_CTX *new_ctx = NULL;
2054 EC_POINT *generator = NULL;
2055 felem tmp_felems[16];
2057 /* throw away old precomputation */
2058 EC_pre_comp_free(group);
2060 if ((ctx = new_ctx = BN_CTX_new()) == NULL)
2063 x = BN_CTX_get(ctx);
2064 y = BN_CTX_get(ctx);
2067 /* get the generator */
2068 if (group->generator == NULL)
2070 generator = EC_POINT_new(group);
2071 if (generator == NULL)
2073 BN_bin2bn(nistp521_curve_params[3], sizeof(felem_bytearray), x);
2074 BN_bin2bn(nistp521_curve_params[4], sizeof(felem_bytearray), y);
2075 if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
2077 if ((pre = nistp521_pre_comp_new()) == NULL)
2080 * if the generator is the standard one, use built-in precomputation
2082 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) {
2083 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
2086 if ((!BN_to_felem(pre->g_pre_comp[1][0], group->generator->X)) ||
2087 (!BN_to_felem(pre->g_pre_comp[1][1], group->generator->Y)) ||
2088 (!BN_to_felem(pre->g_pre_comp[1][2], group->generator->Z)))
2090 /* compute 2^130*G, 2^260*G, 2^390*G */
2091 for (i = 1; i <= 4; i <<= 1) {
2092 point_double(pre->g_pre_comp[2 * i][0], pre->g_pre_comp[2 * i][1],
2093 pre->g_pre_comp[2 * i][2], pre->g_pre_comp[i][0],
2094 pre->g_pre_comp[i][1], pre->g_pre_comp[i][2]);
2095 for (j = 0; j < 129; ++j) {
2096 point_double(pre->g_pre_comp[2 * i][0],
2097 pre->g_pre_comp[2 * i][1],
2098 pre->g_pre_comp[2 * i][2],
2099 pre->g_pre_comp[2 * i][0],
2100 pre->g_pre_comp[2 * i][1],
2101 pre->g_pre_comp[2 * i][2]);
2104 /* g_pre_comp[0] is the point at infinity */
2105 memset(pre->g_pre_comp[0], 0, sizeof(pre->g_pre_comp[0]));
2106 /* the remaining multiples */
2107 /* 2^130*G + 2^260*G */
2108 point_add(pre->g_pre_comp[6][0], pre->g_pre_comp[6][1],
2109 pre->g_pre_comp[6][2], pre->g_pre_comp[4][0],
2110 pre->g_pre_comp[4][1], pre->g_pre_comp[4][2],
2111 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2112 pre->g_pre_comp[2][2]);
2113 /* 2^130*G + 2^390*G */
2114 point_add(pre->g_pre_comp[10][0], pre->g_pre_comp[10][1],
2115 pre->g_pre_comp[10][2], pre->g_pre_comp[8][0],
2116 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
2117 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2118 pre->g_pre_comp[2][2]);
2119 /* 2^260*G + 2^390*G */
2120 point_add(pre->g_pre_comp[12][0], pre->g_pre_comp[12][1],
2121 pre->g_pre_comp[12][2], pre->g_pre_comp[8][0],
2122 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
2123 0, pre->g_pre_comp[4][0], pre->g_pre_comp[4][1],
2124 pre->g_pre_comp[4][2]);
2125 /* 2^130*G + 2^260*G + 2^390*G */
2126 point_add(pre->g_pre_comp[14][0], pre->g_pre_comp[14][1],
2127 pre->g_pre_comp[14][2], pre->g_pre_comp[12][0],
2128 pre->g_pre_comp[12][1], pre->g_pre_comp[12][2],
2129 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2130 pre->g_pre_comp[2][2]);
2131 for (i = 1; i < 8; ++i) {
2132 /* odd multiples: add G */
2133 point_add(pre->g_pre_comp[2 * i + 1][0],
2134 pre->g_pre_comp[2 * i + 1][1],
2135 pre->g_pre_comp[2 * i + 1][2], pre->g_pre_comp[2 * i][0],
2136 pre->g_pre_comp[2 * i][1], pre->g_pre_comp[2 * i][2], 0,
2137 pre->g_pre_comp[1][0], pre->g_pre_comp[1][1],
2138 pre->g_pre_comp[1][2]);
2140 make_points_affine(15, &(pre->g_pre_comp[1]), tmp_felems);
2143 SETPRECOMP(group, nistp521, pre);
2148 EC_POINT_free(generator);
2149 BN_CTX_free(new_ctx);
2150 EC_nistp521_pre_comp_free(pre);
2154 int ec_GFp_nistp521_have_precompute_mult(const EC_GROUP *group)
2156 return HAVEPRECOMP(group, nistp521);