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