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 * ECDSA low level APIs are deprecated for public use, but still ok for
30 #include "internal/deprecated.h"
33 * A 64-bit implementation of the NIST P-521 elliptic curve point multiplication
35 * OpenSSL integration was taken from Emilia Kasper's work in ecp_nistp224.c.
36 * Otherwise based on Emilia's P224 work, which was inspired by my curve25519
37 * work which got its smarts from Daniel J. Bernstein's work on the same.
40 #include <openssl/e_os2.h>
41 #ifdef OPENSSL_NO_EC_NISTP_64_GCC_128
42 NON_EMPTY_TRANSLATION_UNIT
46 # include <openssl/err.h>
47 # include "ec_local.h"
49 # if defined(__SIZEOF_INT128__) && __SIZEOF_INT128__==16
50 /* even with gcc, the typedef won't work for 32-bit platforms */
51 typedef __uint128_t uint128_t; /* nonstandard; implemented by gcc on 64-bit
54 # error "Your compiler doesn't appear to support 128-bit integer types"
61 * The underlying field. P521 operates over GF(2^521-1). We can serialise an
62 * element of this field into 66 bytes where the most significant byte
63 * contains only a single bit. We call this an felem_bytearray.
66 typedef u8 felem_bytearray[66];
69 * These are the parameters of P521, taken from FIPS 186-3, section D.1.2.5.
70 * These values are big-endian.
72 static const felem_bytearray nistp521_curve_params[5] = {
73 {0x01, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, /* p */
74 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
75 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
76 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
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,
82 {0x01, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, /* a = -3 */
83 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
84 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
85 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
86 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
87 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
88 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
89 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
91 {0x00, 0x51, 0x95, 0x3e, 0xb9, 0x61, 0x8e, 0x1c, /* b */
92 0x9a, 0x1f, 0x92, 0x9a, 0x21, 0xa0, 0xb6, 0x85,
93 0x40, 0xee, 0xa2, 0xda, 0x72, 0x5b, 0x99, 0xb3,
94 0x15, 0xf3, 0xb8, 0xb4, 0x89, 0x91, 0x8e, 0xf1,
95 0x09, 0xe1, 0x56, 0x19, 0x39, 0x51, 0xec, 0x7e,
96 0x93, 0x7b, 0x16, 0x52, 0xc0, 0xbd, 0x3b, 0xb1,
97 0xbf, 0x07, 0x35, 0x73, 0xdf, 0x88, 0x3d, 0x2c,
98 0x34, 0xf1, 0xef, 0x45, 0x1f, 0xd4, 0x6b, 0x50,
100 {0x00, 0xc6, 0x85, 0x8e, 0x06, 0xb7, 0x04, 0x04, /* x */
101 0xe9, 0xcd, 0x9e, 0x3e, 0xcb, 0x66, 0x23, 0x95,
102 0xb4, 0x42, 0x9c, 0x64, 0x81, 0x39, 0x05, 0x3f,
103 0xb5, 0x21, 0xf8, 0x28, 0xaf, 0x60, 0x6b, 0x4d,
104 0x3d, 0xba, 0xa1, 0x4b, 0x5e, 0x77, 0xef, 0xe7,
105 0x59, 0x28, 0xfe, 0x1d, 0xc1, 0x27, 0xa2, 0xff,
106 0xa8, 0xde, 0x33, 0x48, 0xb3, 0xc1, 0x85, 0x6a,
107 0x42, 0x9b, 0xf9, 0x7e, 0x7e, 0x31, 0xc2, 0xe5,
109 {0x01, 0x18, 0x39, 0x29, 0x6a, 0x78, 0x9a, 0x3b, /* y */
110 0xc0, 0x04, 0x5c, 0x8a, 0x5f, 0xb4, 0x2c, 0x7d,
111 0x1b, 0xd9, 0x98, 0xf5, 0x44, 0x49, 0x57, 0x9b,
112 0x44, 0x68, 0x17, 0xaf, 0xbd, 0x17, 0x27, 0x3e,
113 0x66, 0x2c, 0x97, 0xee, 0x72, 0x99, 0x5e, 0xf4,
114 0x26, 0x40, 0xc5, 0x50, 0xb9, 0x01, 0x3f, 0xad,
115 0x07, 0x61, 0x35, 0x3c, 0x70, 0x86, 0xa2, 0x72,
116 0xc2, 0x40, 0x88, 0xbe, 0x94, 0x76, 0x9f, 0xd1,
121 * The representation of field elements.
122 * ------------------------------------
124 * We represent field elements with nine values. These values are either 64 or
125 * 128 bits and the field element represented is:
126 * v[0]*2^0 + v[1]*2^58 + v[2]*2^116 + ... + v[8]*2^464 (mod p)
127 * Each of the nine values is called a 'limb'. Since the limbs are spaced only
128 * 58 bits apart, but are greater than 58 bits in length, the most significant
129 * bits of each limb overlap with the least significant bits of the next.
131 * A field element with 64-bit limbs is an 'felem'. One with 128-bit limbs is a
136 typedef uint64_t limb;
137 typedef limb felem[NLIMBS];
138 typedef uint128_t largefelem[NLIMBS];
140 static const limb bottom57bits = 0x1ffffffffffffff;
141 static const limb bottom58bits = 0x3ffffffffffffff;
144 * bin66_to_felem takes a little-endian byte array and converts it into felem
145 * form. This assumes that the CPU is little-endian.
147 static void bin66_to_felem(felem out, const u8 in[66])
149 out[0] = (*((limb *) & in[0])) & bottom58bits;
150 out[1] = (*((limb *) & in[7]) >> 2) & bottom58bits;
151 out[2] = (*((limb *) & in[14]) >> 4) & bottom58bits;
152 out[3] = (*((limb *) & in[21]) >> 6) & bottom58bits;
153 out[4] = (*((limb *) & in[29])) & bottom58bits;
154 out[5] = (*((limb *) & in[36]) >> 2) & bottom58bits;
155 out[6] = (*((limb *) & in[43]) >> 4) & bottom58bits;
156 out[7] = (*((limb *) & in[50]) >> 6) & bottom58bits;
157 out[8] = (*((limb *) & in[58])) & bottom57bits;
161 * felem_to_bin66 takes an felem and serialises into a little endian, 66 byte
162 * array. This assumes that the CPU is little-endian.
164 static void felem_to_bin66(u8 out[66], const felem in)
167 (*((limb *) & out[0])) = in[0];
168 (*((limb *) & out[7])) |= in[1] << 2;
169 (*((limb *) & out[14])) |= in[2] << 4;
170 (*((limb *) & out[21])) |= in[3] << 6;
171 (*((limb *) & out[29])) = in[4];
172 (*((limb *) & out[36])) |= in[5] << 2;
173 (*((limb *) & out[43])) |= in[6] << 4;
174 (*((limb *) & out[50])) |= in[7] << 6;
175 (*((limb *) & out[58])) = in[8];
178 /* BN_to_felem converts an OpenSSL BIGNUM into an felem */
179 static int BN_to_felem(felem out, const BIGNUM *bn)
181 felem_bytearray b_out;
184 if (BN_is_negative(bn)) {
185 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
188 num_bytes = BN_bn2lebinpad(bn, b_out, sizeof(b_out));
190 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
193 bin66_to_felem(out, b_out);
197 /* felem_to_BN converts an felem into an OpenSSL BIGNUM */
198 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
200 felem_bytearray b_out;
201 felem_to_bin66(b_out, in);
202 return BN_lebin2bn(b_out, sizeof(b_out), out);
210 static void felem_one(felem out)
223 static void felem_assign(felem out, const felem in)
236 /* felem_sum64 sets out = out + in. */
237 static void felem_sum64(felem out, const felem in)
250 /* felem_scalar sets out = in * scalar */
251 static void felem_scalar(felem out, const felem in, limb scalar)
253 out[0] = in[0] * scalar;
254 out[1] = in[1] * scalar;
255 out[2] = in[2] * scalar;
256 out[3] = in[3] * scalar;
257 out[4] = in[4] * scalar;
258 out[5] = in[5] * scalar;
259 out[6] = in[6] * scalar;
260 out[7] = in[7] * scalar;
261 out[8] = in[8] * scalar;
264 /* felem_scalar64 sets out = out * scalar */
265 static void felem_scalar64(felem out, limb scalar)
278 /* felem_scalar128 sets out = out * scalar */
279 static void felem_scalar128(largefelem out, limb scalar)
293 * felem_neg sets |out| to |-in|
295 * in[i] < 2^59 + 2^14
299 static void felem_neg(felem out, const felem in)
301 /* In order to prevent underflow, we subtract from 0 mod p. */
302 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5);
303 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4);
305 out[0] = two62m3 - in[0];
306 out[1] = two62m2 - in[1];
307 out[2] = two62m2 - in[2];
308 out[3] = two62m2 - in[3];
309 out[4] = two62m2 - in[4];
310 out[5] = two62m2 - in[5];
311 out[6] = two62m2 - in[6];
312 out[7] = two62m2 - in[7];
313 out[8] = two62m2 - in[8];
317 * felem_diff64 subtracts |in| from |out|
319 * in[i] < 2^59 + 2^14
321 * out[i] < out[i] + 2^62
323 static void felem_diff64(felem out, const felem in)
326 * In order to prevent underflow, we add 0 mod p before subtracting.
328 static const limb two62m3 = (((limb) 1) << 62) - (((limb) 1) << 5);
329 static const limb two62m2 = (((limb) 1) << 62) - (((limb) 1) << 4);
331 out[0] += two62m3 - in[0];
332 out[1] += two62m2 - in[1];
333 out[2] += two62m2 - in[2];
334 out[3] += two62m2 - in[3];
335 out[4] += two62m2 - in[4];
336 out[5] += two62m2 - in[5];
337 out[6] += two62m2 - in[6];
338 out[7] += two62m2 - in[7];
339 out[8] += two62m2 - in[8];
343 * felem_diff_128_64 subtracts |in| from |out|
345 * in[i] < 2^62 + 2^17
347 * out[i] < out[i] + 2^63
349 static void felem_diff_128_64(largefelem out, const felem in)
352 * In order to prevent underflow, we add 64p mod p (which is equivalent
353 * to 0 mod p) before subtracting. p is 2^521 - 1, i.e. in binary a 521
354 * digit number with all bits set to 1. See "The representation of field
355 * elements" comment above for a description of how limbs are used to
356 * represent a number. 64p is represented with 8 limbs containing a number
357 * with 58 bits set and one limb with a number with 57 bits set.
359 static const limb two63m6 = (((limb) 1) << 63) - (((limb) 1) << 6);
360 static const limb two63m5 = (((limb) 1) << 63) - (((limb) 1) << 5);
362 out[0] += two63m6 - in[0];
363 out[1] += two63m5 - in[1];
364 out[2] += two63m5 - in[2];
365 out[3] += two63m5 - in[3];
366 out[4] += two63m5 - in[4];
367 out[5] += two63m5 - in[5];
368 out[6] += two63m5 - in[6];
369 out[7] += two63m5 - in[7];
370 out[8] += two63m5 - in[8];
374 * felem_diff_128_64 subtracts |in| from |out|
378 * out[i] < out[i] + 2^127 - 2^69
380 static void felem_diff128(largefelem out, const largefelem in)
383 * In order to prevent underflow, we add 0 mod p before subtracting.
385 static const uint128_t two127m70 =
386 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 70);
387 static const uint128_t two127m69 =
388 (((uint128_t) 1) << 127) - (((uint128_t) 1) << 69);
390 out[0] += (two127m70 - in[0]);
391 out[1] += (two127m69 - in[1]);
392 out[2] += (two127m69 - in[2]);
393 out[3] += (two127m69 - in[3]);
394 out[4] += (two127m69 - in[4]);
395 out[5] += (two127m69 - in[5]);
396 out[6] += (two127m69 - in[6]);
397 out[7] += (two127m69 - in[7]);
398 out[8] += (two127m69 - in[8]);
402 * felem_square sets |out| = |in|^2
406 * out[i] < 17 * max(in[i]) * max(in[i])
408 static void felem_square(largefelem out, const felem in)
411 felem_scalar(inx2, in, 2);
412 felem_scalar(inx4, in, 4);
415 * We have many cases were we want to do
418 * This is obviously just
420 * However, rather than do the doubling on the 128 bit result, we
421 * double one of the inputs to the multiplication by reading from
425 out[0] = ((uint128_t) in[0]) * in[0];
426 out[1] = ((uint128_t) in[0]) * inx2[1];
427 out[2] = ((uint128_t) in[0]) * inx2[2] + ((uint128_t) in[1]) * in[1];
428 out[3] = ((uint128_t) in[0]) * inx2[3] + ((uint128_t) in[1]) * inx2[2];
429 out[4] = ((uint128_t) in[0]) * inx2[4] +
430 ((uint128_t) in[1]) * inx2[3] + ((uint128_t) in[2]) * in[2];
431 out[5] = ((uint128_t) in[0]) * inx2[5] +
432 ((uint128_t) in[1]) * inx2[4] + ((uint128_t) in[2]) * inx2[3];
433 out[6] = ((uint128_t) in[0]) * inx2[6] +
434 ((uint128_t) in[1]) * inx2[5] +
435 ((uint128_t) in[2]) * inx2[4] + ((uint128_t) in[3]) * in[3];
436 out[7] = ((uint128_t) in[0]) * inx2[7] +
437 ((uint128_t) in[1]) * inx2[6] +
438 ((uint128_t) in[2]) * inx2[5] + ((uint128_t) in[3]) * inx2[4];
439 out[8] = ((uint128_t) in[0]) * inx2[8] +
440 ((uint128_t) in[1]) * inx2[7] +
441 ((uint128_t) in[2]) * inx2[6] +
442 ((uint128_t) in[3]) * inx2[5] + ((uint128_t) in[4]) * in[4];
445 * The remaining limbs fall above 2^521, with the first falling at 2^522.
446 * They correspond to locations one bit up from the limbs produced above
447 * so we would have to multiply by two to align them. Again, rather than
448 * operate on the 128-bit result, we double one of the inputs to the
449 * multiplication. If we want to double for both this reason, and the
450 * reason above, then we end up multiplying by four.
454 out[0] += ((uint128_t) in[1]) * inx4[8] +
455 ((uint128_t) in[2]) * inx4[7] +
456 ((uint128_t) in[3]) * inx4[6] + ((uint128_t) in[4]) * inx4[5];
459 out[1] += ((uint128_t) in[2]) * inx4[8] +
460 ((uint128_t) in[3]) * inx4[7] +
461 ((uint128_t) in[4]) * inx4[6] + ((uint128_t) in[5]) * inx2[5];
464 out[2] += ((uint128_t) in[3]) * inx4[8] +
465 ((uint128_t) in[4]) * inx4[7] + ((uint128_t) in[5]) * inx4[6];
468 out[3] += ((uint128_t) in[4]) * inx4[8] +
469 ((uint128_t) in[5]) * inx4[7] + ((uint128_t) in[6]) * inx2[6];
472 out[4] += ((uint128_t) in[5]) * inx4[8] + ((uint128_t) in[6]) * inx4[7];
475 out[5] += ((uint128_t) in[6]) * inx4[8] + ((uint128_t) in[7]) * inx2[7];
478 out[6] += ((uint128_t) in[7]) * inx4[8];
481 out[7] += ((uint128_t) in[8]) * inx2[8];
485 * felem_mul sets |out| = |in1| * |in2|
490 * out[i] < 17 * max(in1[i]) * max(in2[i])
492 static void felem_mul(largefelem out, const felem in1, const felem in2)
495 felem_scalar(in2x2, in2, 2);
497 out[0] = ((uint128_t) in1[0]) * in2[0];
499 out[1] = ((uint128_t) in1[0]) * in2[1] +
500 ((uint128_t) in1[1]) * in2[0];
502 out[2] = ((uint128_t) in1[0]) * in2[2] +
503 ((uint128_t) in1[1]) * in2[1] +
504 ((uint128_t) in1[2]) * in2[0];
506 out[3] = ((uint128_t) in1[0]) * in2[3] +
507 ((uint128_t) in1[1]) * in2[2] +
508 ((uint128_t) in1[2]) * in2[1] +
509 ((uint128_t) in1[3]) * in2[0];
511 out[4] = ((uint128_t) in1[0]) * in2[4] +
512 ((uint128_t) in1[1]) * in2[3] +
513 ((uint128_t) in1[2]) * in2[2] +
514 ((uint128_t) in1[3]) * in2[1] +
515 ((uint128_t) in1[4]) * in2[0];
517 out[5] = ((uint128_t) in1[0]) * in2[5] +
518 ((uint128_t) in1[1]) * in2[4] +
519 ((uint128_t) in1[2]) * in2[3] +
520 ((uint128_t) in1[3]) * in2[2] +
521 ((uint128_t) in1[4]) * in2[1] +
522 ((uint128_t) in1[5]) * in2[0];
524 out[6] = ((uint128_t) in1[0]) * in2[6] +
525 ((uint128_t) in1[1]) * in2[5] +
526 ((uint128_t) in1[2]) * in2[4] +
527 ((uint128_t) in1[3]) * in2[3] +
528 ((uint128_t) in1[4]) * in2[2] +
529 ((uint128_t) in1[5]) * in2[1] +
530 ((uint128_t) in1[6]) * in2[0];
532 out[7] = ((uint128_t) in1[0]) * in2[7] +
533 ((uint128_t) in1[1]) * in2[6] +
534 ((uint128_t) in1[2]) * in2[5] +
535 ((uint128_t) in1[3]) * in2[4] +
536 ((uint128_t) in1[4]) * in2[3] +
537 ((uint128_t) in1[5]) * in2[2] +
538 ((uint128_t) in1[6]) * in2[1] +
539 ((uint128_t) in1[7]) * in2[0];
541 out[8] = ((uint128_t) in1[0]) * in2[8] +
542 ((uint128_t) in1[1]) * in2[7] +
543 ((uint128_t) in1[2]) * in2[6] +
544 ((uint128_t) in1[3]) * in2[5] +
545 ((uint128_t) in1[4]) * in2[4] +
546 ((uint128_t) in1[5]) * in2[3] +
547 ((uint128_t) in1[6]) * in2[2] +
548 ((uint128_t) in1[7]) * in2[1] +
549 ((uint128_t) in1[8]) * in2[0];
551 /* See comment in felem_square about the use of in2x2 here */
553 out[0] += ((uint128_t) in1[1]) * in2x2[8] +
554 ((uint128_t) in1[2]) * in2x2[7] +
555 ((uint128_t) in1[3]) * in2x2[6] +
556 ((uint128_t) in1[4]) * in2x2[5] +
557 ((uint128_t) in1[5]) * in2x2[4] +
558 ((uint128_t) in1[6]) * in2x2[3] +
559 ((uint128_t) in1[7]) * in2x2[2] +
560 ((uint128_t) in1[8]) * in2x2[1];
562 out[1] += ((uint128_t) in1[2]) * in2x2[8] +
563 ((uint128_t) in1[3]) * in2x2[7] +
564 ((uint128_t) in1[4]) * in2x2[6] +
565 ((uint128_t) in1[5]) * in2x2[5] +
566 ((uint128_t) in1[6]) * in2x2[4] +
567 ((uint128_t) in1[7]) * in2x2[3] +
568 ((uint128_t) in1[8]) * in2x2[2];
570 out[2] += ((uint128_t) in1[3]) * in2x2[8] +
571 ((uint128_t) in1[4]) * in2x2[7] +
572 ((uint128_t) in1[5]) * in2x2[6] +
573 ((uint128_t) in1[6]) * in2x2[5] +
574 ((uint128_t) in1[7]) * in2x2[4] +
575 ((uint128_t) in1[8]) * in2x2[3];
577 out[3] += ((uint128_t) in1[4]) * in2x2[8] +
578 ((uint128_t) in1[5]) * in2x2[7] +
579 ((uint128_t) in1[6]) * in2x2[6] +
580 ((uint128_t) in1[7]) * in2x2[5] +
581 ((uint128_t) in1[8]) * in2x2[4];
583 out[4] += ((uint128_t) in1[5]) * in2x2[8] +
584 ((uint128_t) in1[6]) * in2x2[7] +
585 ((uint128_t) in1[7]) * in2x2[6] +
586 ((uint128_t) in1[8]) * in2x2[5];
588 out[5] += ((uint128_t) in1[6]) * in2x2[8] +
589 ((uint128_t) in1[7]) * in2x2[7] +
590 ((uint128_t) in1[8]) * in2x2[6];
592 out[6] += ((uint128_t) in1[7]) * in2x2[8] +
593 ((uint128_t) in1[8]) * in2x2[7];
595 out[7] += ((uint128_t) in1[8]) * in2x2[8];
598 static const limb bottom52bits = 0xfffffffffffff;
601 * felem_reduce converts a largefelem to an felem.
605 * out[i] < 2^59 + 2^14
607 static void felem_reduce(felem out, const largefelem in)
609 u64 overflow1, overflow2;
611 out[0] = ((limb) in[0]) & bottom58bits;
612 out[1] = ((limb) in[1]) & bottom58bits;
613 out[2] = ((limb) in[2]) & bottom58bits;
614 out[3] = ((limb) in[3]) & bottom58bits;
615 out[4] = ((limb) in[4]) & bottom58bits;
616 out[5] = ((limb) in[5]) & bottom58bits;
617 out[6] = ((limb) in[6]) & bottom58bits;
618 out[7] = ((limb) in[7]) & bottom58bits;
619 out[8] = ((limb) in[8]) & bottom58bits;
623 out[1] += ((limb) in[0]) >> 58;
624 out[1] += (((limb) (in[0] >> 64)) & bottom52bits) << 6;
626 * out[1] < 2^58 + 2^6 + 2^58
629 out[2] += ((limb) (in[0] >> 64)) >> 52;
631 out[2] += ((limb) in[1]) >> 58;
632 out[2] += (((limb) (in[1] >> 64)) & bottom52bits) << 6;
633 out[3] += ((limb) (in[1] >> 64)) >> 52;
635 out[3] += ((limb) in[2]) >> 58;
636 out[3] += (((limb) (in[2] >> 64)) & bottom52bits) << 6;
637 out[4] += ((limb) (in[2] >> 64)) >> 52;
639 out[4] += ((limb) in[3]) >> 58;
640 out[4] += (((limb) (in[3] >> 64)) & bottom52bits) << 6;
641 out[5] += ((limb) (in[3] >> 64)) >> 52;
643 out[5] += ((limb) in[4]) >> 58;
644 out[5] += (((limb) (in[4] >> 64)) & bottom52bits) << 6;
645 out[6] += ((limb) (in[4] >> 64)) >> 52;
647 out[6] += ((limb) in[5]) >> 58;
648 out[6] += (((limb) (in[5] >> 64)) & bottom52bits) << 6;
649 out[7] += ((limb) (in[5] >> 64)) >> 52;
651 out[7] += ((limb) in[6]) >> 58;
652 out[7] += (((limb) (in[6] >> 64)) & bottom52bits) << 6;
653 out[8] += ((limb) (in[6] >> 64)) >> 52;
655 out[8] += ((limb) in[7]) >> 58;
656 out[8] += (((limb) (in[7] >> 64)) & bottom52bits) << 6;
658 * out[x > 1] < 2^58 + 2^6 + 2^58 + 2^12
661 overflow1 = ((limb) (in[7] >> 64)) >> 52;
663 overflow1 += ((limb) in[8]) >> 58;
664 overflow1 += (((limb) (in[8] >> 64)) & bottom52bits) << 6;
665 overflow2 = ((limb) (in[8] >> 64)) >> 52;
667 overflow1 <<= 1; /* overflow1 < 2^13 + 2^7 + 2^59 */
668 overflow2 <<= 1; /* overflow2 < 2^13 */
670 out[0] += overflow1; /* out[0] < 2^60 */
671 out[1] += overflow2; /* out[1] < 2^59 + 2^6 + 2^13 */
673 out[1] += out[0] >> 58;
674 out[0] &= bottom58bits;
677 * out[1] < 2^59 + 2^6 + 2^13 + 2^2
682 static void felem_square_reduce(felem out, const felem in)
685 felem_square(tmp, in);
686 felem_reduce(out, tmp);
689 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
692 felem_mul(tmp, in1, in2);
693 felem_reduce(out, tmp);
697 * felem_inv calculates |out| = |in|^{-1}
699 * Based on Fermat's Little Theorem:
701 * a^{p-1} = 1 (mod p)
702 * a^{p-2} = a^{-1} (mod p)
704 static void felem_inv(felem out, const felem in)
706 felem ftmp, ftmp2, ftmp3, ftmp4;
710 felem_square(tmp, in);
711 felem_reduce(ftmp, tmp); /* 2^1 */
712 felem_mul(tmp, in, ftmp);
713 felem_reduce(ftmp, tmp); /* 2^2 - 2^0 */
714 felem_assign(ftmp2, ftmp);
715 felem_square(tmp, ftmp);
716 felem_reduce(ftmp, tmp); /* 2^3 - 2^1 */
717 felem_mul(tmp, in, ftmp);
718 felem_reduce(ftmp, tmp); /* 2^3 - 2^0 */
719 felem_square(tmp, ftmp);
720 felem_reduce(ftmp, tmp); /* 2^4 - 2^1 */
722 felem_square(tmp, ftmp2);
723 felem_reduce(ftmp3, tmp); /* 2^3 - 2^1 */
724 felem_square(tmp, ftmp3);
725 felem_reduce(ftmp3, tmp); /* 2^4 - 2^2 */
726 felem_mul(tmp, ftmp3, ftmp2);
727 felem_reduce(ftmp3, tmp); /* 2^4 - 2^0 */
729 felem_assign(ftmp2, ftmp3);
730 felem_square(tmp, ftmp3);
731 felem_reduce(ftmp3, tmp); /* 2^5 - 2^1 */
732 felem_square(tmp, ftmp3);
733 felem_reduce(ftmp3, tmp); /* 2^6 - 2^2 */
734 felem_square(tmp, ftmp3);
735 felem_reduce(ftmp3, tmp); /* 2^7 - 2^3 */
736 felem_square(tmp, ftmp3);
737 felem_reduce(ftmp3, tmp); /* 2^8 - 2^4 */
738 felem_assign(ftmp4, ftmp3);
739 felem_mul(tmp, ftmp3, ftmp);
740 felem_reduce(ftmp4, tmp); /* 2^8 - 2^1 */
741 felem_square(tmp, ftmp4);
742 felem_reduce(ftmp4, tmp); /* 2^9 - 2^2 */
743 felem_mul(tmp, ftmp3, ftmp2);
744 felem_reduce(ftmp3, tmp); /* 2^8 - 2^0 */
745 felem_assign(ftmp2, ftmp3);
747 for (i = 0; i < 8; i++) {
748 felem_square(tmp, ftmp3);
749 felem_reduce(ftmp3, tmp); /* 2^16 - 2^8 */
751 felem_mul(tmp, ftmp3, ftmp2);
752 felem_reduce(ftmp3, tmp); /* 2^16 - 2^0 */
753 felem_assign(ftmp2, ftmp3);
755 for (i = 0; i < 16; i++) {
756 felem_square(tmp, ftmp3);
757 felem_reduce(ftmp3, tmp); /* 2^32 - 2^16 */
759 felem_mul(tmp, ftmp3, ftmp2);
760 felem_reduce(ftmp3, tmp); /* 2^32 - 2^0 */
761 felem_assign(ftmp2, ftmp3);
763 for (i = 0; i < 32; i++) {
764 felem_square(tmp, ftmp3);
765 felem_reduce(ftmp3, tmp); /* 2^64 - 2^32 */
767 felem_mul(tmp, ftmp3, ftmp2);
768 felem_reduce(ftmp3, tmp); /* 2^64 - 2^0 */
769 felem_assign(ftmp2, ftmp3);
771 for (i = 0; i < 64; i++) {
772 felem_square(tmp, ftmp3);
773 felem_reduce(ftmp3, tmp); /* 2^128 - 2^64 */
775 felem_mul(tmp, ftmp3, ftmp2);
776 felem_reduce(ftmp3, tmp); /* 2^128 - 2^0 */
777 felem_assign(ftmp2, ftmp3);
779 for (i = 0; i < 128; i++) {
780 felem_square(tmp, ftmp3);
781 felem_reduce(ftmp3, tmp); /* 2^256 - 2^128 */
783 felem_mul(tmp, ftmp3, ftmp2);
784 felem_reduce(ftmp3, tmp); /* 2^256 - 2^0 */
785 felem_assign(ftmp2, ftmp3);
787 for (i = 0; i < 256; i++) {
788 felem_square(tmp, ftmp3);
789 felem_reduce(ftmp3, tmp); /* 2^512 - 2^256 */
791 felem_mul(tmp, ftmp3, ftmp2);
792 felem_reduce(ftmp3, tmp); /* 2^512 - 2^0 */
794 for (i = 0; i < 9; i++) {
795 felem_square(tmp, ftmp3);
796 felem_reduce(ftmp3, tmp); /* 2^521 - 2^9 */
798 felem_mul(tmp, ftmp3, ftmp4);
799 felem_reduce(ftmp3, tmp); /* 2^512 - 2^2 */
800 felem_mul(tmp, ftmp3, in);
801 felem_reduce(out, tmp); /* 2^512 - 3 */
804 /* This is 2^521-1, expressed as an felem */
805 static const felem kPrime = {
806 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff,
807 0x03ffffffffffffff, 0x03ffffffffffffff, 0x03ffffffffffffff,
808 0x03ffffffffffffff, 0x03ffffffffffffff, 0x01ffffffffffffff
812 * felem_is_zero returns a limb with all bits set if |in| == 0 (mod p) and 0
815 * in[i] < 2^59 + 2^14
817 static limb felem_is_zero(const felem in)
821 felem_assign(ftmp, in);
823 ftmp[0] += ftmp[8] >> 57;
824 ftmp[8] &= bottom57bits;
826 ftmp[1] += ftmp[0] >> 58;
827 ftmp[0] &= bottom58bits;
828 ftmp[2] += ftmp[1] >> 58;
829 ftmp[1] &= bottom58bits;
830 ftmp[3] += ftmp[2] >> 58;
831 ftmp[2] &= bottom58bits;
832 ftmp[4] += ftmp[3] >> 58;
833 ftmp[3] &= bottom58bits;
834 ftmp[5] += ftmp[4] >> 58;
835 ftmp[4] &= bottom58bits;
836 ftmp[6] += ftmp[5] >> 58;
837 ftmp[5] &= bottom58bits;
838 ftmp[7] += ftmp[6] >> 58;
839 ftmp[6] &= bottom58bits;
840 ftmp[8] += ftmp[7] >> 58;
841 ftmp[7] &= bottom58bits;
842 /* ftmp[8] < 2^57 + 4 */
845 * The ninth limb of 2*(2^521-1) is 0x03ffffffffffffff, which is greater
846 * than our bound for ftmp[8]. Therefore we only have to check if the
847 * zero is zero or 2^521-1.
863 * We know that ftmp[i] < 2^63, therefore the only way that the top bit
864 * can be set is if is_zero was 0 before the decrement.
866 is_zero = 0 - (is_zero >> 63);
868 is_p = ftmp[0] ^ kPrime[0];
869 is_p |= ftmp[1] ^ kPrime[1];
870 is_p |= ftmp[2] ^ kPrime[2];
871 is_p |= ftmp[3] ^ kPrime[3];
872 is_p |= ftmp[4] ^ kPrime[4];
873 is_p |= ftmp[5] ^ kPrime[5];
874 is_p |= ftmp[6] ^ kPrime[6];
875 is_p |= ftmp[7] ^ kPrime[7];
876 is_p |= ftmp[8] ^ kPrime[8];
879 is_p = 0 - (is_p >> 63);
885 static int felem_is_zero_int(const void *in)
887 return (int)(felem_is_zero(in) & ((limb) 1));
891 * felem_contract converts |in| to its unique, minimal representation.
893 * in[i] < 2^59 + 2^14
895 static void felem_contract(felem out, const felem in)
897 limb is_p, is_greater, sign;
898 static const limb two58 = ((limb) 1) << 58;
900 felem_assign(out, in);
902 out[0] += out[8] >> 57;
903 out[8] &= bottom57bits;
905 out[1] += out[0] >> 58;
906 out[0] &= bottom58bits;
907 out[2] += out[1] >> 58;
908 out[1] &= bottom58bits;
909 out[3] += out[2] >> 58;
910 out[2] &= bottom58bits;
911 out[4] += out[3] >> 58;
912 out[3] &= bottom58bits;
913 out[5] += out[4] >> 58;
914 out[4] &= bottom58bits;
915 out[6] += out[5] >> 58;
916 out[5] &= bottom58bits;
917 out[7] += out[6] >> 58;
918 out[6] &= bottom58bits;
919 out[8] += out[7] >> 58;
920 out[7] &= bottom58bits;
921 /* out[8] < 2^57 + 4 */
924 * If the value is greater than 2^521-1 then we have to subtract 2^521-1
925 * out. See the comments in felem_is_zero regarding why we don't test for
926 * other multiples of the prime.
930 * First, if |out| is equal to 2^521-1, we subtract it out to get zero.
933 is_p = out[0] ^ kPrime[0];
934 is_p |= out[1] ^ kPrime[1];
935 is_p |= out[2] ^ kPrime[2];
936 is_p |= out[3] ^ kPrime[3];
937 is_p |= out[4] ^ kPrime[4];
938 is_p |= out[5] ^ kPrime[5];
939 is_p |= out[6] ^ kPrime[6];
940 is_p |= out[7] ^ kPrime[7];
941 is_p |= out[8] ^ kPrime[8];
950 is_p = 0 - (is_p >> 63);
953 /* is_p is 0 iff |out| == 2^521-1 and all ones otherwise */
966 * In order to test that |out| >= 2^521-1 we need only test if out[8] >>
967 * 57 is greater than zero as (2^521-1) + x >= 2^522
969 is_greater = out[8] >> 57;
970 is_greater |= is_greater << 32;
971 is_greater |= is_greater << 16;
972 is_greater |= is_greater << 8;
973 is_greater |= is_greater << 4;
974 is_greater |= is_greater << 2;
975 is_greater |= is_greater << 1;
976 is_greater = 0 - (is_greater >> 63);
978 out[0] -= kPrime[0] & is_greater;
979 out[1] -= kPrime[1] & is_greater;
980 out[2] -= kPrime[2] & is_greater;
981 out[3] -= kPrime[3] & is_greater;
982 out[4] -= kPrime[4] & is_greater;
983 out[5] -= kPrime[5] & is_greater;
984 out[6] -= kPrime[6] & is_greater;
985 out[7] -= kPrime[7] & is_greater;
986 out[8] -= kPrime[8] & is_greater;
988 /* Eliminate negative coefficients */
989 sign = -(out[0] >> 63);
990 out[0] += (two58 & sign);
991 out[1] -= (1 & sign);
992 sign = -(out[1] >> 63);
993 out[1] += (two58 & sign);
994 out[2] -= (1 & sign);
995 sign = -(out[2] >> 63);
996 out[2] += (two58 & sign);
997 out[3] -= (1 & sign);
998 sign = -(out[3] >> 63);
999 out[3] += (two58 & sign);
1000 out[4] -= (1 & sign);
1001 sign = -(out[4] >> 63);
1002 out[4] += (two58 & sign);
1003 out[5] -= (1 & sign);
1004 sign = -(out[0] >> 63);
1005 out[5] += (two58 & sign);
1006 out[6] -= (1 & sign);
1007 sign = -(out[6] >> 63);
1008 out[6] += (two58 & sign);
1009 out[7] -= (1 & sign);
1010 sign = -(out[7] >> 63);
1011 out[7] += (two58 & sign);
1012 out[8] -= (1 & sign);
1013 sign = -(out[5] >> 63);
1014 out[5] += (two58 & sign);
1015 out[6] -= (1 & sign);
1016 sign = -(out[6] >> 63);
1017 out[6] += (two58 & sign);
1018 out[7] -= (1 & sign);
1019 sign = -(out[7] >> 63);
1020 out[7] += (two58 & sign);
1021 out[8] -= (1 & sign);
1028 * Building on top of the field operations we have the operations on the
1029 * elliptic curve group itself. Points on the curve are represented in Jacobian
1033 * point_double calculates 2*(x_in, y_in, z_in)
1035 * The method is taken from:
1036 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
1038 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed.
1039 * while x_out == y_in is not (maybe this works, but it's not tested). */
1041 point_double(felem x_out, felem y_out, felem z_out,
1042 const felem x_in, const felem y_in, const felem z_in)
1044 largefelem tmp, tmp2;
1045 felem delta, gamma, beta, alpha, ftmp, ftmp2;
1047 felem_assign(ftmp, x_in);
1048 felem_assign(ftmp2, x_in);
1051 felem_square(tmp, z_in);
1052 felem_reduce(delta, tmp); /* delta[i] < 2^59 + 2^14 */
1055 felem_square(tmp, y_in);
1056 felem_reduce(gamma, tmp); /* gamma[i] < 2^59 + 2^14 */
1058 /* beta = x*gamma */
1059 felem_mul(tmp, x_in, gamma);
1060 felem_reduce(beta, tmp); /* beta[i] < 2^59 + 2^14 */
1062 /* alpha = 3*(x-delta)*(x+delta) */
1063 felem_diff64(ftmp, delta);
1064 /* ftmp[i] < 2^61 */
1065 felem_sum64(ftmp2, delta);
1066 /* ftmp2[i] < 2^60 + 2^15 */
1067 felem_scalar64(ftmp2, 3);
1068 /* ftmp2[i] < 3*2^60 + 3*2^15 */
1069 felem_mul(tmp, ftmp, ftmp2);
1071 * tmp[i] < 17(3*2^121 + 3*2^76)
1072 * = 61*2^121 + 61*2^76
1073 * < 64*2^121 + 64*2^76
1077 felem_reduce(alpha, tmp);
1079 /* x' = alpha^2 - 8*beta */
1080 felem_square(tmp, alpha);
1082 * tmp[i] < 17*2^120 < 2^125
1084 felem_assign(ftmp, beta);
1085 felem_scalar64(ftmp, 8);
1086 /* ftmp[i] < 2^62 + 2^17 */
1087 felem_diff_128_64(tmp, ftmp);
1088 /* tmp[i] < 2^125 + 2^63 + 2^62 + 2^17 */
1089 felem_reduce(x_out, tmp);
1091 /* z' = (y + z)^2 - gamma - delta */
1092 felem_sum64(delta, gamma);
1093 /* delta[i] < 2^60 + 2^15 */
1094 felem_assign(ftmp, y_in);
1095 felem_sum64(ftmp, z_in);
1096 /* ftmp[i] < 2^60 + 2^15 */
1097 felem_square(tmp, ftmp);
1099 * tmp[i] < 17(2^122) < 2^127
1101 felem_diff_128_64(tmp, delta);
1102 /* tmp[i] < 2^127 + 2^63 */
1103 felem_reduce(z_out, tmp);
1105 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
1106 felem_scalar64(beta, 4);
1107 /* beta[i] < 2^61 + 2^16 */
1108 felem_diff64(beta, x_out);
1109 /* beta[i] < 2^61 + 2^60 + 2^16 */
1110 felem_mul(tmp, alpha, beta);
1112 * tmp[i] < 17*((2^59 + 2^14)(2^61 + 2^60 + 2^16))
1113 * = 17*(2^120 + 2^75 + 2^119 + 2^74 + 2^75 + 2^30)
1114 * = 17*(2^120 + 2^119 + 2^76 + 2^74 + 2^30)
1117 felem_square(tmp2, gamma);
1119 * tmp2[i] < 17*(2^59 + 2^14)^2
1120 * = 17*(2^118 + 2^74 + 2^28)
1122 felem_scalar128(tmp2, 8);
1124 * tmp2[i] < 8*17*(2^118 + 2^74 + 2^28)
1125 * = 2^125 + 2^121 + 2^81 + 2^77 + 2^35 + 2^31
1128 felem_diff128(tmp, tmp2);
1130 * tmp[i] < 2^127 - 2^69 + 17(2^120 + 2^119 + 2^76 + 2^74 + 2^30)
1131 * = 2^127 + 2^124 + 2^122 + 2^120 + 2^118 + 2^80 + 2^78 + 2^76 +
1132 * 2^74 + 2^69 + 2^34 + 2^30
1135 felem_reduce(y_out, tmp);
1138 /* copy_conditional copies in to out iff mask is all ones. */
1139 static void copy_conditional(felem out, const felem in, limb mask)
1142 for (i = 0; i < NLIMBS; ++i) {
1143 const limb tmp = mask & (in[i] ^ out[i]);
1149 * point_add calculates (x1, y1, z1) + (x2, y2, z2)
1151 * The method is taken from
1152 * http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl,
1153 * adapted for mixed addition (z2 = 1, or z2 = 0 for the point at infinity).
1155 * This function includes a branch for checking whether the two input points
1156 * are equal (while not equal to the point at infinity). See comment below
1159 static void point_add(felem x3, felem y3, felem z3,
1160 const felem x1, const felem y1, const felem z1,
1161 const int mixed, const felem x2, const felem y2,
1164 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, ftmp6, x_out, y_out, z_out;
1165 largefelem tmp, tmp2;
1166 limb x_equal, y_equal, z1_is_zero, z2_is_zero;
1169 z1_is_zero = felem_is_zero(z1);
1170 z2_is_zero = felem_is_zero(z2);
1172 /* ftmp = z1z1 = z1**2 */
1173 felem_square(tmp, z1);
1174 felem_reduce(ftmp, tmp);
1177 /* ftmp2 = z2z2 = z2**2 */
1178 felem_square(tmp, z2);
1179 felem_reduce(ftmp2, tmp);
1181 /* u1 = ftmp3 = x1*z2z2 */
1182 felem_mul(tmp, x1, ftmp2);
1183 felem_reduce(ftmp3, tmp);
1185 /* ftmp5 = z1 + z2 */
1186 felem_assign(ftmp5, z1);
1187 felem_sum64(ftmp5, z2);
1188 /* ftmp5[i] < 2^61 */
1190 /* ftmp5 = (z1 + z2)**2 - z1z1 - z2z2 = 2*z1z2 */
1191 felem_square(tmp, ftmp5);
1192 /* tmp[i] < 17*2^122 */
1193 felem_diff_128_64(tmp, ftmp);
1194 /* tmp[i] < 17*2^122 + 2^63 */
1195 felem_diff_128_64(tmp, ftmp2);
1196 /* tmp[i] < 17*2^122 + 2^64 */
1197 felem_reduce(ftmp5, tmp);
1199 /* ftmp2 = z2 * z2z2 */
1200 felem_mul(tmp, ftmp2, z2);
1201 felem_reduce(ftmp2, tmp);
1203 /* s1 = ftmp6 = y1 * z2**3 */
1204 felem_mul(tmp, y1, ftmp2);
1205 felem_reduce(ftmp6, tmp);
1208 * We'll assume z2 = 1 (special case z2 = 0 is handled later)
1211 /* u1 = ftmp3 = x1*z2z2 */
1212 felem_assign(ftmp3, x1);
1214 /* ftmp5 = 2*z1z2 */
1215 felem_scalar(ftmp5, z1, 2);
1217 /* s1 = ftmp6 = y1 * z2**3 */
1218 felem_assign(ftmp6, y1);
1222 felem_mul(tmp, x2, ftmp);
1223 /* tmp[i] < 17*2^120 */
1225 /* h = ftmp4 = u2 - u1 */
1226 felem_diff_128_64(tmp, ftmp3);
1227 /* tmp[i] < 17*2^120 + 2^63 */
1228 felem_reduce(ftmp4, tmp);
1230 x_equal = felem_is_zero(ftmp4);
1232 /* z_out = ftmp5 * h */
1233 felem_mul(tmp, ftmp5, ftmp4);
1234 felem_reduce(z_out, tmp);
1236 /* ftmp = z1 * z1z1 */
1237 felem_mul(tmp, ftmp, z1);
1238 felem_reduce(ftmp, tmp);
1240 /* s2 = tmp = y2 * z1**3 */
1241 felem_mul(tmp, y2, ftmp);
1242 /* tmp[i] < 17*2^120 */
1244 /* r = ftmp5 = (s2 - s1)*2 */
1245 felem_diff_128_64(tmp, ftmp6);
1246 /* tmp[i] < 17*2^120 + 2^63 */
1247 felem_reduce(ftmp5, tmp);
1248 y_equal = felem_is_zero(ftmp5);
1249 felem_scalar64(ftmp5, 2);
1250 /* ftmp5[i] < 2^61 */
1253 * The formulae are incorrect if the points are equal, in affine coordinates
1254 * (X_1, Y_1) == (X_2, Y_2), so we check for this and do doubling if this
1257 * We use bitwise operations to avoid potential side-channels introduced by
1258 * the short-circuiting behaviour of boolean operators.
1260 * The special case of either point being the point at infinity (z1 and/or
1261 * z2 are zero), is handled separately later on in this function, so we
1262 * avoid jumping to point_double here in those special cases.
1264 * Notice the comment below on the implications of this branching for timing
1265 * leaks and why it is considered practically irrelevant.
1267 points_equal = (x_equal & y_equal & (~z1_is_zero) & (~z2_is_zero));
1271 * This is obviously not constant-time but it will almost-never happen
1272 * for ECDH / ECDSA. The case where it can happen is during scalar-mult
1273 * where the intermediate value gets very close to the group order.
1274 * Since |ec_GFp_nistp_recode_scalar_bits| produces signed digits for
1275 * the scalar, it's possible for the intermediate value to be a small
1276 * negative multiple of the base point, and for the final signed digit
1277 * to be the same value. We believe that this only occurs for the scalar
1278 * 1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
1279 * ffffffa51868783bf2f966b7fcc0148f709a5d03bb5c9b8899c47aebb6fb
1280 * 71e913863f7, in that case the penultimate intermediate is -9G and
1281 * the final digit is also -9G. Since this only happens for a single
1282 * scalar, the timing leak is irrelevant. (Any attacker who wanted to
1283 * check whether a secret scalar was that exact value, can already do
1286 point_double(x3, y3, z3, x1, y1, z1);
1290 /* I = ftmp = (2h)**2 */
1291 felem_assign(ftmp, ftmp4);
1292 felem_scalar64(ftmp, 2);
1293 /* ftmp[i] < 2^61 */
1294 felem_square(tmp, ftmp);
1295 /* tmp[i] < 17*2^122 */
1296 felem_reduce(ftmp, tmp);
1298 /* J = ftmp2 = h * I */
1299 felem_mul(tmp, ftmp4, ftmp);
1300 felem_reduce(ftmp2, tmp);
1302 /* V = ftmp4 = U1 * I */
1303 felem_mul(tmp, ftmp3, ftmp);
1304 felem_reduce(ftmp4, tmp);
1306 /* x_out = r**2 - J - 2V */
1307 felem_square(tmp, ftmp5);
1308 /* tmp[i] < 17*2^122 */
1309 felem_diff_128_64(tmp, ftmp2);
1310 /* tmp[i] < 17*2^122 + 2^63 */
1311 felem_assign(ftmp3, ftmp4);
1312 felem_scalar64(ftmp4, 2);
1313 /* ftmp4[i] < 2^61 */
1314 felem_diff_128_64(tmp, ftmp4);
1315 /* tmp[i] < 17*2^122 + 2^64 */
1316 felem_reduce(x_out, tmp);
1318 /* y_out = r(V-x_out) - 2 * s1 * J */
1319 felem_diff64(ftmp3, x_out);
1321 * ftmp3[i] < 2^60 + 2^60 = 2^61
1323 felem_mul(tmp, ftmp5, ftmp3);
1324 /* tmp[i] < 17*2^122 */
1325 felem_mul(tmp2, ftmp6, ftmp2);
1326 /* tmp2[i] < 17*2^120 */
1327 felem_scalar128(tmp2, 2);
1328 /* tmp2[i] < 17*2^121 */
1329 felem_diff128(tmp, tmp2);
1331 * tmp[i] < 2^127 - 2^69 + 17*2^122
1332 * = 2^126 - 2^122 - 2^6 - 2^2 - 1
1335 felem_reduce(y_out, tmp);
1337 copy_conditional(x_out, x2, z1_is_zero);
1338 copy_conditional(x_out, x1, z2_is_zero);
1339 copy_conditional(y_out, y2, z1_is_zero);
1340 copy_conditional(y_out, y1, z2_is_zero);
1341 copy_conditional(z_out, z2, z1_is_zero);
1342 copy_conditional(z_out, z1, z2_is_zero);
1343 felem_assign(x3, x_out);
1344 felem_assign(y3, y_out);
1345 felem_assign(z3, z_out);
1349 * Base point pre computation
1350 * --------------------------
1352 * Two different sorts of precomputed tables are used in the following code.
1353 * Each contain various points on the curve, where each point is three field
1354 * elements (x, y, z).
1356 * For the base point table, z is usually 1 (0 for the point at infinity).
1357 * This table has 16 elements:
1358 * index | bits | point
1359 * ------+---------+------------------------------
1362 * 2 | 0 0 1 0 | 2^130G
1363 * 3 | 0 0 1 1 | (2^130 + 1)G
1364 * 4 | 0 1 0 0 | 2^260G
1365 * 5 | 0 1 0 1 | (2^260 + 1)G
1366 * 6 | 0 1 1 0 | (2^260 + 2^130)G
1367 * 7 | 0 1 1 1 | (2^260 + 2^130 + 1)G
1368 * 8 | 1 0 0 0 | 2^390G
1369 * 9 | 1 0 0 1 | (2^390 + 1)G
1370 * 10 | 1 0 1 0 | (2^390 + 2^130)G
1371 * 11 | 1 0 1 1 | (2^390 + 2^130 + 1)G
1372 * 12 | 1 1 0 0 | (2^390 + 2^260)G
1373 * 13 | 1 1 0 1 | (2^390 + 2^260 + 1)G
1374 * 14 | 1 1 1 0 | (2^390 + 2^260 + 2^130)G
1375 * 15 | 1 1 1 1 | (2^390 + 2^260 + 2^130 + 1)G
1377 * The reason for this is so that we can clock bits into four different
1378 * locations when doing simple scalar multiplies against the base point.
1380 * Tables for other points have table[i] = iG for i in 0 .. 16. */
1382 /* gmul is the table of precomputed base points */
1383 static const felem gmul[16][3] = {
1384 {{0, 0, 0, 0, 0, 0, 0, 0, 0},
1385 {0, 0, 0, 0, 0, 0, 0, 0, 0},
1386 {0, 0, 0, 0, 0, 0, 0, 0, 0}},
1387 {{0x017e7e31c2e5bd66, 0x022cf0615a90a6fe, 0x00127a2ffa8de334,
1388 0x01dfbf9d64a3f877, 0x006b4d3dbaa14b5e, 0x014fed487e0a2bd8,
1389 0x015b4429c6481390, 0x03a73678fb2d988e, 0x00c6858e06b70404},
1390 {0x00be94769fd16650, 0x031c21a89cb09022, 0x039013fad0761353,
1391 0x02657bd099031542, 0x03273e662c97ee72, 0x01e6d11a05ebef45,
1392 0x03d1bd998f544495, 0x03001172297ed0b1, 0x011839296a789a3b},
1393 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1394 {{0x0373faacbc875bae, 0x00f325023721c671, 0x00f666fd3dbde5ad,
1395 0x01a6932363f88ea7, 0x01fc6d9e13f9c47b, 0x03bcbffc2bbf734e,
1396 0x013ee3c3647f3a92, 0x029409fefe75d07d, 0x00ef9199963d85e5},
1397 {0x011173743ad5b178, 0x02499c7c21bf7d46, 0x035beaeabb8b1a58,
1398 0x00f989c4752ea0a3, 0x0101e1de48a9c1a3, 0x01a20076be28ba6c,
1399 0x02f8052e5eb2de95, 0x01bfe8f82dea117c, 0x0160074d3c36ddb7},
1400 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1401 {{0x012f3fc373393b3b, 0x03d3d6172f1419fa, 0x02adc943c0b86873,
1402 0x00d475584177952b, 0x012a4d1673750ee2, 0x00512517a0f13b0c,
1403 0x02b184671a7b1734, 0x0315b84236f1a50a, 0x00a4afc472edbdb9},
1404 {0x00152a7077f385c4, 0x03044007d8d1c2ee, 0x0065829d61d52b52,
1405 0x00494ff6b6631d0d, 0x00a11d94d5f06bcf, 0x02d2f89474d9282e,
1406 0x0241c5727c06eeb9, 0x0386928710fbdb9d, 0x01f883f727b0dfbe},
1407 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1408 {{0x019b0c3c9185544d, 0x006243a37c9d97db, 0x02ee3cbe030a2ad2,
1409 0x00cfdd946bb51e0d, 0x0271c00932606b91, 0x03f817d1ec68c561,
1410 0x03f37009806a369c, 0x03c1f30baf184fd5, 0x01091022d6d2f065},
1411 {0x0292c583514c45ed, 0x0316fca51f9a286c, 0x00300af507c1489a,
1412 0x0295f69008298cf1, 0x02c0ed8274943d7b, 0x016509b9b47a431e,
1413 0x02bc9de9634868ce, 0x005b34929bffcb09, 0x000c1a0121681524},
1414 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1415 {{0x0286abc0292fb9f2, 0x02665eee9805b3f7, 0x01ed7455f17f26d6,
1416 0x0346355b83175d13, 0x006284944cd0a097, 0x0191895bcdec5e51,
1417 0x02e288370afda7d9, 0x03b22312bfefa67a, 0x01d104d3fc0613fe},
1418 {0x0092421a12f7e47f, 0x0077a83fa373c501, 0x03bd25c5f696bd0d,
1419 0x035c41e4d5459761, 0x01ca0d1742b24f53, 0x00aaab27863a509c,
1420 0x018b6de47df73917, 0x025c0b771705cd01, 0x01fd51d566d760a7},
1421 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1422 {{0x01dd92ff6b0d1dbd, 0x039c5e2e8f8afa69, 0x0261ed13242c3b27,
1423 0x0382c6e67026e6a0, 0x01d60b10be2089f9, 0x03c15f3dce86723f,
1424 0x03c764a32d2a062d, 0x017307eac0fad056, 0x018207c0b96c5256},
1425 {0x0196a16d60e13154, 0x03e6ce74c0267030, 0x00ddbf2b4e52a5aa,
1426 0x012738241bbf31c8, 0x00ebe8dc04685a28, 0x024c2ad6d380d4a2,
1427 0x035ee062a6e62d0e, 0x0029ed74af7d3a0f, 0x00eef32aec142ebd},
1428 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1429 {{0x00c31ec398993b39, 0x03a9f45bcda68253, 0x00ac733c24c70890,
1430 0x00872b111401ff01, 0x01d178c23195eafb, 0x03bca2c816b87f74,
1431 0x0261a9af46fbad7a, 0x0324b2a8dd3d28f9, 0x00918121d8f24e23},
1432 {0x032bc8c1ca983cd7, 0x00d869dfb08fc8c6, 0x01693cb61fce1516,
1433 0x012a5ea68f4e88a8, 0x010869cab88d7ae3, 0x009081ad277ceee1,
1434 0x033a77166d064cdc, 0x03955235a1fb3a95, 0x01251a4a9b25b65e},
1435 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1436 {{0x00148a3a1b27f40b, 0x0123186df1b31fdc, 0x00026e7beaad34ce,
1437 0x01db446ac1d3dbba, 0x0299c1a33437eaec, 0x024540610183cbb7,
1438 0x0173bb0e9ce92e46, 0x02b937e43921214b, 0x01ab0436a9bf01b5},
1439 {0x0383381640d46948, 0x008dacbf0e7f330f, 0x03602122bcc3f318,
1440 0x01ee596b200620d6, 0x03bd0585fda430b3, 0x014aed77fd123a83,
1441 0x005ace749e52f742, 0x0390fe041da2b842, 0x0189a8ceb3299242},
1442 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1443 {{0x012a19d6b3282473, 0x00c0915918b423ce, 0x023a954eb94405ae,
1444 0x00529f692be26158, 0x0289fa1b6fa4b2aa, 0x0198ae4ceea346ef,
1445 0x0047d8cdfbdedd49, 0x00cc8c8953f0f6b8, 0x001424abbff49203},
1446 {0x0256732a1115a03a, 0x0351bc38665c6733, 0x03f7b950fb4a6447,
1447 0x000afffa94c22155, 0x025763d0a4dab540, 0x000511e92d4fc283,
1448 0x030a7e9eda0ee96c, 0x004c3cd93a28bf0a, 0x017edb3a8719217f},
1449 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1450 {{0x011de5675a88e673, 0x031d7d0f5e567fbe, 0x0016b2062c970ae5,
1451 0x03f4a2be49d90aa7, 0x03cef0bd13822866, 0x03f0923dcf774a6c,
1452 0x0284bebc4f322f72, 0x016ab2645302bb2c, 0x01793f95dace0e2a},
1453 {0x010646e13527a28f, 0x01ca1babd59dc5e7, 0x01afedfd9a5595df,
1454 0x01f15785212ea6b1, 0x0324e5d64f6ae3f4, 0x02d680f526d00645,
1455 0x0127920fadf627a7, 0x03b383f75df4f684, 0x0089e0057e783b0a},
1456 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1457 {{0x00f334b9eb3c26c6, 0x0298fdaa98568dce, 0x01c2d24843a82292,
1458 0x020bcb24fa1b0711, 0x02cbdb3d2b1875e6, 0x0014907598f89422,
1459 0x03abe3aa43b26664, 0x02cbf47f720bc168, 0x0133b5e73014b79b},
1460 {0x034aab5dab05779d, 0x00cdc5d71fee9abb, 0x0399f16bd4bd9d30,
1461 0x03582fa592d82647, 0x02be1cdfb775b0e9, 0x0034f7cea32e94cb,
1462 0x0335a7f08f56f286, 0x03b707e9565d1c8b, 0x0015c946ea5b614f},
1463 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1464 {{0x024676f6cff72255, 0x00d14625cac96378, 0x00532b6008bc3767,
1465 0x01fc16721b985322, 0x023355ea1b091668, 0x029de7afdc0317c3,
1466 0x02fc8a7ca2da037c, 0x02de1217d74a6f30, 0x013f7173175b73bf},
1467 {0x0344913f441490b5, 0x0200f9e272b61eca, 0x0258a246b1dd55d2,
1468 0x03753db9ea496f36, 0x025e02937a09c5ef, 0x030cbd3d14012692,
1469 0x01793a67e70dc72a, 0x03ec1d37048a662e, 0x006550f700c32a8d},
1470 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1471 {{0x00d3f48a347eba27, 0x008e636649b61bd8, 0x00d3b93716778fb3,
1472 0x004d1915757bd209, 0x019d5311a3da44e0, 0x016d1afcbbe6aade,
1473 0x0241bf5f73265616, 0x0384672e5d50d39b, 0x005009fee522b684},
1474 {0x029b4fab064435fe, 0x018868ee095bbb07, 0x01ea3d6936cc92b8,
1475 0x000608b00f78a2f3, 0x02db911073d1c20f, 0x018205938470100a,
1476 0x01f1e4964cbe6ff2, 0x021a19a29eed4663, 0x01414485f42afa81},
1477 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1478 {{0x01612b3a17f63e34, 0x03813992885428e6, 0x022b3c215b5a9608,
1479 0x029b4057e19f2fcb, 0x0384059a587af7e6, 0x02d6400ace6fe610,
1480 0x029354d896e8e331, 0x00c047ee6dfba65e, 0x0037720542e9d49d},
1481 {0x02ce9eed7c5e9278, 0x0374ed703e79643b, 0x01316c54c4072006,
1482 0x005aaa09054b2ee8, 0x002824000c840d57, 0x03d4eba24771ed86,
1483 0x0189c50aabc3bdae, 0x0338c01541e15510, 0x00466d56e38eed42},
1484 {1, 0, 0, 0, 0, 0, 0, 0, 0}},
1485 {{0x007efd8330ad8bd6, 0x02465ed48047710b, 0x0034c6606b215e0c,
1486 0x016ae30c53cbf839, 0x01fa17bd37161216, 0x018ead4e61ce8ab9,
1487 0x005482ed5f5dee46, 0x037543755bba1d7f, 0x005e5ac7e70a9d0f},
1488 {0x0117e1bb2fdcb2a2, 0x03deea36249f40c4, 0x028d09b4a6246cb7,
1489 0x03524b8855bcf756, 0x023d7d109d5ceb58, 0x0178e43e3223ef9c,
1490 0x0154536a0c6e966a, 0x037964d1286ee9fe, 0x0199bcd90e125055},
1491 {1, 0, 0, 0, 0, 0, 0, 0, 0}}
1495 * select_point selects the |idx|th point from a precomputation table and
1498 /* pre_comp below is of the size provided in |size| */
1499 static void select_point(const limb idx, unsigned int size,
1500 const felem pre_comp[][3], felem out[3])
1503 limb *outlimbs = &out[0][0];
1505 memset(out, 0, sizeof(*out) * 3);
1507 for (i = 0; i < size; i++) {
1508 const limb *inlimbs = &pre_comp[i][0][0];
1509 limb mask = i ^ idx;
1515 for (j = 0; j < NLIMBS * 3; j++)
1516 outlimbs[j] |= inlimbs[j] & mask;
1520 /* get_bit returns the |i|th bit in |in| */
1521 static char get_bit(const felem_bytearray in, int i)
1525 return (in[i >> 3] >> (i & 7)) & 1;
1529 * Interleaved point multiplication using precomputed point multiples: The
1530 * small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[], the scalars
1531 * in scalars[]. If g_scalar is non-NULL, we also add this multiple of the
1532 * generator, using certain (large) precomputed multiples in g_pre_comp.
1533 * Output point (X, Y, Z) is stored in x_out, y_out, z_out
1535 static void batch_mul(felem x_out, felem y_out, felem z_out,
1536 const felem_bytearray scalars[],
1537 const unsigned num_points, const u8 *g_scalar,
1538 const int mixed, const felem pre_comp[][17][3],
1539 const felem g_pre_comp[16][3])
1542 unsigned num, gen_mul = (g_scalar != NULL);
1543 felem nq[3], tmp[4];
1547 /* set nq to the point at infinity */
1548 memset(nq, 0, sizeof(nq));
1551 * Loop over all scalars msb-to-lsb, interleaving additions of multiples
1552 * of the generator (last quarter of rounds) and additions of other
1553 * points multiples (every 5th round).
1555 skip = 1; /* save two point operations in the first
1557 for (i = (num_points ? 520 : 130); i >= 0; --i) {
1560 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1562 /* add multiples of the generator */
1563 if (gen_mul && (i <= 130)) {
1564 bits = get_bit(g_scalar, i + 390) << 3;
1566 bits |= get_bit(g_scalar, i + 260) << 2;
1567 bits |= get_bit(g_scalar, i + 130) << 1;
1568 bits |= get_bit(g_scalar, i);
1570 /* select the point to add, in constant time */
1571 select_point(bits, 16, g_pre_comp, tmp);
1573 /* The 1 argument below is for "mixed" */
1574 point_add(nq[0], nq[1], nq[2],
1575 nq[0], nq[1], nq[2], 1, tmp[0], tmp[1], tmp[2]);
1577 memcpy(nq, tmp, 3 * sizeof(felem));
1582 /* do other additions every 5 doublings */
1583 if (num_points && (i % 5 == 0)) {
1584 /* loop over all scalars */
1585 for (num = 0; num < num_points; ++num) {
1586 bits = get_bit(scalars[num], i + 4) << 5;
1587 bits |= get_bit(scalars[num], i + 3) << 4;
1588 bits |= get_bit(scalars[num], i + 2) << 3;
1589 bits |= get_bit(scalars[num], i + 1) << 2;
1590 bits |= get_bit(scalars[num], i) << 1;
1591 bits |= get_bit(scalars[num], i - 1);
1592 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1595 * select the point to add or subtract, in constant time
1597 select_point(digit, 17, pre_comp[num], tmp);
1598 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative
1600 copy_conditional(tmp[1], tmp[3], (-(limb) sign));
1603 point_add(nq[0], nq[1], nq[2],
1604 nq[0], nq[1], nq[2],
1605 mixed, tmp[0], tmp[1], tmp[2]);
1607 memcpy(nq, tmp, 3 * sizeof(felem));
1613 felem_assign(x_out, nq[0]);
1614 felem_assign(y_out, nq[1]);
1615 felem_assign(z_out, nq[2]);
1618 /* Precomputation for the group generator. */
1619 struct nistp521_pre_comp_st {
1620 felem g_pre_comp[16][3];
1621 CRYPTO_REF_COUNT references;
1622 CRYPTO_RWLOCK *lock;
1625 const EC_METHOD *EC_GFp_nistp521_method(void)
1627 static const EC_METHOD ret = {
1628 EC_FLAGS_DEFAULT_OCT,
1629 NID_X9_62_prime_field,
1630 ec_GFp_nistp521_group_init,
1631 ec_GFp_simple_group_finish,
1632 ec_GFp_simple_group_clear_finish,
1633 ec_GFp_nist_group_copy,
1634 ec_GFp_nistp521_group_set_curve,
1635 ec_GFp_simple_group_get_curve,
1636 ec_GFp_simple_group_get_degree,
1637 ec_group_simple_order_bits,
1638 ec_GFp_simple_group_check_discriminant,
1639 ec_GFp_simple_point_init,
1640 ec_GFp_simple_point_finish,
1641 ec_GFp_simple_point_clear_finish,
1642 ec_GFp_simple_point_copy,
1643 ec_GFp_simple_point_set_to_infinity,
1644 ec_GFp_simple_set_Jprojective_coordinates_GFp,
1645 ec_GFp_simple_get_Jprojective_coordinates_GFp,
1646 ec_GFp_simple_point_set_affine_coordinates,
1647 ec_GFp_nistp521_point_get_affine_coordinates,
1648 0 /* point_set_compressed_coordinates */ ,
1653 ec_GFp_simple_invert,
1654 ec_GFp_simple_is_at_infinity,
1655 ec_GFp_simple_is_on_curve,
1657 ec_GFp_simple_make_affine,
1658 ec_GFp_simple_points_make_affine,
1659 ec_GFp_nistp521_points_mul,
1660 ec_GFp_nistp521_precompute_mult,
1661 ec_GFp_nistp521_have_precompute_mult,
1662 ec_GFp_nist_field_mul,
1663 ec_GFp_nist_field_sqr,
1665 ec_GFp_simple_field_inv,
1666 0 /* field_encode */ ,
1667 0 /* field_decode */ ,
1668 0, /* field_set_to_one */
1669 ec_key_simple_priv2oct,
1670 ec_key_simple_oct2priv,
1671 0, /* set private */
1672 ec_key_simple_generate_key,
1673 ec_key_simple_check_key,
1674 ec_key_simple_generate_public_key,
1677 ecdh_simple_compute_key,
1678 ecdsa_simple_sign_setup,
1679 ecdsa_simple_sign_sig,
1680 ecdsa_simple_verify_sig,
1681 0, /* field_inverse_mod_ord */
1682 0, /* blind_coordinates */
1684 0, /* ladder_step */
1691 /******************************************************************************/
1693 * FUNCTIONS TO MANAGE PRECOMPUTATION
1696 static NISTP521_PRE_COMP *nistp521_pre_comp_new(void)
1698 NISTP521_PRE_COMP *ret = OPENSSL_zalloc(sizeof(*ret));
1701 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1705 ret->references = 1;
1707 ret->lock = CRYPTO_THREAD_lock_new();
1708 if (ret->lock == NULL) {
1709 ECerr(EC_F_NISTP521_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1716 NISTP521_PRE_COMP *EC_nistp521_pre_comp_dup(NISTP521_PRE_COMP *p)
1720 CRYPTO_UP_REF(&p->references, &i, p->lock);
1724 void EC_nistp521_pre_comp_free(NISTP521_PRE_COMP *p)
1731 CRYPTO_DOWN_REF(&p->references, &i, p->lock);
1732 REF_PRINT_COUNT("EC_nistp521", x);
1735 REF_ASSERT_ISNT(i < 0);
1737 CRYPTO_THREAD_lock_free(p->lock);
1741 /******************************************************************************/
1743 * OPENSSL EC_METHOD FUNCTIONS
1746 int ec_GFp_nistp521_group_init(EC_GROUP *group)
1749 ret = ec_GFp_simple_group_init(group);
1750 group->a_is_minus3 = 1;
1754 int ec_GFp_nistp521_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1755 const BIGNUM *a, const BIGNUM *b,
1759 BIGNUM *curve_p, *curve_a, *curve_b;
1761 BN_CTX *new_ctx = NULL;
1764 ctx = new_ctx = BN_CTX_new();
1770 curve_p = BN_CTX_get(ctx);
1771 curve_a = BN_CTX_get(ctx);
1772 curve_b = BN_CTX_get(ctx);
1773 if (curve_b == NULL)
1775 BN_bin2bn(nistp521_curve_params[0], sizeof(felem_bytearray), curve_p);
1776 BN_bin2bn(nistp521_curve_params[1], sizeof(felem_bytearray), curve_a);
1777 BN_bin2bn(nistp521_curve_params[2], sizeof(felem_bytearray), curve_b);
1778 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) || (BN_cmp(curve_b, b))) {
1779 ECerr(EC_F_EC_GFP_NISTP521_GROUP_SET_CURVE,
1780 EC_R_WRONG_CURVE_PARAMETERS);
1783 group->field_mod_func = BN_nist_mod_521;
1784 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1788 BN_CTX_free(new_ctx);
1794 * Takes the Jacobian coordinates (X, Y, Z) of a point and returns (X', Y') =
1797 int ec_GFp_nistp521_point_get_affine_coordinates(const EC_GROUP *group,
1798 const EC_POINT *point,
1799 BIGNUM *x, BIGNUM *y,
1802 felem z1, z2, x_in, y_in, x_out, y_out;
1805 if (EC_POINT_is_at_infinity(group, point)) {
1806 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1807 EC_R_POINT_AT_INFINITY);
1810 if ((!BN_to_felem(x_in, point->X)) || (!BN_to_felem(y_in, point->Y)) ||
1811 (!BN_to_felem(z1, point->Z)))
1814 felem_square(tmp, z2);
1815 felem_reduce(z1, tmp);
1816 felem_mul(tmp, x_in, z1);
1817 felem_reduce(x_in, tmp);
1818 felem_contract(x_out, x_in);
1820 if (!felem_to_BN(x, x_out)) {
1821 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1826 felem_mul(tmp, z1, z2);
1827 felem_reduce(z1, tmp);
1828 felem_mul(tmp, y_in, z1);
1829 felem_reduce(y_in, tmp);
1830 felem_contract(y_out, y_in);
1832 if (!felem_to_BN(y, y_out)) {
1833 ECerr(EC_F_EC_GFP_NISTP521_POINT_GET_AFFINE_COORDINATES,
1841 /* points below is of size |num|, and tmp_felems is of size |num+1/ */
1842 static void make_points_affine(size_t num, felem points[][3],
1846 * Runs in constant time, unless an input is the point at infinity (which
1847 * normally shouldn't happen).
1849 ec_GFp_nistp_points_make_affine_internal(num,
1853 (void (*)(void *))felem_one,
1855 (void (*)(void *, const void *))
1857 (void (*)(void *, const void *))
1858 felem_square_reduce, (void (*)
1865 (void (*)(void *, const void *))
1867 (void (*)(void *, const void *))
1872 * Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL
1873 * values Result is stored in r (r can equal one of the inputs).
1875 int ec_GFp_nistp521_points_mul(const EC_GROUP *group, EC_POINT *r,
1876 const BIGNUM *scalar, size_t num,
1877 const EC_POINT *points[],
1878 const BIGNUM *scalars[], BN_CTX *ctx)
1883 BIGNUM *x, *y, *z, *tmp_scalar;
1884 felem_bytearray g_secret;
1885 felem_bytearray *secrets = NULL;
1886 felem (*pre_comp)[17][3] = NULL;
1887 felem *tmp_felems = NULL;
1890 int have_pre_comp = 0;
1891 size_t num_points = num;
1892 felem x_in, y_in, z_in, x_out, y_out, z_out;
1893 NISTP521_PRE_COMP *pre = NULL;
1894 felem(*g_pre_comp)[3] = NULL;
1895 EC_POINT *generator = NULL;
1896 const EC_POINT *p = NULL;
1897 const BIGNUM *p_scalar = NULL;
1900 x = BN_CTX_get(ctx);
1901 y = BN_CTX_get(ctx);
1902 z = BN_CTX_get(ctx);
1903 tmp_scalar = BN_CTX_get(ctx);
1904 if (tmp_scalar == NULL)
1907 if (scalar != NULL) {
1908 pre = group->pre_comp.nistp521;
1910 /* we have precomputation, try to use it */
1911 g_pre_comp = &pre->g_pre_comp[0];
1913 /* try to use the standard precomputation */
1914 g_pre_comp = (felem(*)[3]) gmul;
1915 generator = EC_POINT_new(group);
1916 if (generator == NULL)
1918 /* get the generator from precomputation */
1919 if (!felem_to_BN(x, g_pre_comp[1][0]) ||
1920 !felem_to_BN(y, g_pre_comp[1][1]) ||
1921 !felem_to_BN(z, g_pre_comp[1][2])) {
1922 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1925 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1929 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1930 /* precomputation matches generator */
1934 * we don't have valid precomputation: treat the generator as a
1940 if (num_points > 0) {
1941 if (num_points >= 2) {
1943 * unless we precompute multiples for just one point, converting
1944 * those into affine form is time well spent
1948 secrets = OPENSSL_zalloc(sizeof(*secrets) * num_points);
1949 pre_comp = OPENSSL_zalloc(sizeof(*pre_comp) * num_points);
1952 OPENSSL_malloc(sizeof(*tmp_felems) * (num_points * 17 + 1));
1953 if ((secrets == NULL) || (pre_comp == NULL)
1954 || (mixed && (tmp_felems == NULL))) {
1955 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1960 * we treat NULL scalars as 0, and NULL points as points at infinity,
1961 * i.e., they contribute nothing to the linear combination
1963 for (i = 0; i < num_points; ++i) {
1966 * we didn't have a valid precomputation, so we pick the
1969 p = EC_GROUP_get0_generator(group);
1972 /* the i^th point */
1974 p_scalar = scalars[i];
1976 if ((p_scalar != NULL) && (p != NULL)) {
1977 /* reduce scalar to 0 <= scalar < 2^521 */
1978 if ((BN_num_bits(p_scalar) > 521)
1979 || (BN_is_negative(p_scalar))) {
1981 * this is an unusual input, and we don't guarantee
1984 if (!BN_nnmod(tmp_scalar, p_scalar, group->order, ctx)) {
1985 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1988 num_bytes = BN_bn2lebinpad(tmp_scalar,
1989 secrets[i], sizeof(secrets[i]));
1991 num_bytes = BN_bn2lebinpad(p_scalar,
1992 secrets[i], sizeof(secrets[i]));
1994 if (num_bytes < 0) {
1995 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
1998 /* precompute multiples */
1999 if ((!BN_to_felem(x_out, p->X)) ||
2000 (!BN_to_felem(y_out, p->Y)) ||
2001 (!BN_to_felem(z_out, p->Z)))
2003 memcpy(pre_comp[i][1][0], x_out, sizeof(felem));
2004 memcpy(pre_comp[i][1][1], y_out, sizeof(felem));
2005 memcpy(pre_comp[i][1][2], z_out, sizeof(felem));
2006 for (j = 2; j <= 16; ++j) {
2008 point_add(pre_comp[i][j][0], pre_comp[i][j][1],
2009 pre_comp[i][j][2], pre_comp[i][1][0],
2010 pre_comp[i][1][1], pre_comp[i][1][2], 0,
2011 pre_comp[i][j - 1][0],
2012 pre_comp[i][j - 1][1],
2013 pre_comp[i][j - 1][2]);
2015 point_double(pre_comp[i][j][0], pre_comp[i][j][1],
2016 pre_comp[i][j][2], pre_comp[i][j / 2][0],
2017 pre_comp[i][j / 2][1],
2018 pre_comp[i][j / 2][2]);
2024 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
2027 /* the scalar for the generator */
2028 if ((scalar != NULL) && (have_pre_comp)) {
2029 memset(g_secret, 0, sizeof(g_secret));
2030 /* reduce scalar to 0 <= scalar < 2^521 */
2031 if ((BN_num_bits(scalar) > 521) || (BN_is_negative(scalar))) {
2033 * this is an unusual input, and we don't guarantee
2036 if (!BN_nnmod(tmp_scalar, scalar, group->order, ctx)) {
2037 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
2040 num_bytes = BN_bn2lebinpad(tmp_scalar, g_secret, sizeof(g_secret));
2042 num_bytes = BN_bn2lebinpad(scalar, g_secret, sizeof(g_secret));
2044 /* do the multiplication with generator precomputation */
2045 batch_mul(x_out, y_out, z_out,
2046 (const felem_bytearray(*))secrets, num_points,
2048 mixed, (const felem(*)[17][3])pre_comp,
2049 (const felem(*)[3])g_pre_comp);
2051 /* do the multiplication without generator precomputation */
2052 batch_mul(x_out, y_out, z_out,
2053 (const felem_bytearray(*))secrets, num_points,
2054 NULL, mixed, (const felem(*)[17][3])pre_comp, NULL);
2056 /* reduce the output to its unique minimal representation */
2057 felem_contract(x_in, x_out);
2058 felem_contract(y_in, y_out);
2059 felem_contract(z_in, z_out);
2060 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
2061 (!felem_to_BN(z, z_in))) {
2062 ECerr(EC_F_EC_GFP_NISTP521_POINTS_MUL, ERR_R_BN_LIB);
2065 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
2069 EC_POINT_free(generator);
2070 OPENSSL_free(secrets);
2071 OPENSSL_free(pre_comp);
2072 OPENSSL_free(tmp_felems);
2076 int ec_GFp_nistp521_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
2079 NISTP521_PRE_COMP *pre = NULL;
2082 EC_POINT *generator = NULL;
2083 felem tmp_felems[16];
2085 BN_CTX *new_ctx = NULL;
2088 /* throw away old precomputation */
2089 EC_pre_comp_free(group);
2093 ctx = new_ctx = BN_CTX_new();
2099 x = BN_CTX_get(ctx);
2100 y = BN_CTX_get(ctx);
2103 /* get the generator */
2104 if (group->generator == NULL)
2106 generator = EC_POINT_new(group);
2107 if (generator == NULL)
2109 BN_bin2bn(nistp521_curve_params[3], sizeof(felem_bytearray), x);
2110 BN_bin2bn(nistp521_curve_params[4], sizeof(felem_bytearray), y);
2111 if (!EC_POINT_set_affine_coordinates(group, generator, x, y, ctx))
2113 if ((pre = nistp521_pre_comp_new()) == NULL)
2116 * if the generator is the standard one, use built-in precomputation
2118 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) {
2119 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
2122 if ((!BN_to_felem(pre->g_pre_comp[1][0], group->generator->X)) ||
2123 (!BN_to_felem(pre->g_pre_comp[1][1], group->generator->Y)) ||
2124 (!BN_to_felem(pre->g_pre_comp[1][2], group->generator->Z)))
2126 /* compute 2^130*G, 2^260*G, 2^390*G */
2127 for (i = 1; i <= 4; i <<= 1) {
2128 point_double(pre->g_pre_comp[2 * i][0], pre->g_pre_comp[2 * i][1],
2129 pre->g_pre_comp[2 * i][2], pre->g_pre_comp[i][0],
2130 pre->g_pre_comp[i][1], pre->g_pre_comp[i][2]);
2131 for (j = 0; j < 129; ++j) {
2132 point_double(pre->g_pre_comp[2 * i][0],
2133 pre->g_pre_comp[2 * i][1],
2134 pre->g_pre_comp[2 * i][2],
2135 pre->g_pre_comp[2 * i][0],
2136 pre->g_pre_comp[2 * i][1],
2137 pre->g_pre_comp[2 * i][2]);
2140 /* g_pre_comp[0] is the point at infinity */
2141 memset(pre->g_pre_comp[0], 0, sizeof(pre->g_pre_comp[0]));
2142 /* the remaining multiples */
2143 /* 2^130*G + 2^260*G */
2144 point_add(pre->g_pre_comp[6][0], pre->g_pre_comp[6][1],
2145 pre->g_pre_comp[6][2], pre->g_pre_comp[4][0],
2146 pre->g_pre_comp[4][1], pre->g_pre_comp[4][2],
2147 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2148 pre->g_pre_comp[2][2]);
2149 /* 2^130*G + 2^390*G */
2150 point_add(pre->g_pre_comp[10][0], pre->g_pre_comp[10][1],
2151 pre->g_pre_comp[10][2], pre->g_pre_comp[8][0],
2152 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
2153 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2154 pre->g_pre_comp[2][2]);
2155 /* 2^260*G + 2^390*G */
2156 point_add(pre->g_pre_comp[12][0], pre->g_pre_comp[12][1],
2157 pre->g_pre_comp[12][2], pre->g_pre_comp[8][0],
2158 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
2159 0, pre->g_pre_comp[4][0], pre->g_pre_comp[4][1],
2160 pre->g_pre_comp[4][2]);
2161 /* 2^130*G + 2^260*G + 2^390*G */
2162 point_add(pre->g_pre_comp[14][0], pre->g_pre_comp[14][1],
2163 pre->g_pre_comp[14][2], pre->g_pre_comp[12][0],
2164 pre->g_pre_comp[12][1], pre->g_pre_comp[12][2],
2165 0, pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
2166 pre->g_pre_comp[2][2]);
2167 for (i = 1; i < 8; ++i) {
2168 /* odd multiples: add G */
2169 point_add(pre->g_pre_comp[2 * i + 1][0],
2170 pre->g_pre_comp[2 * i + 1][1],
2171 pre->g_pre_comp[2 * i + 1][2], pre->g_pre_comp[2 * i][0],
2172 pre->g_pre_comp[2 * i][1], pre->g_pre_comp[2 * i][2], 0,
2173 pre->g_pre_comp[1][0], pre->g_pre_comp[1][1],
2174 pre->g_pre_comp[1][2]);
2176 make_points_affine(15, &(pre->g_pre_comp[1]), tmp_felems);
2179 SETPRECOMP(group, nistp521, pre);
2184 EC_POINT_free(generator);
2186 BN_CTX_free(new_ctx);
2188 EC_nistp521_pre_comp_free(pre);
2192 int ec_GFp_nistp521_have_precompute_mult(const EC_GROUP *group)
2194 return HAVEPRECOMP(group, nistp521);