2 * Copyright 2011-2018 The OpenSSL Project Authors. All Rights Reserved.
4 * Licensed under the Apache License 2.0 (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(__SIZEOF_INT128__) && __SIZEOF_INT128__==16
44 /* even with gcc, the typedef won't work for 32-bit platforms */
45 typedef __uint128_t uint128_t; /* nonstandard; implemented by gcc on 64-bit
48 # error "Your compiler doesn't appear to support 128-bit integer types"
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.
60 typedef u8 felem_bytearray[66];
63 * These are the parameters of P521, taken from FIPS 186-3, section D.1.2.5.
64 * These values are big-endian.
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,
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,
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,
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,
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,
115 * The representation of field elements.
116 * ------------------------------------
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.
125 * A field element with 64-bit limbs is an 'felem'. One with 128-bit limbs is a
130 typedef uint64_t limb;
131 typedef limb felem[NLIMBS];
132 typedef uint128_t largefelem[NLIMBS];
134 static const limb bottom57bits = 0x1ffffffffffffff;
135 static const limb bottom58bits = 0x3ffffffffffffff;
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.
141 static void bin66_to_felem(felem out, const u8 in[66])
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;
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.
158 static void felem_to_bin66(u8 out[66], const felem in)
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];
172 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */
173 static void flip_endian(u8 *out, const u8 *in, unsigned len)
176 for (i = 0; i < len; ++i)
177 out[i] = in[len - 1 - i];
180 /* BN_to_felem converts an OpenSSL BIGNUM into an felem */
181 static int BN_to_felem(felem out, const BIGNUM *bn)
183 felem_bytearray b_in;
184 felem_bytearray b_out;
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);
194 if (BN_is_negative(bn)) {
195 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
198 num_bytes = BN_bn2bin(bn, b_in);
199 flip_endian(b_out, b_in, num_bytes);
200 bin66_to_felem(out, b_out);
204 /* felem_to_BN converts an felem into an OpenSSL BIGNUM */
205 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
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);
218 static void felem_one(felem out)
231 static void felem_assign(felem out, const felem in)
244 /* felem_sum64 sets out = out + in. */
245 static void felem_sum64(felem out, const felem in)
258 /* felem_scalar sets out = in * scalar */
259 static void felem_scalar(felem out, const felem in, limb scalar)
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;
272 /* felem_scalar64 sets out = out * scalar */
273 static void felem_scalar64(felem out, limb scalar)
286 /* felem_scalar128 sets out = out * scalar */
287 static void felem_scalar128(largefelem out, limb scalar)
301 * felem_neg sets |out| to |-in|
303 * in[i] < 2^59 + 2^14
307 static void felem_neg(felem out, const felem in)
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);
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];
325 * felem_diff64 subtracts |in| from |out|
327 * in[i] < 2^59 + 2^14
329 * out[i] < out[i] + 2^62
331 static void felem_diff64(felem out, const felem in)
334 * In order to prevent underflow, we add 0 mod p before subtracting.
336 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5);
337 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4);
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];
351 * felem_diff_128_64 subtracts |in| from |out|
353 * in[i] < 2^62 + 2^17
355 * out[i] < out[i] + 2^63
357 static void felem_diff_128_64(largefelem out, const felem in)
360 * In order to prevent underflow, we add 64p mod p (which is equivalent
361 * to 0 mod p) before subtracting. p is 2^521 - 1, i.e. in binary a 521
362 * digit number with all bits set to 1. See "The representation of field
363 * elements" comment above for a description of how limbs are used to
364 * represent a number. 64p is represented with 8 limbs containing a number
365 * with 58 bits set and one limb with a number with 57 bits set.
367 static const limb two63m6 = (((limb) 1) << 63) - (((limb) 1) << 6);
368 static const limb two63m5 = (((limb) 1) << 63) - (((limb) 1) << 5);
370 out[0] += two63m6 - in[0];
371 out[1] += two63m5 - in[1];
372 out[2] += two63m5 - in[2];
373 out[3] += two63m5 - in[3];
374 out[4] += two63m5 - in[4];
375 out[5] += two63m5 - in[5];
376 out[6] += two63m5 - in[6];
377 out[7] += two63m5 - in[7];
378 out[8] += two63m5 - in[8];
382 * felem_diff_128_64 subtracts |in| from |out|
386 * out[i] < out[i] + 2^127 - 2^69
388 static void felem_diff128(largefelem out, const largefelem in)
391 * In order to prevent underflow, we add 0 mod p before subtracting.
393 static const uint128_t two127m70 =
394 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 70);
395 static const uint128_t two127m69 =
396 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 69);
398 out[0] += (two127m70 - in[0]);
399 out[1] += (two127m69 - in[1]);
400 out[2] += (two127m69 - in[2]);
401 out[3] += (two127m69 - in[3]);
402 out[4] += (two127m69 - in[4]);
403 out[5] += (two127m69 - in[5]);
404 out[6] += (two127m69 - in[6]);
405 out[7] += (two127m69 - in[7]);
406 out[8] += (two127m69 - in[8]);
410 * felem_square sets |out| = |in|^2
414 * out[i] < 17 * max(in[i]) * max(in[i])
416 static void felem_square(largefelem out, const felem in)
419 felem_scalar(inx2, in, 2);
420 felem_scalar(inx4, in, 4);
423 * We have many cases were we want to do
426 * This is obviously just
428 * However, rather than do the doubling on the 128 bit result, we
429 * double one of the inputs to the multiplication by reading from
433 out[0] = ((uint128_t) in[0]) * in[0];
434 out[1] = ((uint128_t) in[0]) * inx2[1];
435 out[2] = ((uint128_t) in[0]) * inx2[2] + ((uint128_t) in[1]) * in[1];
436 out[3] = ((uint128_t) in[0]) * inx2[3] + ((uint128_t) in[1]) * inx2[2];
437 out[4] = ((uint128_t) in[0]) * inx2[4] +
438 ((uint128_t) in[1]) * inx2[3] + ((uint128_t) in[2]) * in[2];
439 out[5] = ((uint128_t) in[0]) * inx2[5] +
440 ((uint128_t) in[1]) * inx2[4] + ((uint128_t) in[2]) * inx2[3];
441 out[6] = ((uint128_t) in[0]) * inx2[6] +
442 ((uint128_t) in[1]) * inx2[5] +
443 ((uint128_t) in[2]) * inx2[4] + ((uint128_t) in[3]) * in[3];
444 out[7] = ((uint128_t) in[0]) * inx2[7] +
445 ((uint128_t) in[1]) * inx2[6] +
446 ((uint128_t) in[2]) * inx2[5] + ((uint128_t) in[3]) * inx2[4];
447 out[8] = ((uint128_t) in[0]) * inx2[8] +
448 ((uint128_t) in[1]) * inx2[7] +
449 ((uint128_t) in[2]) * inx2[6] +
450 ((uint128_t) in[3]) * inx2[5] + ((uint128_t) in[4]) * in[4];
453 * The remaining limbs fall above 2^521, with the first falling at 2^522.
454 * They correspond to locations one bit up from the limbs produced above
455 * so we would have to multiply by two to align them. Again, rather than
456 * operate on the 128-bit result, we double one of the inputs to the
457 * multiplication. If we want to double for both this reason, and the
458 * reason above, then we end up multiplying by four.
462 out[0] += ((uint128_t) in[1]) * inx4[8] +
463 ((uint128_t) in[2]) * inx4[7] +
464 ((uint128_t) in[3]) * inx4[6] + ((uint128_t) in[4]) * inx4[5];
467 out[1] += ((uint128_t) in[2]) * inx4[8] +
468 ((uint128_t) in[3]) * inx4[7] +
469 ((uint128_t) in[4]) * inx4[6] + ((uint128_t) in[5]) * inx2[5];
472 out[2] += ((uint128_t) in[3]) * inx4[8] +
473 ((uint128_t) in[4]) * inx4[7] + ((uint128_t) in[5]) * inx4[6];
476 out[3] += ((uint128_t) in[4]) * inx4[8] +
477 ((uint128_t) in[5]) * inx4[7] + ((uint128_t) in[6]) * inx2[6];
480 out[4] += ((uint128_t) in[5]) * inx4[8] + ((uint128_t) in[6]) * inx4[7];
483 out[5] += ((uint128_t) in[6]) * inx4[8] + ((uint128_t) in[7]) * inx2[7];
486 out[6] += ((uint128_t) in[7]) * inx4[8];
489 out[7] += ((uint128_t) in[8]) * inx2[8];
493 * felem_mul sets |out| = |in1| * |in2|
498 * out[i] < 17 * max(in1[i]) * max(in2[i])
500 static void felem_mul(largefelem out, const felem in1, const felem in2)
503 felem_scalar(in2x2, in2, 2);
505 out[0] = ((uint128_t) in1[0]) * in2[0];
507 out[1] = ((uint128_t) in1[0]) * in2[1] +
508 ((uint128_t) in1[1]) * in2[0];
510 out[2] = ((uint128_t) in1[0]) * in2[2] +
511 ((uint128_t) in1[1]) * in2[1] +
512 ((uint128_t) in1[2]) * in2[0];
514 out[3] = ((uint128_t) in1[0]) * in2[3] +
515 ((uint128_t) in1[1]) * in2[2] +
516 ((uint128_t) in1[2]) * in2[1] +
517 ((uint128_t) in1[3]) * in2[0];
519 out[4] = ((uint128_t) in1[0]) * in2[4] +
520 ((uint128_t) in1[1]) * in2[3] +
521 ((uint128_t) in1[2]) * in2[2] +
522 ((uint128_t) in1[3]) * in2[1] +
523 ((uint128_t) in1[4]) * in2[0];
525 out[5] = ((uint128_t) in1[0]) * in2[5] +
526 ((uint128_t) in1[1]) * in2[4] +
527 ((uint128_t) in1[2]) * in2[3] +
528 ((uint128_t) in1[3]) * in2[2] +
529 ((uint128_t) in1[4]) * in2[1] +
530 ((uint128_t) in1[5]) * in2[0];
532 out[6] = ((uint128_t) in1[0]) * in2[6] +
533 ((uint128_t) in1[1]) * in2[5] +
534 ((uint128_t) in1[2]) * in2[4] +
535 ((uint128_t) in1[3]) * in2[3] +
536 ((uint128_t) in1[4]) * in2[2] +
537 ((uint128_t) in1[5]) * in2[1] +
538 ((uint128_t) in1[6]) * in2[0];
540 out[7] = ((uint128_t) in1[0]) * in2[7] +
541 ((uint128_t) in1[1]) * in2[6] +
542 ((uint128_t) in1[2]) * in2[5] +
543 ((uint128_t) in1[3]) * in2[4] +
544 ((uint128_t) in1[4]) * in2[3] +
545 ((uint128_t) in1[5]) * in2[2] +
546 ((uint128_t) in1[6]) * in2[1] +
547 ((uint128_t) in1[7]) * in2[0];
549 out[8] = ((uint128_t) in1[0]) * in2[8] +
550 ((uint128_t) in1[1]) * in2[7] +
551 ((uint128_t) in1[2]) * in2[6] +
552 ((uint128_t) in1[3]) * in2[5] +
553 ((uint128_t) in1[4]) * in2[4] +
554 ((uint128_t) in1[5]) * in2[3] +
555 ((uint128_t) in1[6]) * in2[2] +
556 ((uint128_t) in1[7]) * in2[1] +
557 ((uint128_t) in1[8]) * in2[0];
559 /* See comment in felem_square about the use of in2x2 here */
561 out[0] += ((uint128_t) in1[1]) * in2x2[8] +
562 ((uint128_t) in1[2]) * in2x2[7] +
563 ((uint128_t) in1[3]) * in2x2[6] +
564 ((uint128_t) in1[4]) * in2x2[5] +
565 ((uint128_t) in1[5]) * in2x2[4] +
566 ((uint128_t) in1[6]) * in2x2[3] +
567 ((uint128_t) in1[7]) * in2x2[2] +
568 ((uint128_t) in1[8]) * in2x2[1];
570 out[1] += ((uint128_t) in1[2]) * in2x2[8] +
571 ((uint128_t) in1[3]) * in2x2[7] +
572 ((uint128_t) in1[4]) * in2x2[6] +
573 ((uint128_t) in1[5]) * in2x2[5] +
574 ((uint128_t) in1[6]) * in2x2[4] +
575 ((uint128_t) in1[7]) * in2x2[3] +
576 ((uint128_t) in1[8]) * in2x2[2];
578 out[2] += ((uint128_t) in1[3]) * in2x2[8] +
579 ((uint128_t) in1[4]) * in2x2[7] +
580 ((uint128_t) in1[5]) * in2x2[6] +
581 ((uint128_t) in1[6]) * in2x2[5] +
582 ((uint128_t) in1[7]) * in2x2[4] +
583 ((uint128_t) in1[8]) * in2x2[3];
585 out[3] += ((uint128_t) in1[4]) * in2x2[8] +
586 ((uint128_t) in1[5]) * in2x2[7] +
587 ((uint128_t) in1[6]) * in2x2[6] +
588 ((uint128_t) in1[7]) * in2x2[5] +
589 ((uint128_t) in1[8]) * in2x2[4];
591 out[4] += ((uint128_t) in1[5]) * in2x2[8] +
592 ((uint128_t) in1[6]) * in2x2[7] +
593 ((uint128_t) in1[7]) * in2x2[6] +
594 ((uint128_t) in1[8]) * in2x2[5];
596 out[5] += ((uint128_t) in1[6]) * in2x2[8] +
597 ((uint128_t) in1[7]) * in2x2[7] +
598 ((uint128_t) in1[8]) * in2x2[6];
600 out[6] += ((uint128_t) in1[7]) * in2x2[8] +
601 ((uint128_t) in1[8]) * in2x2[7];
603 out[7] += ((uint128_t) in1[8]) * in2x2[8];
606 static const limb bottom52bits = 0xfffffffffffff;
609 * felem_reduce converts a largefelem to an felem.
613 * out[i] < 2^59 + 2^14
615 static void felem_reduce(felem out, const largefelem in)
617 u64 overflow1, overflow2;
619 out[0] = ((limb) in[0]) & bottom58bits;
620 out[1] = ((limb) in[1]) & bottom58bits;
621 out[2] = ((limb) in[2]) & bottom58bits;
622 out[3] = ((limb) in[3]) & bottom58bits;
623 out[4] = ((limb) in[4]) & bottom58bits;
624 out[5] = ((limb) in[5]) & bottom58bits;
625 out[6] = ((limb) in[6]) & bottom58bits;
626 out[7] = ((limb) in[7]) & bottom58bits;
627 out[8] = ((limb) in[8]) & bottom58bits;
631 out[1] += ((limb) in[0]) >> 58;
632 out[1] += (((limb) (in[0] >> 64)) & bottom52bits) << 6;
634 * out[1] < 2^58 + 2^6 + 2^58
637 out[2] += ((limb) (in[0] >> 64)) >> 52;
639 out[2] += ((limb) in[1]) >> 58;
640 out[2] += (((limb) (in[1] >> 64)) & bottom52bits) << 6;
641 out[3] += ((limb) (in[1] >> 64)) >> 52;
643 out[3] += ((limb) in[2]) >> 58;
644 out[3] += (((limb) (in[2] >> 64)) & bottom52bits) << 6;
645 out[4] += ((limb) (in[2] >> 64)) >> 52;
647 out[4] += ((limb) in[3]) >> 58;
648 out[4] += (((limb) (in[3] >> 64)) & bottom52bits) << 6;
649 out[5] += ((limb) (in[3] >> 64)) >> 52;
651 out[5] += ((limb) in[4]) >> 58;
652 out[5] += (((limb) (in[4] >> 64)) & bottom52bits) << 6;
653 out[6] += ((limb) (in[4] >> 64)) >> 52;
655 out[6] += ((limb) in[5]) >> 58;
656 out[6] += (((limb) (in[5] >> 64)) & bottom52bits) << 6;
657 out[7] += ((limb) (in[5] >> 64)) >> 52;
659 out[7] += ((limb) in[6]) >> 58;
660 out[7] += (((limb) (in[6] >> 64)) & bottom52bits) << 6;
661 out[8] += ((limb) (in[6] >> 64)) >> 52;
663 out[8] += ((limb) in[7]) >> 58;
664 out[8] += (((limb) (in[7] >> 64)) & bottom52bits) << 6;
666 * out[x > 1] < 2^58 + 2^6 + 2^58 + 2^12
669 overflow1 = ((limb) (in[7] >> 64)) >> 52;
671 overflow1 += ((limb) in[8]) >> 58;
672 overflow1 += (((limb) (in[8] >> 64)) & bottom52bits) << 6;
673 overflow2 = ((limb) (in[8] >> 64)) >> 52;
675 overflow1 <<= 1; /* overflow1 < 2^13 + 2^7 + 2^59 */
676 overflow2 <<= 1; /* overflow2 < 2^13 */
678 out[0] += overflow1; /* out[0] < 2^60 */
679 out[1] += overflow2; /* out[1] < 2^59 + 2^6 + 2^13 */
681 out[1] += out[0] >> 58;
682 out[0] &= bottom58bits;
685 * out[1] < 2^59 + 2^6 + 2^13 + 2^2
690 static void felem_square_reduce(felem out, const felem in)
693 felem_square(tmp, in);
694 felem_reduce(out, tmp);
697 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
700 felem_mul(tmp, in1, in2);
701 felem_reduce(out, tmp);
705 * felem_inv calculates |out| = |in|^{-1}
707 * Based on Fermat's Little Theorem:
709 * a^{p-1} = 1 (mod p)
710 * a^{p-2} = a^{-1} (mod p)
712 static void felem_inv(felem out, const felem in)
714 felem ftmp, ftmp2, ftmp3, ftmp4;
718 felem_square(tmp, in);
719 felem_reduce(ftmp, tmp); /* 2^1 */
720 felem_mul(tmp, in, ftmp);
721 felem_reduce(ftmp, tmp); /* 2^2 - 2^0 */
722 felem_assign(ftmp2, ftmp);
723 felem_square(tmp, ftmp);
724 felem_reduce(ftmp, tmp); /* 2^3 - 2^1 */
725 felem_mul(tmp, in, ftmp);
726 felem_reduce(ftmp, tmp); /* 2^3 - 2^0 */
727 felem_square(tmp, ftmp);
728 felem_reduce(ftmp, tmp); /* 2^4 - 2^1 */
730 felem_square(tmp, ftmp2);
731 felem_reduce(ftmp3, tmp); /* 2^3 - 2^1 */
732 felem_square(tmp, ftmp3);
733 felem_reduce(ftmp3, tmp); /* 2^4 - 2^2 */
734 felem_mul(tmp, ftmp3, ftmp2);
735 felem_reduce(ftmp3, tmp); /* 2^4 - 2^0 */
737 felem_assign(ftmp2, ftmp3);
738 felem_square(tmp, ftmp3);
739 felem_reduce(ftmp3, tmp); /* 2^5 - 2^1 */
740 felem_square(tmp, ftmp3);
741 felem_reduce(ftmp3, tmp); /* 2^6 - 2^2 */
742 felem_square(tmp, ftmp3);
743 felem_reduce(ftmp3, tmp); /* 2^7 - 2^3 */
744 felem_square(tmp, ftmp3);
745 felem_reduce(ftmp3, tmp); /* 2^8 - 2^4 */
746 felem_assign(ftmp4, ftmp3);
747 felem_mul(tmp, ftmp3, ftmp);
748 felem_reduce(ftmp4, tmp); /* 2^8 - 2^1 */
749 felem_square(tmp, ftmp4);
750 felem_reduce(ftmp4, tmp); /* 2^9 - 2^2 */
751 felem_mul(tmp, ftmp3, ftmp2);
752 felem_reduce(ftmp3, tmp); /* 2^8 - 2^0 */
753 felem_assign(ftmp2, ftmp3);
755 for (i = 0; i < 8; i++) {
756 felem_square(tmp, ftmp3);
757 felem_reduce(ftmp3, tmp); /* 2^16 - 2^8 */
759 felem_mul(tmp, ftmp3, ftmp2);
760 felem_reduce(ftmp3, tmp); /* 2^16 - 2^0 */
761 felem_assign(ftmp2, ftmp3);
763 for (i = 0; i < 16; i++) {
764 felem_square(tmp, ftmp3);
765 felem_reduce(ftmp3, tmp); /* 2^32 - 2^16 */
767 felem_mul(tmp, ftmp3, ftmp2);
768 felem_reduce(ftmp3, tmp); /* 2^32 - 2^0 */
769 felem_assign(ftmp2, ftmp3);
771 for (i = 0; i < 32; i++) {
772 felem_square(tmp, ftmp3);
773 felem_reduce(ftmp3, tmp); /* 2^64 - 2^32 */
775 felem_mul(tmp, ftmp3, ftmp2);
776 felem_reduce(ftmp3, tmp); /* 2^64 - 2^0 */
777 felem_assign(ftmp2, ftmp3);
779 for (i = 0; i < 64; i++) {
780 felem_square(tmp, ftmp3);
781 felem_reduce(ftmp3, tmp); /* 2^128 - 2^64 */
783 felem_mul(tmp, ftmp3, ftmp2);
784 felem_reduce(ftmp3, tmp); /* 2^128 - 2^0 */
785 felem_assign(ftmp2, ftmp3);
787 for (i = 0; i < 128; i++) {
788 felem_square(tmp, ftmp3);
789 felem_reduce(ftmp3, tmp); /* 2^256 - 2^128 */
791 felem_mul(tmp, ftmp3, ftmp2);
792 felem_reduce(ftmp3, tmp); /* 2^256 - 2^0 */
793 felem_assign(ftmp2, ftmp3);
795 for (i = 0; i < 256; i++) {
796 felem_square(tmp, ftmp3);
797 felem_reduce(ftmp3, tmp); /* 2^512 - 2^256 */
799 felem_mul(tmp, ftmp3, ftmp2);
800 felem_reduce(ftmp3, tmp); /* 2^512 - 2^0 */
802 for (i = 0; i < 9; i++) {
803 felem_square(tmp, ftmp3);
804 felem_reduce(ftmp3, tmp); /* 2^521 - 2^9 */
806 felem_mul(tmp, ftmp3, ftmp4);
807 felem_reduce(ftmp3, tmp); /* 2^512 - 2^2 */
808 felem_mul(tmp, ftmp3, in);
809 felem_reduce(out, tmp); /* 2^512 - 3 */
812 /* This is 2^521-1, expressed as an felem */
813 static const felem kPrime = {
814 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff,
815 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff,
816 0x03ffffffffffffff, 0x03ffffffffffffff, 0x01ffffffffffffff
820 * felem_is_zero returns a limb with all bits set if |in| == 0 (mod p) and 0
823 * in[i] < 2^59 + 2^14
825 static limb felem_is_zero(const felem in)
829 felem_assign(ftmp, in);
831 ftmp[0] += ftmp[8] >> 57;
832 ftmp[8] &= bottom57bits;
834 ftmp[1] += ftmp[0] >> 58;
835 ftmp[0] &= bottom58bits;
836 ftmp[2] += ftmp[1] >> 58;
837 ftmp[1] &= bottom58bits;
838 ftmp[3] += ftmp[2] >> 58;
839 ftmp[2] &= bottom58bits;
840 ftmp[4] += ftmp[3] >> 58;
841 ftmp[3] &= bottom58bits;
842 ftmp[5] += ftmp[4] >> 58;
843 ftmp[4] &= bottom58bits;
844 ftmp[6] += ftmp[5] >> 58;
845 ftmp[5] &= bottom58bits;
846 ftmp[7] += ftmp[6] >> 58;
847 ftmp[6] &= bottom58bits;
848 ftmp[8] += ftmp[7] >> 58;
849 ftmp[7] &= bottom58bits;
850 /* ftmp[8] < 2^57 + 4 */
853 * The ninth limb of 2*(2^521-1) is 0x03ffffffffffffff, which is greater
854 * than our bound for ftmp[8]. Therefore we only have to check if the
855 * zero is zero or 2^521-1.
871 * We know that ftmp[i] < 2^63, therefore the only way that the top bit
872 * can be set is if is_zero was 0 before the decrement.
874 is_zero = 0 - (is_zero >> 63);
876 is_p = ftmp[0] ^ kPrime[0];
877 is_p |= ftmp[1] ^ kPrime[1];
878 is_p |= ftmp[2] ^ kPrime[2];
879 is_p |= ftmp[3] ^ kPrime[3];
880 is_p |= ftmp[4] ^ kPrime[4];
881 is_p |= ftmp[5] ^ kPrime[5];
882 is_p |= ftmp[6] ^ kPrime[6];
883 is_p |= ftmp[7] ^ kPrime[7];
884 is_p |= ftmp[8] ^ kPrime[8];
887 is_p = 0 - (is_p >> 63);
893 static int felem_is_zero_int(const void *in)
895 return (int)(felem_is_zero(in) & ((limb) 1));
899 * felem_contract converts |in| to its unique, minimal representation.
901 * in[i] < 2^59 + 2^14
903 static void felem_contract(felem out, const felem in)
905 limb is_p, is_greater, sign;
906 static const limb two58 = ((limb) 1) << 58;
908 felem_assign(out, in);
910 out[0] += out[8] >> 57;
911 out[8] &= bottom57bits;
913 out[1] += out[0] >> 58;
914 out[0] &= bottom58bits;
915 out[2] += out[1] >> 58;
916 out[1] &= bottom58bits;
917 out[3] += out[2] >> 58;
918 out[2] &= bottom58bits;
919 out[4] += out[3] >> 58;
920 out[3] &= bottom58bits;
921 out[5] += out[4] >> 58;
922 out[4] &= bottom58bits;
923 out[6] += out[5] >> 58;
924 out[5] &= bottom58bits;
925 out[7] += out[6] >> 58;
926 out[6] &= bottom58bits;
927 out[8] += out[7] >> 58;
928 out[7] &= bottom58bits;
929 /* out[8] < 2^57 + 4 */
932 * If the value is greater than 2^521-1 then we have to subtract 2^521-1
933 * out. See the comments in felem_is_zero regarding why we don't test for
934 * other multiples of the prime.
938 * First, if |out| is equal to 2^521-1, we subtract it out to get zero.
941 is_p = out[0] ^ kPrime[0];
942 is_p |= out[1] ^ kPrime[1];
943 is_p |= out[2] ^ kPrime[2];
944 is_p |= out[3] ^ kPrime[3];
945 is_p |= out[4] ^ kPrime[4];
946 is_p |= out[5] ^ kPrime[5];
947 is_p |= out[6] ^ kPrime[6];
948 is_p |= out[7] ^ kPrime[7];
949 is_p |= out[8] ^ kPrime[8];
958 is_p = 0 - (is_p >> 63);
961 /* is_p is 0 iff |out| == 2^521-1 and all ones otherwise */
974 * In order to test that |out| >= 2^521-1 we need only test if out[8] >>
975 * 57 is greater than zero as (2^521-1) + x >= 2^522
977 is_greater = out[8] >> 57;
978 is_greater |= is_greater << 32;
979 is_greater |= is_greater << 16;
980 is_greater |= is_greater << 8;
981 is_greater |= is_greater << 4;
982 is_greater |= is_greater << 2;
983 is_greater |= is_greater << 1;
984 is_greater = 0 - (is_greater >> 63);
986 out[0] -= kPrime[0] & is_greater;
987 out[1] -= kPrime[1] & is_greater;
988 out[2] -= kPrime[2] & is_greater;
989 out[3] -= kPrime[3] & is_greater;
990 out[4] -= kPrime[4] & is_greater;
991 out[5] -= kPrime[5] & is_greater;
992 out[6] -= kPrime[6] & is_greater;
993 out[7] -= kPrime[7] & is_greater;
994 out[8] -= kPrime[8] & is_greater;
996 /* Eliminate negative coefficients */
997 sign = -(out[0] >> 63);
998 out[0] += (two58 & sign);
999 out[1] -= (1 & sign);
1000 sign = -(out[1] >> 63);
1001 out[1] += (two58 & sign);
1002 out[2] -= (1 & sign);
1003 sign = -(out[2] >> 63);
1004 out[2] += (two58 & sign);
1005 out[3] -= (1 & sign);
1006 sign = -(out[3] >> 63);
1007 out[3] += (two58 & sign);
1008 out[4] -= (1 & sign);
1009 sign = -(out[4] >> 63);
1010 out[4] += (two58 & sign);
1011 out[5] -= (1 & sign);
1012 sign = -(out[0] >> 63);
1013 out[5] += (two58 & sign);
1014 out[6] -= (1 & sign);
1015 sign = -(out[6] >> 63);
1016 out[6] += (two58 & sign);
1017 out[7] -= (1 & sign);
1018 sign = -(out[7] >> 63);
1019 out[7] += (two58 & sign);
1020 out[8] -= (1 & sign);
1021 sign = -(out[5] >> 63);
1022 out[5] += (two58 & sign);
1023 out[6] -= (1 & sign);
1024 sign = -(out[6] >> 63);
1025 out[6] += (two58 & sign);
1026 out[7] -= (1 & sign);
1027 sign = -(out[7] >> 63);
1028 out[7] += (two58 & sign);
1029 out[8] -= (1 & sign);
1036 * Building on top of the field operations we have the operations on the
1037 * elliptic curve group itself. Points on the curve are represented in Jacobian
1041 * point_double calculates 2*(x_in, y_in, z_in)
1043 * The method is taken from:
1044 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
1046 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed.
1047 * while x_out == y_in is not (maybe this works, but it's not tested). */
1049 point_double(felem x_out, felem y_out, felem z_out,
1050 const felem x_in, const felem y_in, const felem z_in)
1052 largefelem tmp, tmp2;
1053 felem delta, gamma, beta, alpha, ftmp, ftmp2;
1055 felem_assign(ftmp, x_in);
1056 felem_assign(ftmp2, x_in);
1059 felem_square(tmp, z_in);
1060 felem_reduce(delta, tmp); /* delta[i] < 2^59 + 2^14 */
1063 felem_square(tmp, y_in);
1064 felem_reduce(gamma, tmp); /* gamma[i] < 2^59 + 2^14 */
1066 /* beta = x*gamma */
1067 felem_mul(tmp, x_in, gamma);
1068 felem_reduce(beta, tmp); /* beta[i] < 2^59 + 2^14 */
1070 /* alpha = 3*(x-delta)*(x+delta) */
1071 felem_diff64(ftmp, delta);
1072 /* ftmp[i] < 2^61 */
1073 felem_sum64(ftmp2, delta);
1074 /* ftmp2[i] < 2^60 + 2^15 */
1075 felem_scalar64(ftmp2, 3);
1076 /* ftmp2[i] < 3*2^60 + 3*2^15 */
1077 felem_mul(tmp, ftmp, ftmp2);
1079 * tmp[i] < 17(3*2^121 + 3*2^76)
1080 * = 61*2^121 + 61*2^76
1081 * < 64*2^121 + 64*2^76
1085 felem_reduce(alpha, tmp);
1087 /* x' = alpha^2 - 8*beta */
1088 felem_square(tmp, alpha);
1090 * tmp[i] < 17*2^120 < 2^125
1092 felem_assign(ftmp, beta);
1093 felem_scalar64(ftmp, 8);
1094 /* ftmp[i] < 2^62 + 2^17 */
1095 felem_diff_128_64(tmp, ftmp);
1096 /* tmp[i] < 2^125 + 2^63 + 2^62 + 2^17 */
1097 felem_reduce(x_out, tmp);
1099 /* z' = (y + z)^2 - gamma - delta */
1100 felem_sum64(delta, gamma);
1101 /* delta[i] < 2^60 + 2^15 */
1102 felem_assign(ftmp, y_in);
1103 felem_sum64(ftmp, z_in);
1104 /* ftmp[i] < 2^60 + 2^15 */
1105 felem_square(tmp, ftmp);
1107 * tmp[i] < 17(2^122) < 2^127
1109 felem_diff_128_64(tmp, delta);
1110 /* tmp[i] < 2^127 + 2^63 */
1111 felem_reduce(z_out, tmp);
1113 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
1114 felem_scalar64(beta, 4);
1115 /* beta[i] < 2^61 + 2^16 */
1116 felem_diff64(beta, x_out);
1117 /* beta[i] < 2^61 + 2^60 + 2^16 */
1118 felem_mul(tmp, alpha, beta);
1120 * tmp[i] < 17*((2^59 + 2^14)(2^61 + 2^60 + 2^16))
1121 * = 17*(2^120 + 2^75 + 2^119 + 2^74 + 2^75 + 2^30)
1122 * = 17*(2^120 + 2^119 + 2^76 + 2^74 + 2^30)
1125 felem_square(tmp2, gamma);
1127 * tmp2[i] < 17*(2^59 + 2^14)^2
1128 * = 17*(2^118 + 2^74 + 2^28)
1130 felem_scalar128(tmp2, 8);
1132 * tmp2[i] < 8*17*(2^118 + 2^74 + 2^28)
1133 * = 2^125 + 2^121 + 2^81 + 2^77 + 2^35 + 2^31
1136 felem_diff128(tmp, tmp2);
1138 * tmp[i] < 2^127 - 2^69 + 17(2^120 + 2^119 + 2^76 + 2^74 + 2^30)
1139 * = 2^127 + 2^124 + 2^122 + 2^120 + 2^118 + 2^80 + 2^78 + 2^76 +
1140 * 2^74 + 2^69 + 2^34 + 2^30
1143 felem_reduce(y_out, tmp);
1146 /* copy_conditional copies in to out iff mask is all ones. */
1147 static void copy_conditional(felem out, const felem in, limb mask)
1150 for (i = 0; i < NLIMBS; ++i) {
1151 const limb tmp = mask & (in[i] ^ out[i]);
1157 * point_add calculates (x1, y1, z1) + (x2, y2, z2)
1159 * The method is taken from
1160 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl,
1161 * adapted for mixed addition (z2 = 1, or z2 = 0 for the point at infinity).
1163 * This function includes a branch for checking whether the two input points
1164 * are equal (while not equal to the point at infinity). See comment below
1167 static void point_add(felem x3, felem y3, felem z3,
1168 const felem x1, const felem y1, const felem z1,
1169 const int mixed, const felem x2, const felem y2,
1172 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, ftmp6, x_out, y_out, z_out;
1173 largefelem tmp, tmp2;
1174 limb x_equal, y_equal, z1_is_zero, z2_is_zero;
1176 z1_is_zero = felem_is_zero(z1);
1177 z2_is_zero = felem_is_zero(z2);
1179 /* ftmp = z1z1 = z1**2 */
1180 felem_square(tmp, z1);
1181 felem_reduce(ftmp, tmp);
1184 /* ftmp2 = z2z2 = z2**2 */
1185 felem_square(tmp, z2);
1186 felem_reduce(ftmp2, tmp);
1188 /* u1 = ftmp3 = x1*z2z2 */
1189 felem_mul(tmp, x1, ftmp2);
1190 felem_reduce(ftmp3, tmp);
1192 /* ftmp5 = z1 + z2 */
1193 felem_assign(ftmp5, z1);
1194 felem_sum64(ftmp5, z2);
1195 /* ftmp5[i] < 2^61 */
1197 /* ftmp5 = (z1 + z2)**2 - z1z1 - z2z2 = 2*z1z2 */
1198 felem_square(tmp, ftmp5);
1199 /* tmp[i] < 17*2^122 */
1200 felem_diff_128_64(tmp, ftmp);
1201 /* tmp[i] < 17*2^122 + 2^63 */
1202 felem_diff_128_64(tmp, ftmp2);
1203 /* tmp[i] < 17*2^122 + 2^64 */
1204 felem_reduce(ftmp5, tmp);
1206 /* ftmp2 = z2 * z2z2 */
1207 felem_mul(tmp, ftmp2, z2);
1208 felem_reduce(ftmp2, tmp);
1210 /* s1 = ftmp6 = y1 * z2**3 */
1211 felem_mul(tmp, y1, ftmp2);
1212 felem_reduce(ftmp6, tmp);
1215 * We'll assume z2 = 1 (special case z2 = 0 is handled later)
1218 /* u1 = ftmp3 = x1*z2z2 */
1219 felem_assign(ftmp3, x1);
1221 /* ftmp5 = 2*z1z2 */
1222 felem_scalar(ftmp5, z1, 2);
1224 /* s1 = ftmp6 = y1 * z2**3 */
1225 felem_assign(ftmp6, y1);
1229 felem_mul(tmp, x2, ftmp);
1230 /* tmp[i] < 17*2^120 */
1232 /* h = ftmp4 = u2 - u1 */
1233 felem_diff_128_64(tmp, ftmp3);
1234 /* tmp[i] < 17*2^120 + 2^63 */
1235 felem_reduce(ftmp4, tmp);
1237 x_equal = felem_is_zero(ftmp4);
1239 /* z_out = ftmp5 * h */
1240 felem_mul(tmp, ftmp5, ftmp4);
1241 felem_reduce(z_out, tmp);
1243 /* ftmp = z1 * z1z1 */
1244 felem_mul(tmp, ftmp, z1);
1245 felem_reduce(ftmp, tmp);
1247 /* s2 = tmp = y2 * z1**3 */
1248 felem_mul(tmp, y2, ftmp);
1249 /* tmp[i] < 17*2^120 */
1251 /* r = ftmp5 = (s2 - s1)*2 */
1252 felem_diff_128_64(tmp, ftmp6);
1253 /* tmp[i] < 17*2^120 + 2^63 */
1254 felem_reduce(ftmp5, tmp);
1255 y_equal = felem_is_zero(ftmp5);
1256 felem_scalar64(ftmp5, 2);
1257 /* ftmp5[i] < 2^61 */
1259 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
1261 * This is obviously not constant-time but it will almost-never happen
1262 * for ECDH / ECDSA. The case where it can happen is during scalar-mult
1263 * where the intermediate value gets very close to the group order.
1264 * Since |ec_GFp_nistp_recode_scalar_bits| produces signed digits for
1265 * the scalar, it's possible for the intermediate value to be a small
1266 * negative multiple of the base point, and for the final signed digit
1267 * to be the same value. We believe that this only occurs for the scalar
1268 * 1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
1269 * ffffffa51868783bf2f966b7fcc0148f709a5d03bb5c9b8899c47aebb6fb
1270 * 71e913863f7, in that case the penultimate intermediate is -9G and
1271 * the final digit is also -9G. Since this only happens for a single
1272 * scalar, the timing leak is irrelevant. (Any attacker who wanted to
1273 * check whether a secret scalar was that exact value, can already do
1276 point_double(x3, y3, z3, x1, y1, z1);
1280 /* I = ftmp = (2h)**2 */
1281 felem_assign(ftmp, ftmp4);
1282 felem_scalar64(ftmp, 2);
1283 /* ftmp[i] < 2^61 */
1284 felem_square(tmp, ftmp);
1285 /* tmp[i] < 17*2^122 */
1286 felem_reduce(ftmp, tmp);
1288 /* J = ftmp2 = h * I */
1289 felem_mul(tmp, ftmp4, ftmp);
1290 felem_reduce(ftmp2, tmp);
1292 /* V = ftmp4 = U1 * I */
1293 felem_mul(tmp, ftmp3, ftmp);
1294 felem_reduce(ftmp4, tmp);
1296 /* x_out = r**2 - J - 2V */
1297 felem_square(tmp, ftmp5);
1298 /* tmp[i] < 17*2^122 */
1299 felem_diff_128_64(tmp, ftmp2);
1300 /* tmp[i] < 17*2^122 + 2^63 */
1301 felem_assign(ftmp3, ftmp4);
1302 felem_scalar64(ftmp4, 2);
1303 /* ftmp4[i] < 2^61 */
1304 felem_diff_128_64(tmp, ftmp4);
1305 /* tmp[i] < 17*2^122 + 2^64 */
1306 felem_reduce(x_out, tmp);
1308 /* y_out = r(V-x_out) - 2 * s1 * J */
1309 felem_diff64(ftmp3, x_out);
1311 * ftmp3[i] < 2^60 + 2^60 = 2^61
1313 felem_mul(tmp, ftmp5, ftmp3);
1314 /* tmp[i] < 17*2^122 */
1315 felem_mul(tmp2, ftmp6, ftmp2);
1316 /* tmp2[i] < 17*2^120 */
1317 felem_scalar128(tmp2, 2);
1318 /* tmp2[i] < 17*2^121 */
1319 felem_diff128(tmp, tmp2);
1321 * tmp[i] < 2^127 - 2^69 + 17*2^122
1322 * = 2^126 - 2^122 - 2^6 - 2^2 - 1
1325 felem_reduce(y_out, tmp);
1327 copy_conditional(x_out, x2, z1_is_zero);
1328 copy_conditional(x_out, x1, z2_is_zero);
1329 copy_conditional(y_out, y2, z1_is_zero);
1330 copy_conditional(y_out, y1, z2_is_zero);
1331 copy_conditional(z_out, z2, z1_is_zero);
1332 copy_conditional(z_out, z1, z2_is_zero);
1333 felem_assign(x3, x_out);
1334 felem_assign(y3, y_out);
1335 felem_assign(z3, z_out);
1339 * Base point pre computation
1340 * --------------------------
1342 * Two different sorts of precomputed tables are used in the following code.
1343 * Each contain various points on the curve, where each point is three field
1344 * elements (x, y, z).
1346 * For the base point table, z is usually 1 (0 for the point at infinity).
1347 * This table has 16 elements:
1348 * index | bits | point
1349 * ------+---------+------------------------------
1352 * 2 | 0 0 1 0 | 2^130G
1353 * 3 | 0 0 1 1 | (2^130 + 1)G
1354 * 4 | 0 1 0 0 | 2^260G
1355 * 5 | 0 1 0 1 | (2^260 + 1)G
1356 * 6 | 0 1 1 0 | (2^260 + 2^130)G
1357 * 7 | 0 1 1 1 | (2^260 + 2^130 + 1)G
1358 * 8 | 1 0 0 0 | 2^390G
1359 * 9 | 1 0 0 1 | (2^390 + 1)G
1360 * 10 | 1 0 1 0 | (2^390 + 2^130)G
1361 * 11 | 1 0 1 1 | (2^390 + 2^130 + 1)G
1362 * 12 | 1 1 0 0 | (2^390 + 2^260)G
1363 * 13 | 1 1 0 1 | (2^390 + 2^260 + 1)G
1364 * 14 | 1 1 1 0 | (2^390 + 2^260 + 2^130)G
1365 * 15 | 1 1 1 1 | (2^390 + 2^260 + 2^130 + 1)G
1367 * The reason for this is so that we can clock bits into four different
1368 * locations when doing simple scalar multiplies against the base point.
1370 * Tables for other points have table[i] = iG for i in 0 .. 16. */
1372 /* gmul is the table of precomputed base points */
1373 static const felem gmul[16][3] = {
1374 {{0, 0, 0, 0, 0, 0, 0, 0, 0},
1375 {0, 0, 0, 0, 0, 0, 0, 0, 0},
1376 {0, 0, 0, 0, 0, 0, 0, 0, 0}},
1377 {{0x017e7e31c2e5bd66, 0x022cf0615a90a6fe, 0x00127a2ffa8de334,
1378 0x01dfbf9d64a3f877, 0x006b4d3dbaa14b5e, 0x014fed487e0a2bd8,
1379 0x015b4429c6481390, 0x03a73678fb2d988e, 0x00c6858e06b70404},
1380 {0x00be94769fd16650, 0x031c21a89cb09022, 0x039013fad0761353,
1381 0x02657bd099031542, 0x03273e662c97ee72, 0x01e6d11a05ebef45,
1382 0x03d1bd998f544495, 0x03001172297ed0b1, 0x011839296a789a3b},
1383 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1384 {{0x0373faacbc875bae, 0x00f325023721c671, 0x00f666fd3dbde5ad,
1385 0x01a6932363f88ea7, 0x01fc6d9e13f9c47b, 0x03bcbffc2bbf734e,
1386 0x013ee3c3647f3a92, 0x029409fefe75d07d, 0x00ef9199963d85e5},
1387 {0x011173743ad5b178, 0x02499c7c21bf7d46, 0x035beaeabb8b1a58,
1388 0x00f989c4752ea0a3, 0x0101e1de48a9c1a3, 0x01a20076be28ba6c,
1389 0x02f8052e5eb2de95, 0x01bfe8f82dea117c, 0x0160074d3c36ddb7},
1390 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1391 {{0x012f3fc373393b3b, 0x03d3d6172f1419fa, 0x02adc943c0b86873,
1392 0x00d475584177952b, 0x012a4d1673750ee2, 0x00512517a0f13b0c,
1393 0x02b184671a7b1734, 0x0315b84236f1a50a, 0x00a4afc472edbdb9},
1394 {0x00152a7077f385c4, 0x03044007d8d1c2ee, 0x0065829d61d52b52,
1395 0x00494ff6b6631d0d, 0x00a11d94d5f06bcf, 0x02d2f89474d9282e,
1396 0x0241c5727c06eeb9, 0x0386928710fbdb9d, 0x01f883f727b0dfbe},
1397 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1398 {{0x019b0c3c9185544d, 0x006243a37c9d97db, 0x02ee3cbe030a2ad2,
1399 0x00cfdd946bb51e0d, 0x0271c00932606b91, 0x03f817d1ec68c561,
1400 0x03f37009806a369c, 0x03c1f30baf184fd5, 0x01091022d6d2f065},
1401 {0x0292c583514c45ed, 0x0316fca51f9a286c, 0x00300af507c1489a,
1402 0x0295f69008298cf1, 0x02c0ed8274943d7b, 0x016509b9b47a431e,
1403 0x02bc9de9634868ce, 0x005b34929bffcb09, 0x000c1a0121681524},
1404 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1405 {{0x0286abc0292fb9f2, 0x02665eee9805b3f7, 0x01ed7455f17f26d6,
1406 0x0346355b83175d13, 0x006284944cd0a097, 0x0191895bcdec5e51,
1407 0x02e288370afda7d9, 0x03b22312bfefa67a, 0x01d104d3fc0613fe},
1408 {0x0092421a12f7e47f, 0x0077a83fa373c501, 0x03bd25c5f696bd0d,
1409 0x035c41e4d5459761, 0x01ca0d1742b24f53, 0x00aaab27863a509c,
1410 0x018b6de47df73917, 0x025c0b771705cd01, 0x01fd51d566d760a7},
1411 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1412 {{0x01dd92ff6b0d1dbd, 0x039c5e2e8f8afa69, 0x0261ed13242c3b27,
1413 0x0382c6e67026e6a0, 0x01d60b10be2089f9, 0x03c15f3dce86723f,
1414 0x03c764a32d2a062d, 0x017307eac0fad056, 0x018207c0b96c5256},
1415 {0x0196a16d60e13154, 0x03e6ce74c0267030, 0x00ddbf2b4e52a5aa,
1416 0x012738241bbf31c8, 0x00ebe8dc04685a28, 0x024c2ad6d380d4a2,
1417 0x035ee062a6e62d0e, 0x0029ed74af7d3a0f, 0x00eef32aec142ebd},
1418 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1419 {{0x00c31ec398993b39, 0x03a9f45bcda68253, 0x00ac733c24c70890,
1420 0x00872b111401ff01, 0x01d178c23195eafb, 0x03bca2c816b87f74,
1421 0x0261a9af46fbad7a, 0x0324b2a8dd3d28f9, 0x00918121d8f24e23},
1422 {0x032bc8c1ca983cd7, 0x00d869dfb08fc8c6, 0x01693cb61fce1516,
1423 0x012a5ea68f4e88a8, 0x010869cab88d7ae3, 0x009081ad277ceee1,
1424 0x033a77166d064cdc, 0x03955235a1fb3a95, 0x01251a4a9b25b65e},
1425 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1426 {{0x00148a3a1b27f40b, 0x0123186df1b31fdc, 0x00026e7beaad34ce,
1427 0x01db446ac1d3dbba, 0x0299c1a33437eaec, 0x024540610183cbb7,
1428 0x0173bb0e9ce92e46, 0x02b937e43921214b, 0x01ab0436a9bf01b5},
1429 {0x0383381640d46948, 0x008dacbf0e7f330f, 0x03602122bcc3f318,
1430 0x01ee596b200620d6, 0x03bd0585fda430b3, 0x014aed77fd123a83,
1431 0x005ace749e52f742, 0x0390fe041da2b842, 0x0189a8ceb3299242},
1432 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1433 {{0x012a19d6b3282473, 0x00c0915918b423ce, 0x023a954eb94405ae,
1434 0x00529f692be26158, 0x0289fa1b6fa4b2aa, 0x0198ae4ceea346ef,
1435 0x0047d8cdfbdedd49, 0x00cc8c8953f0f6b8, 0x001424abbff49203},
1436 {0x0256732a1115a03a, 0x0351bc38665c6733, 0x03f7b950fb4a6447,
1437 0x000afffa94c22155, 0x025763d0a4dab540, 0x000511e92d4fc283,
1438 0x030a7e9eda0ee96c, 0x004c3cd93a28bf0a, 0x017edb3a8719217f},
1439 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1440 {{0x011de5675a88e673, 0x031d7d0f5e567fbe, 0x0016b2062c970ae5,
1441 0x03f4a2be49d90aa7, 0x03cef0bd13822866, 0x03f0923dcf774a6c,
1442 0x0284bebc4f322f72, 0x016ab2645302bb2c, 0x01793f95dace0e2a},
1443 {0x010646e13527a28f, 0x01ca1babd59dc5e7, 0x01afedfd9a5595df,
1444 0x01f15785212ea6b1, 0x0324e5d64f6ae3f4, 0x02d680f526d00645,
1445 0x0127920fadf627a7, 0x03b383f75df4f684, 0x0089e0057e783b0a},
1446 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1447 {{0x00f334b9eb3c26c6, 0x0298fdaa98568dce, 0x01c2d24843a82292,
1448 0x020bcb24fa1b0711, 0x02cbdb3d2b1875e6, 0x0014907598f89422,
1449 0x03abe3aa43b26664, 0x02cbf47f720bc168, 0x0133b5e73014b79b},
1450 {0x034aab5dab05779d, 0x00cdc5d71fee9abb, 0x0399f16bd4bd9d30,
1451 0x03582fa592d82647, 0x02be1cdfb775b0e9, 0x0034f7cea32e94cb,
1452 0x0335a7f08f56f286, 0x03b707e9565d1c8b, 0x0015c946ea5b614f},
1453 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1454 {{0x024676f6cff72255, 0x00d14625cac96378, 0x00532b6008bc3767,
1455 0x01fc16721b985322, 0x023355ea1b091668, 0x029de7afdc0317c3,
1456 0x02fc8a7ca2da037c, 0x02de1217d74a6f30, 0x013f7173175b73bf},
1457 {0x0344913f441490b5, 0x0200f9e272b61eca, 0x0258a246b1dd55d2,
1458 0x03753db9ea496f36, 0x025e02937a09c5ef, 0x030cbd3d14012692,
1459 0x01793a67e70dc72a, 0x03ec1d37048a662e, 0x006550f700c32a8d},
1460 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1461 {{0x00d3f48a347eba27, 0x008e636649b61bd8, 0x00d3b93716778fb3,
1462 0x004d1915757bd209, 0x019d5311a3da44e0, 0x016d1afcbbe6aade,
1463 0x0241bf5f73265616, 0x0384672e5d50d39b, 0x005009fee522b684},
1464 {0x029b4fab064435fe, 0x018868ee095bbb07, 0x01ea3d6936cc92b8,
1465 0x000608b00f78a2f3, 0x02db911073d1c20f, 0x018205938470100a,
1466 0x01f1e4964cbe6ff2, 0x021a19a29eed4663, 0x01414485f42afa81},
1467 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1468 {{0x01612b3a17f63e34, 0x03813992885428e6, 0x022b3c215b5a9608,
1469 0x029b4057e19f2fcb, 0x0384059a587af7e6, 0x02d6400ace6fe610,
1470 0x029354d896e8e331, 0x00c047ee6dfba65e, 0x0037720542e9d49d},
1471 {0x02ce9eed7c5e9278, 0x0374ed703e79643b, 0x01316c54c4072006,
1472 0x005aaa09054b2ee8, 0x002824000c840d57, 0x03d4eba24771ed86,
1473 0x0189c50aabc3bdae, 0x0338c01541e15510, 0x00466d56e38eed42},
1474 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1475 {{0x007efd8330ad8bd6, 0x02465ed48047710b, 0x0034c6606b215e0c,
1476 0x016ae30c53cbf839, 0x01fa17bd37161216, 0x018ead4e61ce8ab9,
1477 0x005482ed5f5dee46, 0x037543755bba1d7f, 0x005e5ac7e70a9d0f},
1478 {0x0117e1bb2fdcb2a2, 0x03deea36249f40c4, 0x028d09b4a6246cb7,
1479 0x03524b8855bcf756, 0x023d7d109d5ceb58, 0x0178e43e3223ef9c,
1480 0x0154536a0c6e966a, 0x037964d1286ee9fe, 0x0199bcd90e125055},
1481 {1, 0, 0, 0, 0, 0, 0, 0, 0}}
1485 * select_point selects the |idx|th point from a precomputation table and
1488 /* pre_comp below is of the size provided in |size| */
1489 static void select_point(const limb idx, unsigned int size,
1490 const felem pre_comp[][3], felem out[3])
1493 limb *outlimbs = &out[0][0];
1495 memset(out, 0, sizeof(*out) * 3);
1497 for (i = 0; i < size; i++) {
1498 const limb *inlimbs = &pre_comp[i][0][0];
1499 limb mask = i ^ idx;
1505 for (j = 0; j < NLIMBS * 3; j++)
1506 outlimbs[j] |= inlimbs[j] & mask;
1510 /* get_bit returns the |i|th bit in |in| */
1511 static char get_bit(const felem_bytearray in, int i)
1515 return (in[i >> 3] >> (i & 7)) & 1;
1519 * Interleaved point multiplication using precomputed point multiples: The
1520 * small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[], the scalars
1521 * in scalars[]. If g_scalar is non-NULL, we also add this multiple of the
1522 * generator, using certain (large) precomputed multiples in g_pre_comp.
1523 * Output point (X, Y, Z) is stored in x_out, y_out, z_out
1525 static void batch_mul(felem x_out, felem y_out, felem z_out,
1526 const felem_bytearray scalars[],
1527 const unsigned num_points, const u8 *g_scalar,
1528 const int mixed, const felem pre_comp[][17][3],
1529 const felem g_pre_comp[16][3])
1532 unsigned num, gen_mul = (g_scalar != NULL);
1533 felem nq[3], tmp[4];
1537 /* set nq to the point at infinity */
1538 memset(nq, 0, sizeof(nq));
1541 * Loop over all scalars msb-to-lsb, interleaving additions of multiples
1542 * of the generator (last quarter of rounds) and additions of other
1543 * points multiples (every 5th round).
1545 skip = 1; /* save two point operations in the first
1547 for (i = (num_points ? 520 : 130); i >= 0; --i) {
1550 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1552 /* add multiples of the generator */
1553 if (gen_mul && (i <= 130)) {
1554 bits = get_bit(g_scalar, i + 390) << 3;
1556 bits |= get_bit(g_scalar, i + 260) << 2;
1557 bits |= get_bit(g_scalar, i + 130) << 1;
1558 bits |= get_bit(g_scalar, i);
1560 /* select the point to add, in constant time */
1561 select_point(bits, 16, g_pre_comp, tmp);
1563 /* The 1 argument below is for "mixed" */
1564 point_add(nq[0], nq[1], nq[2],
1565 nq[0], nq[1], nq[2], 1, tmp[0], tmp[1], tmp[2]);
1567 memcpy(nq, tmp, 3 * sizeof(felem));
1572 /* do other additions every 5 doublings */
1573 if (num_points && (i % 5 == 0)) {
1574 /* loop over all scalars */
1575 for (num = 0; num < num_points; ++num) {
1576 bits = get_bit(scalars[num], i + 4) << 5;
1577 bits |= get_bit(scalars[num], i + 3) << 4;
1578 bits |= get_bit(scalars[num], i + 2) << 3;
1579 bits |= get_bit(scalars[num], i + 1) << 2;
1580 bits |= get_bit(scalars[num], i) << 1;
1581 bits |= get_bit(scalars[num], i - 1);
1582 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1585 * select the point to add or subtract, in constant time
1587 select_point(digit, 17, pre_comp[num], tmp);
1588 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative
1590 copy_conditional(tmp[1], tmp[3], (-(limb) sign));
1593 point_add(nq[0], nq[1], nq[2],
1594 nq[0], nq[1], nq[2],
1595 mixed, tmp[0], tmp[1], tmp[2]);
1597 memcpy(nq, tmp, 3 * sizeof(felem));
1603 felem_assign(x_out, nq[0]);
1604 felem_assign(y_out, nq[1]);
1605 felem_assign(z_out, nq[2]);
1608 /* Precomputation for the group generator. */
1609 struct nistp521_pre_comp_st {
1610 felem g_pre_comp[16][3];
1611 CRYPTO_REF_COUNT references;
1612 CRYPTO_RWLOCK *lock;
1615 const EC_METHOD *EC_GFp_nistp521_method(void)
1617 static const EC_METHOD ret = {
1618 EC_FLAGS_DEFAULT_OCT,
1619 NID_X9_62_prime_field,
1620 ec_GFp_nistp521_group_init,
1621 ec_GFp_simple_group_finish,
1622 ec_GFp_simple_group_clear_finish,
1623 ec_GFp_nist_group_copy,
1624 ec_GFp_nistp521_group_set_curve,
1625 ec_GFp_simple_group_get_curve,
1626 ec_GFp_simple_group_get_degree,
1627 ec_group_simple_order_bits,
1628 ec_GFp_simple_group_check_discriminant,
1629 ec_GFp_simple_point_init,
1630 ec_GFp_simple_point_finish,
1631 ec_GFp_simple_point_clear_finish,
1632 ec_GFp_simple_point_copy,
1633 ec_GFp_simple_point_set_to_infinity,
1634 ec_GFp_simple_set_Jprojective_coordinates_GFp,
1635 ec_GFp_simple_get_Jprojective_coordinates_GFp,
1636 ec_GFp_simple_point_set_affine_coordinates,
1637 ec_GFp_nistp521_point_get_affine_coordinates,
1638 0 /* point_set_compressed_coordinates */ ,
1643 ec_GFp_simple_invert,
1644 ec_GFp_simple_is_at_infinity,
1645 ec_GFp_simple_is_on_curve,
1647 ec_GFp_simple_make_affine,
1648 ec_GFp_simple_points_make_affine,
1649 ec_GFp_nistp521_points_mul,
1650 ec_GFp_nistp521_precompute_mult,
1651 ec_GFp_nistp521_have_precompute_mult,
1652 ec_GFp_nist_field_mul,
1653 ec_GFp_nist_field_sqr,
1655 ec_GFp_simple_field_inv,
1656 0 /* field_encode */ ,
1657 0 /* field_decode */ ,
1658 0, /* field_set_to_one */
1659 ec_key_simple_priv2oct,
1660 ec_key_simple_oct2priv,
1661 0, /* set private */
1662 ec_key_simple_generate_key,
1663 ec_key_simple_check_key,
1664 ec_key_simple_generate_public_key,
1667 ecdh_simple_compute_key,
1668 0, /* field_inverse_mod_ord */
1669 0, /* blind_coordinates */
1671 0, /* ladder_step */
1678 /******************************************************************************/
1680 * FUNCTIONS TO MANAGE PRECOMPUTATION
1683 static NISTP521_PRE_COMP *nistp521_pre_comp_new(void)
1685 NISTP521_PRE_COMP *ret = OPENSSL_zalloc(sizeof(*ret));
1688 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1692 ret->references = 1;
1694 ret->lock = CRYPTO_THREAD_lock_new();
1695 if (ret->lock == NULL) {
1696 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1703 NISTP521_PRE_COMP *EC_nistp521_pre_comp_dup(NISTP521_PRE_COMP *p)
1707 CRYPTO_UP_REF(&p->references, &i, p->lock);
1711 void EC_nistp521_pre_comp_free(NISTP521_PRE_COMP *p)
1718 CRYPTO_DOWN_REF(&p->references, &i, p->lock);
1719 REF_PRINT_COUNT("EC_nistp521", x);
1722 REF_ASSERT_ISNT(i < 0);
1724 CRYPTO_THREAD_lock_free(p->lock);
1728 /******************************************************************************/
1730 * OPENSSL EC_METHOD FUNCTIONS
1733 int ec_GFp_nistp521_group_init(EC_GROUP *group)
1736 ret = ec_GFp_simple_group_init(group);
1737 group->a_is_minus3 = 1;
1741 int ec_GFp_nistp521_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1742 const BIGNUM *a, const BIGNUM *b,
1746 BIGNUM *curve_p, *curve_a, *curve_b;
1748 BN_CTX *new_ctx = NULL;
1751 ctx = new_ctx = BN_CTX_new();
1757 curve_p = BN_CTX_get(ctx);
1758 curve_a = BN_CTX_get(ctx);
1759 curve_b = BN_CTX_get(ctx);
1760 if (curve_b == NULL)
1762 BN_bin2bn(nistp521_curve_params[0], sizeof(felem_bytearray), curve_p);
1763 BN_bin2bn(nistp521_curve_params[1], sizeof(felem_bytearray), curve_a);
1764 BN_bin2bn(nistp521_curve_params[2], sizeof(felem_bytearray), curve_b);
1765 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) || (BN_cmp(curve_b, b))) {
1766 ECerr(EC_F_EC_GFP_NISTP521_GROUP_SET_CURVE,
1767 EC_R_WRONG_CURVE_PARAMETERS);
1770 group->field_mod_func = BN_nist_mod_521;
1771 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1775 BN_CTX_free(new_ctx);
1781 * Takes the Jacobian coordinates (X, Y, Z) of a point and returns (X', Y') =
1784 int ec_GFp_nistp521_point_get_affine_coordinates(const EC_GROUP *group,
1785 const EC_POINT *point,
1786 BIGNUM *x, BIGNUM *y,
1789 felem z1, z2, x_in, y_in, x_out, y_out;
1792 if (EC_POINT_is_at_infinity(group, point)) {
1793 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1794 EC_R_POINT_AT_INFINITY);
1797 if ((!BN_to_felem(x_in, point->X)) || (!BN_to_felem(y_in, point->Y)) ||
1798 (!BN_to_felem(z1, point->Z)))
1801 felem_square(tmp, z2);
1802 felem_reduce(z1, tmp);
1803 felem_mul(tmp, x_in, z1);
1804 felem_reduce(x_in, tmp);
1805 felem_contract(x_out, x_in);
1807 if (!felem_to_BN(x, x_out)) {
1808 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1813 felem_mul(tmp, z1, z2);
1814 felem_reduce(z1, tmp);
1815 felem_mul(tmp, y_in, z1);
1816 felem_reduce(y_in, tmp);
1817 felem_contract(y_out, y_in);
1819 if (!felem_to_BN(y, y_out)) {
1820 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1828 /* points below is of size |num|, and tmp_felems is of size |num+1/ */
1829 static void make_points_affine(size_t num, felem points[][3],
1833 * Runs in constant time, unless an input is the point at infinity (which
1834 * normally shouldn't happen).
1836 ec_GFp_nistp_points_make_affine_internal(num,
1840 (void (*)(void *))felem_one,
1842 (void (*)(void *, const void *))
1844 (void (*)(void *, const void *))
1845 felem_square_reduce, (void (*)
1852 (void (*)(void *, const void *))
1854 (void (*)(void *, const void *))
1859 * Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL
1860 * values Result is stored in r (r can equal one of the inputs).
1862 int ec_GFp_nistp521_points_mul(const EC_GROUP *group, EC_POINT *r,
1863 const BIGNUM *scalar, size_t num,
1864 const EC_POINT *points[],
1865 const BIGNUM *scalars[], BN_CTX *ctx)
1870 BIGNUM *x, *y, *z, *tmp_scalar;
1871 felem_bytearray g_secret;
1872 felem_bytearray *secrets = NULL;
1873 felem (*pre_comp)[17][3] = NULL;
1874 felem *tmp_felems = NULL;
1875 felem_bytearray tmp;
1876 unsigned i, num_bytes;
1877 int have_pre_comp = 0;
1878 size_t num_points = num;
1879 felem x_in, y_in, z_in, x_out, y_out, z_out;
1880 NISTP521_PRE_COMP *pre = NULL;
1881 felem(*g_pre_comp)[3] = NULL;
1882 EC_POINT *generator = NULL;
1883 const EC_POINT *p = NULL;
1884 const BIGNUM *p_scalar = NULL;
1887 x = BN_CTX_get(ctx);
1888 y = BN_CTX_get(ctx);
1889 z = BN_CTX_get(ctx);
1890 tmp_scalar = BN_CTX_get(ctx);
1891 if (tmp_scalar == NULL)
1894 if (scalar != NULL) {
1895 pre = group->pre_comp.nistp521;
1897 /* we have precomputation, try to use it */
1898 g_pre_comp = &pre->g_pre_comp[0];
1900 /* try to use the standard precomputation */
1901 g_pre_comp = (felem(*)[3]) gmul;
1902 generator = EC_POINT_new(group);
1903 if (generator == NULL)
1905 /* get the generator from precomputation */
1906 if (!felem_to_BN(x, g_pre_comp[1][0]) ||
1907 !felem_to_BN(y, g_pre_comp[1][1]) ||
1908 !felem_to_BN(z, g_pre_comp[1][2])) {
1909 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1912 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1916 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1917 /* precomputation matches generator */
1921 * we don't have valid precomputation: treat the generator as a
1927 if (num_points > 0) {
1928 if (num_points >= 2) {
1930 * unless we precompute multiples for just one point, converting
1931 * those into affine form is time well spent
1935 secrets = OPENSSL_zalloc(sizeof(*secrets) * num_points);
1936 pre_comp = OPENSSL_zalloc(sizeof(*pre_comp) * num_points);
1939 OPENSSL_malloc(sizeof(*tmp_felems) * (num_points * 17 + 1));
1940 if ((secrets == NULL) || (pre_comp == NULL)
1941 || (mixed && (tmp_felems == NULL))) {
1942 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1947 * we treat NULL scalars as 0, and NULL points as points at infinity,
1948 * i.e., they contribute nothing to the linear combination
1950 for (i = 0; i < num_points; ++i) {
1953 * we didn't have a valid precomputation, so we pick the
1957 p = EC_GROUP_get0_generator(group);
1960 /* the i^th point */
1963 p_scalar = scalars[i];
1965 if ((p_scalar != NULL) && (p != NULL)) {
1966 /* reduce scalar to 0 <= scalar < 2^521 */
1967 if ((BN_num_bits(p_scalar) > 521)
1968 || (BN_is_negative(p_scalar))) {
1970 * this is an unusual input, and we don't guarantee
1973 if (!BN_nnmod(tmp_scalar, p_scalar, group->order, ctx)) {
1974 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1977 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1979 num_bytes = BN_bn2bin(p_scalar, tmp);
1980 flip_endian(secrets[i], tmp, num_bytes);
1981 /* precompute multiples */
1982 if ((!BN_to_felem(x_out, p->X)) ||
1983 (!BN_to_felem(y_out, p->Y)) ||
1984 (!BN_to_felem(z_out, p->Z)))
1986 memcpy(pre_comp[i][1][0], x_out, sizeof(felem));
1987 memcpy(pre_comp[i][1][1], y_out, sizeof(felem));
1988 memcpy(pre_comp[i][1][2], z_out, sizeof(felem));
1989 for (j = 2; j <= 16; ++j) {
1991 point_add(pre_comp[i][j][0], pre_comp[i][j][1],
1992 pre_comp[i][j][2], pre_comp[i][1][0],
1993 pre_comp[i][1][1], pre_comp[i][1][2], 0,
1994 pre_comp[i][j - 1][0],
1995 pre_comp[i][j - 1][1],
1996 pre_comp[i][j - 1][2]);
1998 point_double(pre_comp[i][j][0], pre_comp[i][j][1],
1999 pre_comp[i][j][2], pre_comp[i][j / 2][0],
2000 pre_comp[i][j / 2][1],
2001 pre_comp[i][j / 2][2]);
2007 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
2010 /* the scalar for the generator */
2011 if ((scalar != NULL) && (have_pre_comp)) {
2012 memset(g_secret, 0, sizeof(g_secret));
2013 /* reduce scalar to 0 <= scalar < 2^521 */
2014 if ((BN_num_bits(scalar) > 521) || (BN_is_negative(scalar))) {
2016 * this is an unusual input, and we don't guarantee
2019 if (!BN_nnmod(tmp_scalar, scalar, group->order, ctx)) {
2020 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
2023 num_bytes = BN_bn2bin(tmp_scalar, tmp);
2025 num_bytes = BN_bn2bin(scalar, tmp);
2026 flip_endian(g_secret, tmp, num_bytes);
2027 /* do the multiplication with generator precomputation */
2028 batch_mul(x_out, y_out, z_out,
2029 (const felem_bytearray(*))secrets, num_points,
2031 mixed, (const felem(*)[17][3])pre_comp,
2032 (const felem(*)[3])g_pre_comp);
2034 /* do the multiplication without generator precomputation */
2035 batch_mul(x_out, y_out, z_out,
2036 (const felem_bytearray(*))secrets, num_points,
2037 NULL, mixed, (const felem(*)[17][3])pre_comp, NULL);
2038 /* reduce the output to its unique minimal representation */
2039 felem_contract(x_in, x_out);
2040 felem_contract(y_in, y_out);
2041 felem_contract(z_in, z_out);
2042 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
2043 (!felem_to_BN(z, z_in))) {
2044 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
2047 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
2051 EC_POINT_free(generator);
2052 OPENSSL_free(secrets);
2053 OPENSSL_free(pre_comp);
2054 OPENSSL_free(tmp_felems);
2058 int ec_GFp_nistp521_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
2061 NISTP521_PRE_COMP *pre = NULL;
2064 EC_POINT *generator = NULL;
2065 felem tmp_felems[16];
2067 BN_CTX *new_ctx = NULL;
2070 /* throw away old precomputation */
2071 EC_pre_comp_free(group);
2075 ctx = new_ctx = BN_CTX_new();
2081 x = BN_CTX_get(ctx);
2082 y = BN_CTX_get(ctx);
2085 /* get the generator */
2086 if (group->generator == NULL)
2088 generator = EC_POINT_new(group);
2089 if (generator == NULL)
2091 BN_bin2bn(nistp521_curve_params[3], sizeof(felem_bytearray), x);
2092 BN_bin2bn(nistp521_curve_params[4], sizeof(felem_bytearray), y);
2093 if (!EC_POINT_set_affine_coordinates(group, generator, x, y, ctx))
2095 if ((pre = nistp521_pre_comp_new()) == NULL)
2098 * if the generator is the standard one, use built-in precomputation
2100 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) {
2101 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
2104 if ((!BN_to_felem(pre->g_pre_comp[1][0], group->generator->X)) ||
2105 (!BN_to_felem(pre->g_pre_comp[1][1], group->generator->Y)) ||
2106 (!BN_to_felem(pre->g_pre_comp[1][2], group->generator->Z)))
2108 /* compute 2^130*G, 2^260*G, 2^390*G */
2109 for (i = 1; i <= 4; i <<= 1) {
2110 point_double(pre->g_pre_comp[2 * i][0], pre->g_pre_comp[2 * i][1],
2111 pre->g_pre_comp[2 * i][2], pre->g_pre_comp[i][0],
2112 pre->g_pre_comp[i][1], pre->g_pre_comp[i][2]);
2113 for (j = 0; j < 129; ++j) {
2114 point_double(pre->g_pre_comp[2 * i][0],
2115 pre->g_pre_comp[2 * i][1],
2116 pre->g_pre_comp[2 * i][2],
2117 pre->g_pre_comp[2 * i][0],
2118 pre->g_pre_comp[2 * i][1],
2119 pre->g_pre_comp[2 * i][2]);
2122 /* g_pre_comp[0] is the point at infinity */
2123 memset(pre->g_pre_comp[0], 0, sizeof(pre->g_pre_comp[0]));
2124 /* the remaining multiples */
2125 /* 2^130*G + 2^260*G */
2126 point_add(pre->g_pre_comp[6][0], pre->g_pre_comp[6][1],
2127 pre->g_pre_comp[6][2], pre->g_pre_comp[4][0],
2128 pre->g_pre_comp[4][1], pre->g_pre_comp[4][2],
2129 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2130 pre->g_pre_comp[2][2]);
2131 /* 2^130*G + 2^390*G */
2132 point_add(pre->g_pre_comp[10][0], pre->g_pre_comp[10][1],
2133 pre->g_pre_comp[10][2], pre->g_pre_comp[8][0],
2134 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
2135 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2136 pre->g_pre_comp[2][2]);
2137 /* 2^260*G + 2^390*G */
2138 point_add(pre->g_pre_comp[12][0], pre->g_pre_comp[12][1],
2139 pre->g_pre_comp[12][2], pre->g_pre_comp[8][0],
2140 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
2141 0, pre->g_pre_comp[4][0], pre->g_pre_comp[4][1],
2142 pre->g_pre_comp[4][2]);
2143 /* 2^130*G + 2^260*G + 2^390*G */
2144 point_add(pre->g_pre_comp[14][0], pre->g_pre_comp[14][1],
2145 pre->g_pre_comp[14][2], pre->g_pre_comp[12][0],
2146 pre->g_pre_comp[12][1], pre->g_pre_comp[12][2],
2147 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2148 pre->g_pre_comp[2][2]);
2149 for (i = 1; i < 8; ++i) {
2150 /* odd multiples: add G */
2151 point_add(pre->g_pre_comp[2 * i + 1][0],
2152 pre->g_pre_comp[2 * i + 1][1],
2153 pre->g_pre_comp[2 * i + 1][2], pre->g_pre_comp[2 * i][0],
2154 pre->g_pre_comp[2 * i][1], pre->g_pre_comp[2 * i][2], 0,
2155 pre->g_pre_comp[1][0], pre->g_pre_comp[1][1],
2156 pre->g_pre_comp[1][2]);
2158 make_points_affine(15, &(pre->g_pre_comp[1]), tmp_felems);
2161 SETPRECOMP(group, nistp521, pre);
2166 EC_POINT_free(generator);
2168 BN_CTX_free(new_ctx);
2170 EC_nistp521_pre_comp_free(pre);
2174 int ec_GFp_nistp521_have_precompute_mult(const EC_GROUP *group)
2176 return HAVEPRECOMP(group, nistp521);