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