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