1 /* eccdata.c
2
3 Generate compile time constant (but machine dependent) tables.
4
5 Copyright (C) 2013, 2014 Niels Möller
6
7 This file is part of GNU Nettle.
8
9 GNU Nettle is free software: you can redistribute it and/or
10 modify it under the terms of either:
11
12 * the GNU Lesser General Public License as published by the Free
13 Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15
16 or
17
18 * the GNU General Public License as published by the Free
19 Software Foundation; either version 2 of the License, or (at your
20 option) any later version.
21
22 or both in parallel, as here.
23
24 GNU Nettle is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 General Public License for more details.
28
29 You should have received copies of the GNU General Public License and
30 the GNU Lesser General Public License along with this program. If
31 not, see http://www.gnu.org/licenses/.
32 */
33
34 /* Development of Nettle's ECC support was funded by the .SE Internet Fund. */
35
36 #include <assert.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
40
41 #include "mini-gmp.c"
42
43 /* Affine coordinates, for simplicity. Infinity point, i.e., te
44 neutral group element, is represented using the is_zero flag. */
45 struct ecc_point
46 {
47 int is_zero;
48 mpz_t x;
49 mpz_t y;
50 };
51
52 enum ecc_type
53 {
54 /* y^2 = x^3 - 3x + b (mod p) */
55 ECC_TYPE_WEIERSTRASS,
56 /* y^2 = x^3 + b x^2 + x */
57 ECC_TYPE_MONTGOMERY
58 };
59
60 struct ecc_curve
61 {
62 unsigned bit_size;
63 unsigned pippenger_k;
64 unsigned pippenger_c;
65
66 enum ecc_type type;
67
68 /* Prime */
69 mpz_t p;
70 mpz_t b;
71
72 /* Curve order */
73 mpz_t q;
74 struct ecc_point g;
75
76 /* Non-zero if we want elements represented as point s(u, v) on an
77 equivalent Edwards curve, using
78
79 u = t x / y
80 v = (x-1) / (x+1)
81 */
82 int use_edwards;
83 mpz_t d;
84 mpz_t t;
85
86 /* Table for pippenger's algorithm.
87 Element
88
89 i 2^c + j_0 + j_1 2 + j_2 2^2 + ... + j_{c-1} 2^{c-1}
90
91 holds
92
93 2^{ikc} ( j_0 + j_1 2^k + j_2 2^{2k} + ... + j_{c-1} 2^{(c-1)k}) g
94 */
95 mp_size_t table_size;
96 struct ecc_point *table;
97
98 /* If non-NULL, holds 2g, 3g, 4g */
99 struct ecc_point *ref;
100 };
101
102 static void
ecc_init(struct ecc_point * p)103 ecc_init (struct ecc_point *p)
104 {
105 mpz_init (p->x);
106 mpz_init (p->y);
107 }
108
109 static void
ecc_clear(struct ecc_point * p)110 ecc_clear (struct ecc_point *p)
111 {
112 mpz_clear (p->x);
113 mpz_clear (p->y);
114 }
115
116 static int
ecc_zero_p(const struct ecc_point * p)117 ecc_zero_p (const struct ecc_point *p)
118 {
119 return p->is_zero;
120 }
121
122 static int
ecc_equal_p(const struct ecc_point * p,const struct ecc_point * q)123 ecc_equal_p (const struct ecc_point *p, const struct ecc_point *q)
124 {
125 return p->is_zero ? q->is_zero
126 : !q->is_zero && mpz_cmp (p->x, q->x) == 0 && mpz_cmp (p->y, q->y) == 0;
127 }
128
129 static void
ecc_set_zero(struct ecc_point * r)130 ecc_set_zero (struct ecc_point *r)
131 {
132 r->is_zero = 1;
133 }
134
135 static void
ecc_set(struct ecc_point * r,const struct ecc_point * p)136 ecc_set (struct ecc_point *r, const struct ecc_point *p)
137 {
138 r->is_zero = p->is_zero;
139 mpz_set (r->x, p->x);
140 mpz_set (r->y, p->y);
141 }
142
143 /* Needs to support in-place operation. */
144 static void
ecc_dup(const struct ecc_curve * ecc,struct ecc_point * r,const struct ecc_point * p)145 ecc_dup (const struct ecc_curve *ecc,
146 struct ecc_point *r, const struct ecc_point *p)
147 {
148 if (ecc_zero_p (p))
149 ecc_set_zero (r);
150
151 else
152 {
153 mpz_t m, t, x, y;
154
155 mpz_init (m);
156 mpz_init (t);
157 mpz_init (x);
158 mpz_init (y);
159
160 /* m = (2 y)^-1 */
161 mpz_mul_ui (m, p->y, 2);
162 mpz_invert (m, m, ecc->p);
163
164 switch (ecc->type)
165 {
166 case ECC_TYPE_WEIERSTRASS:
167 /* t = 3 (x^2 - 1) * m */
168 mpz_mul (t, p->x, p->x);
169 mpz_mod (t, t, ecc->p);
170 mpz_sub_ui (t, t, 1);
171 mpz_mul_ui (t, t, 3);
172 break;
173 case ECC_TYPE_MONTGOMERY:
174 /* t = (3 x^2 + 2 b x + 1) m = [x(3x+2b)+1] m */
175 mpz_mul_ui (t, ecc->b, 2);
176 mpz_addmul_ui (t, p->x, 3);
177 mpz_mul (t, t, p->x);
178 mpz_mod (t, t, ecc->p);
179 mpz_add_ui (t, t, 1);
180 break;
181 }
182 mpz_mul (t, t, m);
183 mpz_mod (t, t, ecc->p);
184
185 /* x' = t^2 - 2 x */
186 mpz_mul (x, t, t);
187 mpz_submul_ui (x, p->x, 2);
188 if (ecc->type == ECC_TYPE_MONTGOMERY)
189 mpz_sub (x, x, ecc->b);
190
191 mpz_mod (x, x, ecc->p);
192
193 /* y' = (x - x') * t - y */
194 mpz_sub (y, p->x, x);
195 mpz_mul (y, y, t);
196 mpz_sub (y, y, p->y);
197 mpz_mod (y, y, ecc->p);
198
199 r->is_zero = 0;
200 mpz_swap (x, r->x);
201 mpz_swap (y, r->y);
202
203 mpz_clear (m);
204 mpz_clear (t);
205 mpz_clear (x);
206 mpz_clear (y);
207 }
208 }
209
210 static void
ecc_add(const struct ecc_curve * ecc,struct ecc_point * r,const struct ecc_point * p,const struct ecc_point * q)211 ecc_add (const struct ecc_curve *ecc,
212 struct ecc_point *r, const struct ecc_point *p, const struct ecc_point *q)
213 {
214 if (ecc_zero_p (p))
215 ecc_set (r, q);
216
217 else if (ecc_zero_p (q))
218 ecc_set (r, p);
219
220 else if (mpz_cmp (p->x, q->x) == 0)
221 {
222 if (mpz_cmp (p->y, q->y) == 0)
223 ecc_dup (ecc, r, p);
224 else
225 ecc_set_zero (r);
226 }
227 else
228 {
229 mpz_t s, t, x, y;
230 mpz_init (s);
231 mpz_init (t);
232 mpz_init (x);
233 mpz_init (y);
234
235 /* t = (q_y - p_y) / (q_x - p_x) */
236 mpz_sub (t, q->x, p->x);
237 mpz_invert (t, t, ecc->p);
238 mpz_sub (s, q->y, p->y);
239 mpz_mul (t, t, s);
240 mpz_mod (t, t, ecc->p);
241
242 /* x' = t^2 - p_x - q_x */
243 mpz_mul (x, t, t);
244 mpz_sub (x, x, p->x);
245 mpz_sub (x, x, q->x);
246 /* This appears to be the only difference between formulas. */
247 if (ecc->type == ECC_TYPE_MONTGOMERY)
248 mpz_sub (x, x, ecc->b);
249 mpz_mod (x, x, ecc->p);
250
251 /* y' = (x - x') * t - y */
252 mpz_sub (y, p->x, x);
253 mpz_mul (y, y, t);
254 mpz_sub (y, y, p->y);
255 mpz_mod (y, y, ecc->p);
256
257 r->is_zero = 0;
258 mpz_swap (x, r->x);
259 mpz_swap (y, r->y);
260
261 mpz_clear (s);
262 mpz_clear (t);
263 mpz_clear (x);
264 mpz_clear (y);
265 }
266 }
267
268 static void
ecc_mul_binary(const struct ecc_curve * ecc,struct ecc_point * r,const mpz_t n,const struct ecc_point * p)269 ecc_mul_binary (const struct ecc_curve *ecc,
270 struct ecc_point *r, const mpz_t n, const struct ecc_point *p)
271 {
272 /* Avoid the mp_bitcnt_t type for compatibility with older GMP
273 versions. */
274 unsigned k;
275
276 assert (r != p);
277 assert (mpz_sgn (n) > 0);
278
279 ecc_set (r, p);
280
281 /* Index of highest one bit */
282 for (k = mpz_sizeinbase (n, 2) - 1; k-- > 0; )
283 {
284 ecc_dup (ecc, r, r);
285 if (mpz_tstbit (n, k))
286 ecc_add (ecc, r, r, p);
287 }
288 }
289
290 static struct ecc_point *
ecc_alloc(size_t n)291 ecc_alloc (size_t n)
292 {
293 struct ecc_point *p = malloc (n * sizeof(*p));
294 size_t i;
295
296 if (!p)
297 {
298 fprintf (stderr, "Virtual memory exhausted.\n");
299 exit (EXIT_FAILURE);
300 }
301 for (i = 0; i < n; i++)
302 ecc_init (&p[i]);
303
304 return p;
305 }
306
307 static void
ecc_set_str(struct ecc_point * p,const char * x,const char * y)308 ecc_set_str (struct ecc_point *p,
309 const char *x, const char *y)
310 {
311 p->is_zero = 0;
312 mpz_set_str (p->x, x, 16);
313 mpz_set_str (p->y, y, 16);
314 }
315
316 static void
ecc_curve_init_str(struct ecc_curve * ecc,enum ecc_type type,const char * p,const char * b,const char * q,const char * gx,const char * gy,const char * d,const char * t)317 ecc_curve_init_str (struct ecc_curve *ecc, enum ecc_type type,
318 const char *p, const char *b, const char *q,
319 const char *gx, const char *gy,
320 const char *d, const char *t)
321 {
322 ecc->type = type;
323
324 mpz_init_set_str (ecc->p, p, 16);
325 mpz_init_set_str (ecc->b, b, 16);
326 mpz_init_set_str (ecc->q, q, 16);
327 ecc_init (&ecc->g);
328 ecc_set_str (&ecc->g, gx, gy);
329
330 ecc->pippenger_k = 0;
331 ecc->pippenger_c = 0;
332 ecc->table = NULL;
333
334 ecc->ref = NULL;
335
336 mpz_init (ecc->d);
337 mpz_init (ecc->t);
338
339 ecc->use_edwards = (t != NULL);
340 if (ecc->use_edwards)
341 {
342 mpz_set_str (ecc->t, t, 16);
343 mpz_set_str (ecc->d, d, 16);
344 }
345 }
346
347 static void
ecc_curve_init(struct ecc_curve * ecc,unsigned bit_size)348 ecc_curve_init (struct ecc_curve *ecc, unsigned bit_size)
349 {
350 switch (bit_size)
351 {
352 case 192:
353 ecc_curve_init_str (ecc, ECC_TYPE_WEIERSTRASS,
354 /* p = 2^{192} - 2^{64} - 1 */
355 "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE"
356 "FFFFFFFFFFFFFFFF",
357
358 "64210519e59c80e70fa7e9ab72243049"
359 "feb8deecc146b9b1",
360
361 "ffffffffffffffffffffffff99def836"
362 "146bc9b1b4d22831",
363
364 "188da80eb03090f67cbf20eb43a18800"
365 "f4ff0afd82ff1012",
366
367 "07192b95ffc8da78631011ed6b24cdd5"
368 "73f977a11e794811",
369 NULL, NULL);
370 ecc->ref = ecc_alloc (3);
371 ecc_set_str (&ecc->ref[0], /* 2 g */
372 "dafebf5828783f2ad35534631588a3f629a70fb16982a888",
373 "dd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab");
374
375 ecc_set_str (&ecc->ref[1], /* 3 g */
376 "76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da",
377 "782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd");
378
379 ecc_set_str (&ecc->ref[2], /* 4 g */
380 "35433907297cc378b0015703374729d7a4fe46647084e4ba",
381 "a2649984f2135c301ea3acb0776cd4f125389b311db3be32");
382
383 break;
384 case 224:
385 ecc_curve_init_str (ecc, ECC_TYPE_WEIERSTRASS,
386 /* p = 2^{224} - 2^{96} + 1 */
387 "ffffffffffffffffffffffffffffffff"
388 "000000000000000000000001",
389
390 "b4050a850c04b3abf54132565044b0b7"
391 "d7bfd8ba270b39432355ffb4",
392
393 "ffffffffffffffffffffffffffff16a2"
394 "e0b8f03e13dd29455c5c2a3d",
395
396 "b70e0cbd6bb4bf7f321390b94a03c1d3"
397 "56c21122343280d6115c1d21",
398
399 "bd376388b5f723fb4c22dfe6cd4375a0"
400 "5a07476444d5819985007e34",
401 NULL, NULL);
402
403 ecc->ref = ecc_alloc (3);
404 ecc_set_str (&ecc->ref[0], /* 2 g */
405 "706a46dc76dcb76798e60e6d89474788d16dc18032d268fd1a704fa6",
406 "1c2b76a7bc25e7702a704fa986892849fca629487acf3709d2e4e8bb");
407
408 ecc_set_str (&ecc->ref[1], /* 3 g */
409 "df1b1d66a551d0d31eff822558b9d2cc75c2180279fe0d08fd896d04",
410 "a3f7f03cadd0be444c0aa56830130ddf77d317344e1af3591981a925");
411
412 ecc_set_str (&ecc->ref[2], /* 4 g */
413 "ae99feebb5d26945b54892092a8aee02912930fa41cd114e40447301",
414 "482580a0ec5bc47e88bc8c378632cd196cb3fa058a7114eb03054c9");
415
416 break;
417 case 256:
418 ecc_curve_init_str (ecc, ECC_TYPE_WEIERSTRASS,
419 /* p = 2^{256} - 2^{224} + 2^{192} + 2^{96} - 1 */
420 "FFFFFFFF000000010000000000000000"
421 "00000000FFFFFFFFFFFFFFFFFFFFFFFF",
422
423 "5AC635D8AA3A93E7B3EBBD55769886BC"
424 "651D06B0CC53B0F63BCE3C3E27D2604B",
425
426 "FFFFFFFF00000000FFFFFFFFFFFFFFFF"
427 "BCE6FAADA7179E84F3B9CAC2FC632551",
428
429 "6B17D1F2E12C4247F8BCE6E563A440F2"
430 "77037D812DEB33A0F4A13945D898C296",
431
432 "4FE342E2FE1A7F9B8EE7EB4A7C0F9E16"
433 "2BCE33576B315ECECBB6406837BF51F5",
434 NULL, NULL);
435
436 ecc->ref = ecc_alloc (3);
437 ecc_set_str (&ecc->ref[0], /* 2 g */
438 "7cf27b188d034f7e8a52380304b51ac3c08969e277f21b35a60b48fc47669978",
439 "7775510db8ed040293d9ac69f7430dbba7dade63ce982299e04b79d227873d1");
440
441 ecc_set_str (&ecc->ref[1], /* 3 g */
442 "5ecbe4d1a6330a44c8f7ef951d4bf165e6c6b721efada985fb41661bc6e7fd6c",
443 "8734640c4998ff7e374b06ce1a64a2ecd82ab036384fb83d9a79b127a27d5032");
444
445 ecc_set_str (&ecc->ref[2], /* 4 g */
446 "e2534a3532d08fbba02dde659ee62bd0031fe2db785596ef509302446b030852",
447 "e0f1575a4c633cc719dfee5fda862d764efc96c3f30ee0055c42c23f184ed8c6");
448
449 break;
450 case 384:
451 ecc_curve_init_str (ecc, ECC_TYPE_WEIERSTRASS,
452 /* p = 2^{384} - 2^{128} - 2^{96} + 2^{32} - 1 */
453 "ffffffffffffffffffffffffffffffff"
454 "fffffffffffffffffffffffffffffffe"
455 "ffffffff0000000000000000ffffffff",
456
457 "b3312fa7e23ee7e4988e056be3f82d19"
458 "181d9c6efe8141120314088f5013875a"
459 "c656398d8a2ed19d2a85c8edd3ec2aef",
460
461 "ffffffffffffffffffffffffffffffff"
462 "ffffffffffffffffc7634d81f4372ddf"
463 "581a0db248b0a77aecec196accc52973",
464
465 "aa87ca22be8b05378eb1c71ef320ad74"
466 "6e1d3b628ba79b9859f741e082542a38"
467 "5502f25dbf55296c3a545e3872760ab7",
468
469 "3617de4a96262c6f5d9e98bf9292dc29"
470 "f8f41dbd289a147ce9da3113b5f0b8c0"
471 "0a60b1ce1d7e819d7a431d7c90ea0e5f",
472 NULL, NULL);
473
474 ecc->ref = ecc_alloc (3);
475 ecc_set_str (&ecc->ref[0], /* 2 g */
476 "8d999057ba3d2d969260045c55b97f089025959a6f434d651d207d19fb96e9e4fe0e86ebe0e64f85b96a9c75295df61",
477 "8e80f1fa5b1b3cedb7bfe8dffd6dba74b275d875bc6cc43e904e505f256ab4255ffd43e94d39e22d61501e700a940e80");
478
479 ecc_set_str (&ecc->ref[1], /* 3 g */
480 "77a41d4606ffa1464793c7e5fdc7d98cb9d3910202dcd06bea4f240d3566da6b408bbae5026580d02d7e5c70500c831",
481 "c995f7ca0b0c42837d0bbe9602a9fc998520b41c85115aa5f7684c0edc111eacc24abd6be4b5d298b65f28600a2f1df1");
482
483 ecc_set_str (&ecc->ref[2], /* 4 g */
484 "138251cd52ac9298c1c8aad977321deb97e709bd0b4ca0aca55dc8ad51dcfc9d1589a1597e3a5120e1efd631c63e1835",
485 "cacae29869a62e1631e8a28181ab56616dc45d918abc09f3ab0e63cf792aa4dced7387be37bba569549f1c02b270ed67");
486
487 break;
488 case 521:
489 ecc_curve_init_str (ecc, ECC_TYPE_WEIERSTRASS,
490 "1ff" /* p = 2^{521} - 1 */
491 "ffffffffffffffffffffffffffffffff"
492 "ffffffffffffffffffffffffffffffff"
493 "ffffffffffffffffffffffffffffffff"
494 "ffffffffffffffffffffffffffffffff",
495
496 "051"
497 "953eb9618e1c9a1f929a21a0b68540ee"
498 "a2da725b99b315f3b8b489918ef109e1"
499 "56193951ec7e937b1652c0bd3bb1bf07"
500 "3573df883d2c34f1ef451fd46b503f00",
501
502 "1ff"
503 "ffffffffffffffffffffffffffffffff"
504 "fffffffffffffffffffffffffffffffa"
505 "51868783bf2f966b7fcc0148f709a5d0"
506 "3bb5c9b8899c47aebb6fb71e91386409",
507
508 "c6"
509 "858e06b70404e9cd9e3ecb662395b442"
510 "9c648139053fb521f828af606b4d3dba"
511 "a14b5e77efe75928fe1dc127a2ffa8de"
512 "3348b3c1856a429bf97e7e31c2e5bd66",
513
514 "118"
515 "39296a789a3bc0045c8a5fb42c7d1bd9"
516 "98f54449579b446817afbd17273e662c"
517 "97ee72995ef42640c550b9013fad0761"
518 "353c7086a272c24088be94769fd16650",
519 NULL, NULL);
520
521 ecc->ref = ecc_alloc (3);
522 ecc_set_str (&ecc->ref[0], /* 2 g */
523 "433c219024277e7e682fcb288148c282747403279b1ccc06352c6e5505d769be97b3b204da6ef55507aa104a3a35c5af41cf2fa364d60fd967f43e3933ba6d783d",
524 "f4bb8cc7f86db26700a7f3eceeeed3f0b5c6b5107c4da97740ab21a29906c42dbbb3e377de9f251f6b93937fa99a3248f4eafcbe95edc0f4f71be356d661f41b02");
525
526 ecc_set_str (&ecc->ref[1], /* 3 g */
527 "1a73d352443de29195dd91d6a64b5959479b52a6e5b123d9ab9e5ad7a112d7a8dd1ad3f164a3a4832051da6bd16b59fe21baeb490862c32ea05a5919d2ede37ad7d",
528 "13e9b03b97dfa62ddd9979f86c6cab814f2f1557fa82a9d0317d2f8ab1fa355ceec2e2dd4cf8dc575b02d5aced1dec3c70cf105c9bc93a590425f588ca1ee86c0e5");
529
530 ecc_set_str (&ecc->ref[2], /* 4 g */
531 "35b5df64ae2ac204c354b483487c9070cdc61c891c5ff39afc06c5d55541d3ceac8659e24afe3d0750e8b88e9f078af066a1d5025b08e5a5e2fbc87412871902f3",
532 "82096f84261279d2b673e0178eb0b4abb65521aef6e6e32e1b5ae63fe2f19907f279f283e54ba385405224f750a95b85eebb7faef04699d1d9e21f47fc346e4d0d");
533
534 break;
535 case 255:
536 /* curve25519, y^2 = x^3 + 486662 x^2 + x (mod p), with p = 2^{255} - 19.
537
538 According to http://cr.yp.to/papers.html#newelliptic, this
539 is birationally equivalent to the Edwards curve
540
541 x^2 + y^2 = 1 + (121665/121666) x^2 y^2 (mod p).
542
543 And since the constant is not a square, the Edwards formulas
544 should be "complete", with no special cases needed for
545 doubling, neutral element, negatives, etc.
546
547 Generator is x = 9, with y coordinate
548 14781619447589544791020593568409986887264606134616475288964881837755586237401,
549 according to
550
551 x = Mod(9, 2^255-19); sqrt(x^3 + 486662*x^2 + x)
552
553 in PARI/GP. Also, in PARI notation,
554
555 curve25519 = Mod([0, 486662, 0, 1, 0], 2^255-19)
556 */
557 ecc_curve_init_str (ecc, ECC_TYPE_MONTGOMERY,
558 "7fffffffffffffffffffffffffffffff"
559 "ffffffffffffffffffffffffffffffed",
560 "76d06",
561 /* Order of the subgroup is 2^252 + q_0, where
562 q_0 = 27742317777372353535851937790883648493,
563 125 bits.
564 */
565 "10000000000000000000000000000000"
566 "14def9dea2f79cd65812631a5cf5d3ed",
567 "9",
568 /* y coordinate from PARI/GP
569 x = Mod(9, 2^255-19); sqrt(x^3 + 486662*x^2 + x)
570 */
571 "20ae19a1b8a086b4e01edd2c7748d14c"
572 "923d4d7e6d7c61b229e9c5a27eced3d9",
573 /* (121665/121666) mod p, from PARI/GP
574 c = Mod(121665, p); c / (c+1)
575 */
576 "2dfc9311d490018c7338bf8688861767"
577 "ff8ff5b2bebe27548a14b235eca6874a",
578 /* A square root of -486664 mod p, PARI/GP
579 -sqrt(Mod(-486664, p)) in PARI/GP.
580
581 Sign is important to map to the right
582 generator on the twisted edwards curve
583 used for EdDSA. */
584 "70d9120b9f5ff9442d84f723fc03b081"
585 "3a5e2c2eb482e57d3391fb5500ba81e7"
586 );
587 ecc->ref = ecc_alloc (3);
588 ecc_set_str (&ecc->ref[0], /* 2 g */
589 "20d342d51873f1b7d9750c687d157114"
590 "8f3f5ced1e350b5c5cae469cdd684efb",
591 "13b57e011700e8ae050a00945d2ba2f3"
592 "77659eb28d8d391ebcd70465c72df563");
593 ecc_set_str (&ecc->ref[1], /* 3 g */
594 "1c12bc1a6d57abe645534d91c21bba64"
595 "f8824e67621c0859c00a03affb713c12",
596 "2986855cbe387eaeaceea446532c338c"
597 "536af570f71ef7cf75c665019c41222b");
598
599 ecc_set_str (&ecc->ref[2], /* 4 g */
600 "79ce98b7e0689d7de7d1d074a15b315f"
601 "fe1805dfcd5d2a230fee85e4550013ef",
602 "75af5bf4ebdc75c8fe26873427d275d7"
603 "3c0fb13da361077a565539f46de1c30");
604
605 break;
606
607 default:
608 fprintf (stderr, "No known curve for size %d\n", bit_size);
609 exit(EXIT_FAILURE);
610 }
611 ecc->bit_size = bit_size;
612 }
613
614 static void
ecc_curve_clear(struct ecc_curve * ecc)615 ecc_curve_clear (struct ecc_curve *ecc)
616 {
617 mpz_clear (ecc->p);
618 mpz_clear (ecc->b);
619 mpz_clear (ecc->q);
620 ecc_clear (&ecc->g);
621 mpz_clear (ecc->d);
622 mpz_clear (ecc->t);
623 if (ecc->table)
624 {
625 size_t i;
626 for (i = 0; i < ecc->table_size; i++)
627 ecc_clear (&ecc->table[i]);
628 free (ecc->table);
629 }
630 if (ecc->ref)
631 {
632 size_t i;
633 for (i = 0; i < 3; i++)
634 ecc_clear (&ecc->ref[i]);
635 free (ecc->ref);
636 }
637 }
638
639 static unsigned
ecc_table_size(unsigned bits,unsigned k,unsigned c)640 ecc_table_size(unsigned bits, unsigned k, unsigned c)
641 {
642 unsigned p = (bits + k-1) / k;
643 unsigned M = (p + c-1)/c;
644 return M;
645 }
646
647 static void
ecc_pippenger_precompute(struct ecc_curve * ecc,unsigned k,unsigned c)648 ecc_pippenger_precompute (struct ecc_curve *ecc, unsigned k, unsigned c)
649 {
650 unsigned M = ecc_table_size (ecc->bit_size, k, c);
651 unsigned i, j;
652
653 if (M < 2)
654 {
655 fprintf (stderr, "Invalid parameters, implies M = %u\n", M);
656 exit (EXIT_FAILURE);
657 }
658
659 if (M == ecc_table_size (ecc->bit_size, k-1, c))
660 fprintf(stderr,
661 "warn: Parameters k = %u, c = %d are suboptimal, could use smaller k\n",
662 k, c);
663
664 ecc->pippenger_k = k;
665 ecc->pippenger_c = c;
666 ecc->table_size = M << c;
667 assert (ecc->table_size >= 2);
668 ecc->table = ecc_alloc (ecc->table_size);
669
670 /* Compute the first 2^c entries */
671 ecc_set_zero (&ecc->table[0]);
672 ecc_set (&ecc->table[1], &ecc->g);
673
674 for (j = 2; j < (1U<<c); j <<= 1)
675 {
676 /* T[j] = 2^k T[j/2] */
677 assert (j < ecc->table_size);
678 ecc_dup (ecc, &ecc->table[j], &ecc->table[j/2]);
679 for (i = 1; i < k; i++)
680 ecc_dup (ecc, &ecc->table[j], &ecc->table[j]);
681
682 for (i = 1; i < j; i++)
683 {
684 assert (j + i < ecc->table_size);
685 ecc_add (ecc, &ecc->table[j + i], &ecc->table[j], &ecc->table[i]);
686 }
687 }
688 for (j = 1<<c; j < ecc->table_size; j++)
689 {
690 /* T[j] = 2^{kc} T[j-2^c] */
691 ecc_dup (ecc, &ecc->table[j], &ecc->table[j - (1<<c)]);
692 for (i = 1; i < k*c; i++)
693 ecc_dup (ecc, &ecc->table[j], &ecc->table[j]);
694 }
695 }
696
697 static void
ecc_mul_pippenger(const struct ecc_curve * ecc,struct ecc_point * r,const mpz_t n_input)698 ecc_mul_pippenger (const struct ecc_curve *ecc,
699 struct ecc_point *r, const mpz_t n_input)
700 {
701 mpz_t n;
702 unsigned k, c;
703 unsigned i, j;
704 unsigned bit_rows;
705
706 mpz_init (n);
707
708 mpz_mod (n, n_input, ecc->q);
709 ecc_set_zero (r);
710
711 k = ecc->pippenger_k;
712 c = ecc->pippenger_c;
713
714 bit_rows = (ecc->bit_size + k - 1) / k;
715
716 for (i = k; i-- > 0; )
717 {
718 ecc_dup (ecc, r, r);
719 for (j = 0; j * c < bit_rows; j++)
720 {
721 unsigned bits;
722 mp_size_t bit_index;
723
724 /* Extract c bits of the exponent, stride k, starting at i + kcj, ending at
725 i + k (cj + c - 1)*/
726 for (bits = 0, bit_index = i + k*(c*j+c); bit_index > i + k*c*j; )
727 {
728 bit_index -= k;
729 bits = (bits << 1) | mpz_tstbit (n, bit_index);
730 }
731
732 ecc_add (ecc, r, r, &ecc->table[(j << c) | bits]);
733 }
734 }
735 mpz_clear (n);
736 }
737
738 static void
ecc_point_out(FILE * f,const struct ecc_point * p)739 ecc_point_out (FILE *f, const struct ecc_point *p)
740 {
741 if (p->is_zero)
742 fprintf (f, "zero");
743 else
744 {
745 fprintf (f, "(");
746 mpz_out_str (f, 16, p->x);
747 fprintf (f, ",\n ");
748 mpz_out_str (f, 16, (p)->y);
749 fprintf (f, ")");
750 }
751 }
752 #define ASSERT_EQUAL(p, q) do { \
753 if (!ecc_equal_p (p, q)) \
754 { \
755 fprintf (stderr, "%s:%d: ASSERT_EQUAL (%s, %s) failed.\n", \
756 __FILE__, __LINE__, #p, #q); \
757 fprintf (stderr, "p = "); \
758 ecc_point_out (stderr, (p)); \
759 fprintf (stderr, "\nq = "); \
760 ecc_point_out (stderr, (q)); \
761 fprintf (stderr, "\n"); \
762 abort(); \
763 } \
764 } while (0)
765
766 #define ASSERT_ZERO(p) do { \
767 if (!ecc_zero_p (p)) \
768 { \
769 fprintf (stderr, "%s:%d: ASSERT_ZERO (%s) failed.\n", \
770 __FILE__, __LINE__, #p); \
771 fprintf (stderr, "p = "); \
772 ecc_point_out (stderr, (p)); \
773 fprintf (stderr, "\n"); \
774 abort(); \
775 } \
776 } while (0)
777
778 static void
ecc_curve_check(const struct ecc_curve * ecc)779 ecc_curve_check (const struct ecc_curve *ecc)
780 {
781 struct ecc_point p, q;
782 mpz_t n;
783
784 ecc_init (&p);
785 ecc_init (&q);
786 mpz_init (n);
787
788 ecc_dup (ecc, &p, &ecc->g);
789 if (ecc->ref)
790 ASSERT_EQUAL (&p, &ecc->ref[0]);
791 else
792 {
793 fprintf (stderr, "g2 = ");
794 mpz_out_str (stderr, 16, p.x);
795 fprintf (stderr, "\n ");
796 mpz_out_str (stderr, 16, p.y);
797 fprintf (stderr, "\n");
798 }
799 ecc_add (ecc, &q, &p, &ecc->g);
800 if (ecc->ref)
801 ASSERT_EQUAL (&q, &ecc->ref[1]);
802 else
803 {
804 fprintf (stderr, "g3 = ");
805 mpz_out_str (stderr, 16, q.x);
806 fprintf (stderr, "\n ");
807 mpz_out_str (stderr, 16, q.y);
808 fprintf (stderr, "\n");
809 }
810
811 ecc_add (ecc, &q, &q, &ecc->g);
812 if (ecc->ref)
813 ASSERT_EQUAL (&q, &ecc->ref[2]);
814 else
815 {
816 fprintf (stderr, "g4 = ");
817 mpz_out_str (stderr, 16, q.x);
818 fprintf (stderr, "\n ");
819 mpz_out_str (stderr, 16, q.y);
820 fprintf (stderr, "\n");
821 }
822
823 ecc_dup (ecc, &q, &p);
824 if (ecc->ref)
825 ASSERT_EQUAL (&q, &ecc->ref[2]);
826 else
827 {
828 fprintf (stderr, "g4 = ");
829 mpz_out_str (stderr, 16, q.x);
830 fprintf (stderr, "\n ");
831 mpz_out_str (stderr, 16, q.y);
832 fprintf (stderr, "\n");
833 }
834
835 ecc_mul_binary (ecc, &p, ecc->q, &ecc->g);
836 ASSERT_ZERO (&p);
837
838 ecc_mul_pippenger (ecc, &q, ecc->q);
839 ASSERT_ZERO (&q);
840
841 ecc_clear (&p);
842 ecc_clear (&q);
843 mpz_clear (n);
844 }
845
846 static void
output_digits(const mpz_t x,unsigned size,unsigned bits_per_limb)847 output_digits (const mpz_t x,
848 unsigned size, unsigned bits_per_limb)
849 {
850 mpz_t t;
851 mpz_t mask;
852 mpz_t limb;
853 unsigned i;
854 const char *suffix;
855
856 mpz_init (t);
857 mpz_init (mask);
858 mpz_init (limb);
859
860 mpz_setbit (mask, bits_per_limb);
861 mpz_sub_ui (mask, mask, 1);
862
863 suffix = bits_per_limb > 32 ? "ULL" : "UL";
864
865 mpz_init_set (t, x);
866
867 for (i = 0; i < size; i++)
868 {
869 if ( (i % 8) == 0)
870 printf("\n ");
871
872 mpz_and (limb, mask, t);
873 printf (" 0x");
874 mpz_out_str (stdout, 16, limb);
875 printf ("%s,", suffix);
876 mpz_tdiv_q_2exp (t, t, bits_per_limb);
877 }
878
879 mpz_clear (t);
880 mpz_clear (mask);
881 mpz_clear (limb);
882 }
883
884 static void
output_bignum(const char * name,const mpz_t x,unsigned size,unsigned bits_per_limb)885 output_bignum (const char *name, const mpz_t x,
886 unsigned size, unsigned bits_per_limb)
887 {
888 printf ("static const mp_limb_t %s[%d] = {", name, size);
889 output_digits (x, size, bits_per_limb);
890 printf("\n};\n");
891 }
892
893 static void
output_point(const char * name,const struct ecc_curve * ecc,const struct ecc_point * p,int use_redc,unsigned size,unsigned bits_per_limb)894 output_point (const char *name, const struct ecc_curve *ecc,
895 const struct ecc_point *p, int use_redc,
896 unsigned size, unsigned bits_per_limb)
897 {
898 mpz_t x, y, t;
899
900 mpz_init (x);
901 mpz_init (y);
902 mpz_init (t);
903
904 if (name)
905 printf("static const mp_limb_t %s[%u] = {", name, 2*size);
906
907 if (ecc->use_edwards)
908 {
909 if (ecc_zero_p (p))
910 {
911 mpz_set_si (x, 0);
912 mpz_set_si (y, 1);
913 }
914 else if (!mpz_sgn (p->y))
915 {
916 assert (!mpz_sgn (p->x));
917 mpz_set_si (x, 0);
918 mpz_set_si (y, -1);
919 }
920 else
921 {
922 mpz_invert (x, p->y, ecc->p);
923 mpz_mul (x, x, p->x);
924 mpz_mul (x, x, ecc->t);
925 mpz_mod (x, x, ecc->p);
926
927 mpz_sub_ui (y, p->x, 1);
928 mpz_add_ui (t, p->x, 1);
929 mpz_invert (t, t, ecc->p);
930 mpz_mul (y, y, t);
931 mpz_mod (y, y, ecc->p);
932 }
933 }
934 else
935 {
936 mpz_set (x, p->x);
937 mpz_set (y, p->y);
938 }
939 if (use_redc)
940 {
941 mpz_mul_2exp (x, x, size * bits_per_limb);
942 mpz_mod (x, x, ecc->p);
943 mpz_mul_2exp (y, y, size * bits_per_limb);
944 mpz_mod (y, y, ecc->p);
945 }
946
947 output_digits (x, size, bits_per_limb);
948 output_digits (y, size, bits_per_limb);
949
950 if (name)
951 printf("\n};\n");
952
953 mpz_clear (x);
954 mpz_clear (y);
955 mpz_clear (t);
956 }
957
958 static unsigned
output_modulo(const char * name,const mpz_t x,unsigned size,unsigned bits_per_limb)959 output_modulo (const char *name, const mpz_t x,
960 unsigned size, unsigned bits_per_limb)
961 {
962 mpz_t mod;
963 unsigned bits;
964
965 mpz_init (mod);
966
967 mpz_setbit (mod, bits_per_limb * size);
968 mpz_mod (mod, mod, x);
969
970 bits = mpz_sizeinbase (mod, 2);
971 output_bignum (name, mod, size, bits_per_limb);
972
973 mpz_clear (mod);
974 return bits;
975 }
976
977 static void
output_curve(const struct ecc_curve * ecc,unsigned bits_per_limb)978 output_curve (const struct ecc_curve *ecc, unsigned bits_per_limb)
979 {
980 unsigned limb_size = (ecc->bit_size + bits_per_limb - 1)/bits_per_limb;
981 unsigned i;
982 unsigned bits, e;
983 int redc_limbs;
984 mpz_t t;
985
986 mpz_init (t);
987
988 printf ("/* For NULL. */\n#include <stddef.h>\n");
989
990 printf ("#define ECC_LIMB_SIZE %u\n", limb_size);
991 printf ("#define ECC_PIPPENGER_K %u\n", ecc->pippenger_k);
992 printf ("#define ECC_PIPPENGER_C %u\n", ecc->pippenger_c);
993
994 output_bignum ("ecc_p", ecc->p, limb_size, bits_per_limb);
995 output_bignum ("ecc_b", ecc->b, limb_size, bits_per_limb);
996 if (ecc->use_edwards)
997 output_bignum ("ecc_d", ecc->d, limb_size, bits_per_limb);
998 output_bignum ("ecc_q", ecc->q, limb_size, bits_per_limb);
999 output_point ("ecc_g", ecc, &ecc->g, 0, limb_size, bits_per_limb);
1000
1001 bits = output_modulo ("ecc_Bmodp", ecc->p, limb_size, bits_per_limb);
1002 printf ("#define ECC_BMODP_SIZE %u\n",
1003 (bits + bits_per_limb - 1) / bits_per_limb);
1004 bits = output_modulo ("ecc_Bmodq", ecc->q, limb_size, bits_per_limb);
1005 printf ("#define ECC_BMODQ_SIZE %u\n",
1006 (bits + bits_per_limb - 1) / bits_per_limb);
1007 bits = mpz_sizeinbase (ecc->q, 2);
1008 if (bits < ecc->bit_size)
1009 {
1010 /* for curve25519, with q = 2^k + q', with a much smaller q' */
1011 unsigned mbits;
1012 unsigned shift;
1013
1014 /* Shift to align the one bit at B */
1015 shift = bits_per_limb * limb_size + 1 - bits;
1016
1017 mpz_set (t, ecc->q);
1018 mpz_clrbit (t, bits-1);
1019 mbits = mpz_sizeinbase (t, 2);
1020
1021 /* The shifted value must be a limb smaller than q. */
1022 if (mbits + shift + bits_per_limb <= bits)
1023 {
1024 /* q of the form 2^k + q', with q' a limb smaller */
1025 mpz_mul_2exp (t, t, shift);
1026 output_bignum ("ecc_mBmodq_shifted", t, limb_size, bits_per_limb);
1027 }
1028 }
1029
1030 if (ecc->bit_size < limb_size * bits_per_limb)
1031 {
1032 int shift;
1033
1034 mpz_set_ui (t, 0);
1035 mpz_setbit (t, ecc->bit_size);
1036 mpz_sub (t, t, ecc->p);
1037 output_bignum ("ecc_Bmodp_shifted", t, limb_size, bits_per_limb);
1038
1039 shift = limb_size * bits_per_limb - ecc->bit_size;
1040 if (shift > 0)
1041 {
1042 /* Check condition for reducing hi limbs. If s is the
1043 normalization shift and n is the bit size (so that s + n
1044 = limb_size * bite_per_limb), then we need
1045
1046 (2^n - 1) + (2^s - 1) (2^n - p) < 2p
1047
1048 or equivalently,
1049
1050 2^s (2^n - p) <= p
1051
1052 To a allow a carry limb to be added in at the same time,
1053 substitute s+1 for s.
1054 */
1055 /* FIXME: For ecdsa verify, we actually need the stricter
1056 inequality < 2 q. */
1057 mpz_mul_2exp (t, t, shift + 1);
1058 if (mpz_cmp (t, ecc->p) > 0)
1059 {
1060 fprintf (stderr, "Reduction condition failed for %u-bit curve.\n",
1061 ecc->bit_size);
1062 exit (EXIT_FAILURE);
1063 }
1064 }
1065 }
1066 else
1067 printf ("#define ecc_Bmodp_shifted ecc_Bmodp\n");
1068
1069 if (bits < limb_size * bits_per_limb)
1070 {
1071 mpz_set_ui (t, 0);
1072 mpz_setbit (t, bits);
1073 mpz_sub (t, t, ecc->q);
1074 output_bignum ("ecc_Bmodq_shifted", t, limb_size, bits_per_limb);
1075 }
1076 else
1077 printf ("#define ecc_Bmodq_shifted ecc_Bmodq\n");
1078
1079 mpz_add_ui (t, ecc->p, 1);
1080 mpz_fdiv_q_2exp (t, t, 1);
1081 output_bignum ("ecc_pp1h", t, limb_size, bits_per_limb);
1082
1083 mpz_add_ui (t, ecc->q, 1);
1084 mpz_fdiv_q_2exp (t, t, 1);
1085 output_bignum ("ecc_qp1h", t, limb_size, bits_per_limb);
1086
1087 if (ecc->use_edwards)
1088 output_bignum ("ecc_edwards", ecc->t, limb_size, bits_per_limb);
1089
1090 /* Trailing zeros in p+1 correspond to trailing ones in p. */
1091 redc_limbs = mpz_scan0 (ecc->p, 0) / bits_per_limb;
1092 if (redc_limbs > 0)
1093 {
1094 mpz_add_ui (t, ecc->p, 1);
1095 mpz_fdiv_q_2exp (t, t, redc_limbs * bits_per_limb);
1096 output_bignum ("ecc_redc_ppm1", t, limb_size - redc_limbs, bits_per_limb);
1097 }
1098 else
1099 {
1100 /* Trailing zeros in p-1 correspond to zeros just above the low
1101 bit of p */
1102 redc_limbs = mpz_scan1 (ecc->p, 1) / bits_per_limb;
1103 if (redc_limbs > 0)
1104 {
1105 printf ("#define ecc_redc_ppm1 (ecc_p + %d)\n",
1106 redc_limbs);
1107 redc_limbs = -redc_limbs;
1108 }
1109 else
1110 printf ("#define ecc_redc_ppm1 NULL\n");
1111 }
1112 printf ("#define ECC_REDC_SIZE %d\n", redc_limbs);
1113
1114 /* For mod p square root computation. */
1115 if (mpz_fdiv_ui (ecc->p, 4) == 3)
1116 {
1117 /* x = a^{(p+1)/4} gives square root of a (if it exists,
1118 otherwise the square root of -a). */
1119 e = 1;
1120 mpz_add_ui (t, ecc->p, 1);
1121 mpz_fdiv_q_2exp (t, t, 2);
1122 }
1123 else
1124 {
1125 /* p-1 = 2^e s, s odd, t = (s-1)/2*/
1126 unsigned g, i;
1127 mpz_t s;
1128 mpz_t z;
1129
1130 mpz_init (s);
1131 mpz_init (z);
1132
1133 mpz_sub_ui (s, ecc->p, 1);
1134 e = mpz_scan1 (s, 0);
1135 assert (e > 1);
1136
1137 mpz_fdiv_q_2exp (s, s, e);
1138
1139 /* Find a non-square g, g^{(p-1)/2} = -1,
1140 and z = g^{(p-1)/4 */
1141 for (g = 2; ; g++)
1142 {
1143 mpz_set_ui (z, g);
1144 mpz_powm (z, z, s, ecc->p);
1145 mpz_mul (t, z, z);
1146 mpz_mod (t, t, ecc->p);
1147
1148 for (i = 2; i < e; i++)
1149 {
1150 mpz_mul (t, t, t);
1151 mpz_mod (t, t, ecc->p);
1152 }
1153 if (mpz_cmp_ui (t, 1) != 0)
1154 break;
1155 }
1156 mpz_add_ui (t, t, 1);
1157 assert (mpz_cmp (t, ecc->p) == 0);
1158 output_bignum ("ecc_sqrt_z", z, limb_size, bits_per_limb);
1159
1160 mpz_fdiv_q_2exp (t, s, 1);
1161
1162 mpz_clear (s);
1163 mpz_clear (z);
1164 }
1165 printf ("#define ECC_SQRT_E %u\n", e);
1166 printf ("#define ECC_SQRT_T_BITS %u\n",
1167 (unsigned) mpz_sizeinbase (t, 2));
1168 output_bignum ("ecc_sqrt_t", t, limb_size, bits_per_limb);
1169
1170 printf ("#if USE_REDC\n");
1171 printf ("#define ecc_unit ecc_Bmodp\n");
1172
1173 printf ("static const mp_limb_t ecc_table[%lu] = {",
1174 (unsigned long) (2*ecc->table_size * limb_size));
1175 for (i = 0; i < ecc->table_size; i++)
1176 output_point (NULL, ecc, &ecc->table[i], 1, limb_size, bits_per_limb);
1177
1178 printf("\n};\n");
1179
1180 printf ("#else\n");
1181
1182 mpz_set_ui (t, 1);
1183 output_bignum ("ecc_unit", t, limb_size, bits_per_limb);
1184
1185 printf ("static const mp_limb_t ecc_table[%lu] = {",
1186 (unsigned long) (2*ecc->table_size * limb_size));
1187 for (i = 0; i < ecc->table_size; i++)
1188 output_point (NULL, ecc, &ecc->table[i], 0, limb_size, bits_per_limb);
1189
1190 printf("\n};\n");
1191 printf ("#endif\n");
1192
1193 mpz_clear (t);
1194 }
1195
1196 int
main(int argc,char ** argv)1197 main (int argc, char **argv)
1198 {
1199 struct ecc_curve ecc;
1200
1201 if (argc < 4)
1202 {
1203 fprintf (stderr, "Usage: %s CURVE-BITS K C [BITS-PER-LIMB]\n", argv[0]);
1204 return EXIT_FAILURE;
1205 }
1206
1207 ecc_curve_init (&ecc, atoi(argv[1]));
1208
1209 ecc_pippenger_precompute (&ecc, atoi(argv[2]), atoi(argv[3]));
1210
1211 fprintf (stderr, "Table size: %lu entries\n",
1212 (unsigned long) ecc.table_size);
1213
1214 ecc_curve_check (&ecc);
1215
1216 if (argc > 4)
1217 output_curve (&ecc, atoi(argv[4]));
1218
1219 ecc_curve_clear (&ecc);
1220 return EXIT_SUCCESS;
1221 }
1222