Revert "Merge branch 'upstream' into tizen"
[platform/upstream/nettle.git] / eccdata.c
1 /* eccdata.c */
2
3 /* Generate compile time constant (but machine dependent) tables. */
4
5 /* nettle, low-level cryptographics library
6  *
7  * Copyright (C) 2013 Niels Möller
8  *  
9  * The nettle library is free software; you can redistribute it and/or modify
10  * it under the terms of the GNU Lesser General Public License as published by
11  * the Free Software Foundation; either version 2.1 of the License, or (at your
12  * option) any later version.
13  * 
14  * The nettle library is distributed in the hope that it will be useful, but
15  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17  * License for more details.
18  * 
19  * You should have received a copy of the GNU Lesser General Public License
20  * along with the nettle library; see the file COPYING.LIB.  If not, write to
21  * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
22  * MA 02111-1301, USA.
23  */
24
25 /* Development of Nettle's ECC support was funded by the .SE Internet Fund. */
26
27 #include <assert.h>
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31
32 #include "mini-gmp.c"
33
34 /* Affine coordinates, for simplicity. Infinity point represented as x
35    == y == 0. */
36 struct ecc_point
37 {
38   mpz_t x;
39   mpz_t y;
40 };
41
42 /* Represents an elliptic curve of the form
43
44      y^2 = x^3 - 3x + b (mod p)
45 */
46 struct ecc_curve
47 {
48   unsigned bit_size;
49   unsigned pippenger_k;
50   unsigned pippenger_c;
51
52   /* Prime */
53   mpz_t p;
54   mpz_t b;
55
56   /* Curve order */
57   mpz_t q;
58   struct ecc_point g;
59
60   /* Table for pippenger's algorithm.
61      Element
62
63        i 2^c + j_0 + j_1 2 + j_2 2^2 + ... + j_{c-1} 2^{c-1}
64
65      holds
66
67        2^{ikc} ( j_0 + j_1 2^k + j_2 2^{2k} + ... + j_{c-1} 2^{(c-1)k}) g
68    */
69   mp_size_t table_size;
70   struct ecc_point *table;
71
72   /* If non-NULL, holds 2g, 3g, 4g */
73   struct ecc_point *ref;
74 };
75
76 static void
77 ecc_init (struct ecc_point *p)
78 {
79   mpz_init (p->x);
80   mpz_init (p->y);
81 }
82
83 static void
84 ecc_clear (struct ecc_point *p)
85 {
86   mpz_clear (p->x);
87   mpz_clear (p->y);
88 }
89
90 static int
91 ecc_zero_p (const struct ecc_point *p)
92 {
93   return mpz_sgn (p->x) == 0 && mpz_sgn (p->y) == 0;
94 }
95
96 static int
97 ecc_equal_p (const struct ecc_point *p, const struct ecc_point *q)
98 {
99   return mpz_cmp (p->x, q->x) == 0 && mpz_cmp (p->y, q->y) == 0;
100 }
101
102 static void
103 ecc_set_zero (struct ecc_point *r)
104 {
105   mpz_set_ui (r->x, 0);
106   mpz_set_ui (r->y, 0);
107 }
108
109 static void
110 ecc_set (struct ecc_point *r, const struct ecc_point *p)
111 {
112   mpz_set (r->x, p->x);
113   mpz_set (r->y, p->y);
114 }
115
116 static void
117 ecc_dup (const struct ecc_curve *ecc,
118          struct ecc_point *r, const struct ecc_point *p)
119 {
120   if (ecc_zero_p (p))
121     ecc_set_zero (r);
122
123   else
124     {
125       mpz_t m, t, x, y;
126   
127       mpz_init (m);
128       mpz_init (t);
129       mpz_init (x);
130       mpz_init (y);
131
132       /* m = (2 y)^-1 */
133       mpz_mul_ui (m, p->y, 2);
134       mpz_invert (m, m, ecc->p);
135
136       /* t = 3 (x^2 - 1) * m */
137       mpz_mul (t, p->x, p->x);
138       mpz_mod (t, t, ecc->p);
139       mpz_sub_ui (t, t, 1);
140       mpz_mul_ui (t, t, 3);
141       mpz_mul (t, t, m);
142
143       /* x' = t^2 - 2 x */
144       mpz_mul (x, t, t);
145       /* mpz_submul_ui (x, p->x, 2); not available in mini-gmp */
146       mpz_mul_ui (m, p->x, 2);
147       mpz_sub (x, x, m);
148       mpz_mod (x, x, ecc->p);
149
150       /* y' = (x - x') * t - y */
151       mpz_sub (y, p->x, x);
152       mpz_mul (y, y, t);
153       mpz_sub (y, y, p->y);
154       mpz_mod (y, y, ecc->p);
155
156       mpz_swap (x, r->x);
157       mpz_swap (y, r->y);
158   
159       mpz_clear (m);
160       mpz_clear (t);
161       mpz_clear (x);
162       mpz_clear (y);
163     }
164 }
165
166 static void
167 ecc_add (const struct ecc_curve *ecc,
168          struct ecc_point *r, const struct ecc_point *p, const struct ecc_point *q)
169 {
170   if (ecc_zero_p (p))
171     ecc_set (r, q);
172
173   else if (ecc_zero_p (q))
174     ecc_set (r, p);
175
176   else if (mpz_cmp (p->x, q->x) == 0)
177     {
178       if (mpz_cmp (p->y, q->y) == 0)
179         ecc_dup (ecc, r, p);
180       else
181         ecc_set_zero (r);
182     }
183   else
184     {
185       mpz_t s, t, x, y;
186       mpz_init (s);
187       mpz_init (t);
188       mpz_init (x);
189       mpz_init (y);
190
191       /* t = (q_y - p_y) / (q_x - p_x) */
192       mpz_sub (t, q->x, p->x);
193       mpz_invert (t, t, ecc->p);
194       mpz_sub (s, q->y, p->y);
195       mpz_mul (t, t, s);
196       mpz_mod (t, t, ecc->p);
197
198       /* x' = t^2 - p_x - q_x */
199       mpz_mul (x, t, t);
200       mpz_sub (x, x, p->x);
201       mpz_sub (x, x, q->x);
202       mpz_mod (x, x, ecc->p);
203
204       /* y' = (x - x') * t - y */
205       mpz_sub (y, p->x, x);
206       mpz_mul (y, y, t);
207       mpz_sub (y, y, p->y);
208       mpz_mod (y, y, ecc->p);
209
210       mpz_swap (x, r->x);
211       mpz_swap (y, r->y);
212
213       mpz_clear (s);
214       mpz_clear (t);
215       mpz_clear (x);
216       mpz_clear (y);
217     }
218 }
219
220 static void 
221 ecc_mul_binary (const struct ecc_curve *ecc,
222                 struct ecc_point *r, const mpz_t n, const struct ecc_point *p)
223 {
224   /* Avoid the mp_bitcnt_t type for compatibility with older GMP
225      versions. */
226   unsigned k;
227
228   assert (r != p);
229   assert (mpz_sgn (n) > 0);
230
231   ecc_set (r, p);
232
233   /* Index of highest one bit */
234   for (k = mpz_sizeinbase (n, 2) - 1; k-- > 0; )
235     {
236       ecc_dup (ecc, r, r);
237       if (mpz_tstbit (n, k))
238         ecc_add (ecc, r, r, p);
239     }  
240 }
241
242 static struct ecc_point *
243 ecc_alloc (size_t n)
244 {
245   struct ecc_point *p = malloc (n * sizeof(*p));
246   size_t i;
247
248   if (!p)
249     {
250       fprintf (stderr, "Virtual memory exhausted.\n");
251       exit (EXIT_FAILURE);
252     }
253   for (i = 0; i < n; i++)
254     ecc_init (&p[i]);
255
256   return p;
257 }
258
259 static void
260 ecc_set_str (struct ecc_point *p,
261              const char *x, const char *y)
262 {
263   mpz_set_str (p->x, x, 16);
264   mpz_set_str (p->y, y, 16);  
265 }
266
267 static void
268 ecc_curve_init_str (struct ecc_curve *ecc,
269                     const char *p, const char *b, const char *q,
270                     const char *gx, const char *gy)
271 {
272   mpz_init_set_str (ecc->p, p, 16);
273   mpz_init_set_str (ecc->b, b, 16);
274   mpz_init_set_str (ecc->q, q, 16);
275   ecc_init (&ecc->g);
276   ecc_set_str (&ecc->g, gx, gy);
277
278   ecc->pippenger_k = 0;
279   ecc->pippenger_c = 0;
280   ecc->table = NULL;
281
282   ecc->ref = NULL;
283 }
284
285 static void
286 ecc_curve_init (struct ecc_curve *ecc, unsigned bit_size)
287 {
288   switch (bit_size)
289     {
290     case 192:      
291       ecc_curve_init_str (ecc,
292                           /* p = 2^{192} - 2^{64} - 1 */
293                           "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE"
294                           "FFFFFFFFFFFFFFFF",
295
296                           "64210519e59c80e70fa7e9ab72243049"
297                           "feb8deecc146b9b1", 
298
299                           "ffffffffffffffffffffffff99def836"
300                           "146bc9b1b4d22831",
301
302                           "188da80eb03090f67cbf20eb43a18800"
303                           "f4ff0afd82ff1012",
304
305                           "07192b95ffc8da78631011ed6b24cdd5"
306                           "73f977a11e794811");
307       ecc->ref = ecc_alloc (3);
308       ecc_set_str (&ecc->ref[0], /* 2 g */
309                    "dafebf5828783f2ad35534631588a3f629a70fb16982a888",
310                    "dd6bda0d993da0fa46b27bbc141b868f59331afa5c7e93ab");
311       
312       ecc_set_str (&ecc->ref[1], /* 3 g */
313                    "76e32a2557599e6edcd283201fb2b9aadfd0d359cbb263da",
314                    "782c37e372ba4520aa62e0fed121d49ef3b543660cfd05fd");
315
316       ecc_set_str (&ecc->ref[2], /* 4 g */
317                    "35433907297cc378b0015703374729d7a4fe46647084e4ba",
318                    "a2649984f2135c301ea3acb0776cd4f125389b311db3be32");
319
320       break;
321     case 224:
322       ecc_curve_init_str (ecc,
323                           /* p = 2^{224} - 2^{96} + 1 */
324                           "ffffffffffffffffffffffffffffffff"
325                           "000000000000000000000001",
326
327                           "b4050a850c04b3abf54132565044b0b7"
328                           "d7bfd8ba270b39432355ffb4",
329
330                           "ffffffffffffffffffffffffffff16a2"
331                           "e0b8f03e13dd29455c5c2a3d",
332
333                           "b70e0cbd6bb4bf7f321390b94a03c1d3"
334                           "56c21122343280d6115c1d21",
335
336                           "bd376388b5f723fb4c22dfe6cd4375a0"
337                           "5a07476444d5819985007e34");
338
339       ecc->ref = ecc_alloc (3);
340       ecc_set_str (&ecc->ref[0], /* 2 g */
341                    "706a46dc76dcb76798e60e6d89474788d16dc18032d268fd1a704fa6",
342                    "1c2b76a7bc25e7702a704fa986892849fca629487acf3709d2e4e8bb");
343       
344       ecc_set_str (&ecc->ref[1], /* 3 g */
345                    "df1b1d66a551d0d31eff822558b9d2cc75c2180279fe0d08fd896d04",
346                    "a3f7f03cadd0be444c0aa56830130ddf77d317344e1af3591981a925");
347
348       ecc_set_str (&ecc->ref[2], /* 4 g */
349                    "ae99feebb5d26945b54892092a8aee02912930fa41cd114e40447301",
350                    "482580a0ec5bc47e88bc8c378632cd196cb3fa058a7114eb03054c9");
351
352       break;
353     case 256:
354       ecc_curve_init_str (ecc,
355                           /* p = 2^{256} - 2^{224} + 2^{192} + 2^{96} - 1 */
356                           "FFFFFFFF000000010000000000000000"
357                           "00000000FFFFFFFFFFFFFFFFFFFFFFFF",
358
359                           "5AC635D8AA3A93E7B3EBBD55769886BC"
360                           "651D06B0CC53B0F63BCE3C3E27D2604B",
361
362                           "FFFFFFFF00000000FFFFFFFFFFFFFFFF"
363                           "BCE6FAADA7179E84F3B9CAC2FC632551",
364
365                           "6B17D1F2E12C4247F8BCE6E563A440F2"
366                           "77037D812DEB33A0F4A13945D898C296",
367
368                           "4FE342E2FE1A7F9B8EE7EB4A7C0F9E16"
369                           "2BCE33576B315ECECBB6406837BF51F5");
370
371       ecc->ref = ecc_alloc (3);
372       ecc_set_str (&ecc->ref[0], /* 2 g */
373                    "7cf27b188d034f7e8a52380304b51ac3c08969e277f21b35a60b48fc47669978",
374                    "7775510db8ed040293d9ac69f7430dbba7dade63ce982299e04b79d227873d1");
375       
376       ecc_set_str (&ecc->ref[1], /* 3 g */
377                    "5ecbe4d1a6330a44c8f7ef951d4bf165e6c6b721efada985fb41661bc6e7fd6c",
378                    "8734640c4998ff7e374b06ce1a64a2ecd82ab036384fb83d9a79b127a27d5032");
379
380       ecc_set_str (&ecc->ref[2], /* 4 g */
381                    "e2534a3532d08fbba02dde659ee62bd0031fe2db785596ef509302446b030852",
382                    "e0f1575a4c633cc719dfee5fda862d764efc96c3f30ee0055c42c23f184ed8c6");
383
384       break;
385     case 384:
386       ecc_curve_init_str (ecc,
387                           /* p = 2^{384} - 2^{128} - 2^{96} + 2^{32} - 1 */
388                           "ffffffffffffffffffffffffffffffff"
389                           "fffffffffffffffffffffffffffffffe"
390                           "ffffffff0000000000000000ffffffff",
391                           
392                           "b3312fa7e23ee7e4988e056be3f82d19"
393                           "181d9c6efe8141120314088f5013875a"
394                           "c656398d8a2ed19d2a85c8edd3ec2aef",
395                           
396                           "ffffffffffffffffffffffffffffffff"
397                           "ffffffffffffffffc7634d81f4372ddf"
398                           "581a0db248b0a77aecec196accc52973",
399                           
400                           "aa87ca22be8b05378eb1c71ef320ad74"
401                           "6e1d3b628ba79b9859f741e082542a38"
402                           "5502f25dbf55296c3a545e3872760ab7",
403                           
404                           "3617de4a96262c6f5d9e98bf9292dc29"
405                           "f8f41dbd289a147ce9da3113b5f0b8c0"
406                           "0a60b1ce1d7e819d7a431d7c90ea0e5f");
407
408       ecc->ref = ecc_alloc (3);
409       ecc_set_str (&ecc->ref[0], /* 2 g */
410                    "8d999057ba3d2d969260045c55b97f089025959a6f434d651d207d19fb96e9e4fe0e86ebe0e64f85b96a9c75295df61",
411                    "8e80f1fa5b1b3cedb7bfe8dffd6dba74b275d875bc6cc43e904e505f256ab4255ffd43e94d39e22d61501e700a940e80");
412
413       ecc_set_str (&ecc->ref[1], /* 3 g */
414                    "77a41d4606ffa1464793c7e5fdc7d98cb9d3910202dcd06bea4f240d3566da6b408bbae5026580d02d7e5c70500c831",
415                    "c995f7ca0b0c42837d0bbe9602a9fc998520b41c85115aa5f7684c0edc111eacc24abd6be4b5d298b65f28600a2f1df1");
416
417       ecc_set_str (&ecc->ref[2], /* 4 g */
418                    "138251cd52ac9298c1c8aad977321deb97e709bd0b4ca0aca55dc8ad51dcfc9d1589a1597e3a5120e1efd631c63e1835",
419                    "cacae29869a62e1631e8a28181ab56616dc45d918abc09f3ab0e63cf792aa4dced7387be37bba569549f1c02b270ed67");
420
421       break;
422     case 521:
423       ecc_curve_init_str (ecc,                    
424                           "1ff" /* p = 2^{521} - 1 */
425                           "ffffffffffffffffffffffffffffffff"
426                           "ffffffffffffffffffffffffffffffff"
427                           "ffffffffffffffffffffffffffffffff"
428                           "ffffffffffffffffffffffffffffffff",
429
430                           "051"
431                           "953eb9618e1c9a1f929a21a0b68540ee"
432                           "a2da725b99b315f3b8b489918ef109e1"
433                           "56193951ec7e937b1652c0bd3bb1bf07"
434                           "3573df883d2c34f1ef451fd46b503f00",
435
436                           "1ff"
437                           "ffffffffffffffffffffffffffffffff"
438                           "fffffffffffffffffffffffffffffffa"
439                           "51868783bf2f966b7fcc0148f709a5d0"
440                           "3bb5c9b8899c47aebb6fb71e91386409",
441
442                           "c6"
443                           "858e06b70404e9cd9e3ecb662395b442"
444                           "9c648139053fb521f828af606b4d3dba"
445                           "a14b5e77efe75928fe1dc127a2ffa8de"
446                           "3348b3c1856a429bf97e7e31c2e5bd66",
447
448                           "118"
449                           "39296a789a3bc0045c8a5fb42c7d1bd9"
450                           "98f54449579b446817afbd17273e662c"
451                           "97ee72995ef42640c550b9013fad0761"
452                           "353c7086a272c24088be94769fd16650");
453
454       ecc->ref = ecc_alloc (3);
455       ecc_set_str (&ecc->ref[0], /* 2 g */
456                    "433c219024277e7e682fcb288148c282747403279b1ccc06352c6e5505d769be97b3b204da6ef55507aa104a3a35c5af41cf2fa364d60fd967f43e3933ba6d783d",
457                    "f4bb8cc7f86db26700a7f3eceeeed3f0b5c6b5107c4da97740ab21a29906c42dbbb3e377de9f251f6b93937fa99a3248f4eafcbe95edc0f4f71be356d661f41b02");
458       
459       ecc_set_str (&ecc->ref[1], /* 3 g */
460                    "1a73d352443de29195dd91d6a64b5959479b52a6e5b123d9ab9e5ad7a112d7a8dd1ad3f164a3a4832051da6bd16b59fe21baeb490862c32ea05a5919d2ede37ad7d",
461                    "13e9b03b97dfa62ddd9979f86c6cab814f2f1557fa82a9d0317d2f8ab1fa355ceec2e2dd4cf8dc575b02d5aced1dec3c70cf105c9bc93a590425f588ca1ee86c0e5");
462
463       ecc_set_str (&ecc->ref[2], /* 4 g */
464                    "35b5df64ae2ac204c354b483487c9070cdc61c891c5ff39afc06c5d55541d3ceac8659e24afe3d0750e8b88e9f078af066a1d5025b08e5a5e2fbc87412871902f3",
465                    "82096f84261279d2b673e0178eb0b4abb65521aef6e6e32e1b5ae63fe2f19907f279f283e54ba385405224f750a95b85eebb7faef04699d1d9e21f47fc346e4d0d");
466
467       break;
468     default:
469       fprintf (stderr, "No known curve for size %d\n", bit_size);
470       exit(EXIT_FAILURE);     
471     }
472   ecc->bit_size = bit_size;
473 }
474
475 static void
476 ecc_pippenger_precompute (struct ecc_curve *ecc, unsigned k, unsigned c)
477 {
478   unsigned p = (ecc->bit_size + k-1) / k;
479   unsigned M = (p + c-1)/c;
480   unsigned i, j;
481
482   ecc->pippenger_k = k;
483   ecc->pippenger_c = c;
484   ecc->table_size = M << c;
485   ecc->table = ecc_alloc (ecc->table_size);
486   
487   /* Compute the first 2^c entries */
488   ecc_set_zero (&ecc->table[0]);
489   ecc_set (&ecc->table[1], &ecc->g);
490
491   for (j = 2; j < (1U<<c); j <<= 1)
492     {
493       /* T[j] = 2^k T[j/2] */
494       ecc_dup (ecc, &ecc->table[j], &ecc->table[j/2]);
495       for (i = 1; i < k; i++)
496         ecc_dup (ecc, &ecc->table[j], &ecc->table[j]);
497
498       for (i = 1; i < j; i++)
499         ecc_add (ecc, &ecc->table[j + i], &ecc->table[j], &ecc->table[i]);
500     }
501   for (j = 1<<c; j < ecc->table_size; j++)
502     {
503       /* T[j] = 2^{kc} T[j-2^c] */
504       ecc_dup (ecc, &ecc->table[j], &ecc->table[j - (1<<c)]);
505       for (i = 1; i < k*c; i++)
506         ecc_dup (ecc, &ecc->table[j], &ecc->table[j]);
507     }
508 }
509
510 static void
511 ecc_mul_pippenger (const struct ecc_curve *ecc,
512                    struct ecc_point *r, const mpz_t n_input)
513 {
514   mpz_t n;
515   unsigned k, c;
516   unsigned i, j;
517   unsigned bit_rows;
518
519   mpz_init (n);
520   
521   mpz_mod (n, n_input, ecc->q);
522   ecc_set_zero (r);
523
524   k = ecc->pippenger_k;
525   c = ecc->pippenger_c;
526
527   bit_rows = (ecc->bit_size + k - 1) / k;
528
529   for (i = k; i-- > 0; )
530     {
531       ecc_dup (ecc, r, r);
532       for (j = 0; j * c < bit_rows; j++)
533         {
534           unsigned bits;
535           mp_size_t bit_index;
536           
537           /* Extract c bits of the exponent, stride k, starting at i + kcj, ending at
538             i + k (cj + c - 1)*/
539           for (bits = 0, bit_index = i + k*(c*j+c); bit_index > i + k*c*j; )
540             {
541               bit_index -= k;
542               bits = (bits << 1) | mpz_tstbit (n, bit_index);
543             }
544
545           ecc_add (ecc, r, r, &ecc->table[(j << c) | bits]);
546         }
547     }
548   mpz_clear (n);
549 }
550
551 #define ASSERT_EQUAL(p, q) do {                                         \
552     if (!ecc_equal_p (p, q))                                            \
553       {                                                                 \
554         fprintf (stderr, "%s:%d: ASSERT_EQUAL (%s, %s) failed.\n",      \
555                  __FILE__, __LINE__, #p, #q);                           \
556         fprintf (stderr, "p = (");                                      \
557         mpz_out_str (stderr, 16, (p)->x);                               \
558         fprintf (stderr, ",\n     ");                                   \
559         mpz_out_str (stderr, 16, (p)->y);                               \
560         fprintf (stderr, ")\nq = (");                                   \
561         mpz_out_str (stderr, 16, (q)->x);                               \
562         fprintf (stderr, ",\n     ");                                   \
563         mpz_out_str (stderr, 16, (q)->y);                               \
564         fprintf (stderr, ")\n");                                        \
565         abort();                                                        \
566       }                                                                 \
567   } while (0)
568
569 #define ASSERT_ZERO(p) do {                                             \
570     if (!ecc_zero_p (p))                                                \
571       {                                                                 \
572         fprintf (stderr, "%s:%d: ASSERT_ZERO (%s) failed.\n",           \
573                  __FILE__, __LINE__, #p);                               \
574         fprintf (stderr, "p = (");                                      \
575         mpz_out_str (stderr, 16, (p)->x);                               \
576         fprintf (stderr, ",\n     ");                                   \
577         mpz_out_str (stderr, 16, (p)->y);                               \
578         fprintf (stderr, ")\n");                                        \
579         abort();                                                        \
580       }                                                                 \
581   } while (0)
582
583 static void
584 ecc_curve_check (const struct ecc_curve *ecc)
585 {
586   struct ecc_point p, q;
587   mpz_t n;
588
589   ecc_init (&p);
590   ecc_init (&q);
591   mpz_init (n);
592
593   ecc_dup (ecc, &p, &ecc->g);
594   if (ecc->ref)
595     ASSERT_EQUAL (&p, &ecc->ref[0]);
596   else
597     {
598       fprintf (stderr, "g2 = ");
599       mpz_out_str (stderr, 16, p.x);
600       fprintf (stderr, "\n     ");
601       mpz_out_str (stderr, 16, p.y);
602       fprintf (stderr, "\n");
603     }
604   ecc_add (ecc, &q, &p, &ecc->g);
605   if (ecc->ref)
606     ASSERT_EQUAL (&q, &ecc->ref[1]);
607   else
608     {
609       fprintf (stderr, "g3 = ");
610       mpz_out_str (stderr, 16, q.x);
611       fprintf (stderr, "\n     ");
612       mpz_out_str (stderr, 16, q.y);
613       fprintf (stderr, "\n");
614     }
615
616   ecc_add (ecc, &q, &q, &ecc->g);
617   if (ecc->ref)
618     ASSERT_EQUAL (&q, &ecc->ref[2]);
619   else
620     {
621       fprintf (stderr, "g4 = ");
622       mpz_out_str (stderr, 16, q.x);
623       fprintf (stderr, "\n     ");
624       mpz_out_str (stderr, 16, q.y);
625       fprintf (stderr, "\n");
626     }
627
628   ecc_dup (ecc, &q, &p);
629   if (ecc->ref)
630     ASSERT_EQUAL (&q, &ecc->ref[2]);
631   else
632     {
633       fprintf (stderr, "g4 = ");
634       mpz_out_str (stderr, 16, q.x);
635       fprintf (stderr, "\n     ");
636       mpz_out_str (stderr, 16, q.y);
637       fprintf (stderr, "\n");
638     }
639
640   ecc_mul_binary (ecc, &p, ecc->q, &ecc->g);
641   ASSERT_ZERO (&p);
642
643   ecc_mul_pippenger (ecc, &q, ecc->q);
644   ASSERT_ZERO (&q);
645
646   ecc_clear (&p);
647   ecc_clear (&q);
648   mpz_clear (n);
649 }
650
651 static void
652 output_digits (const mpz_t x,
653                unsigned size, unsigned bits_per_limb)
654 {  
655   mpz_t t;
656   mpz_t mask;
657   mpz_t limb;
658   unsigned i;
659   const char *suffix;
660
661   mpz_init (t);
662   mpz_init (mask);
663   mpz_init (limb);
664
665   mpz_setbit (mask, bits_per_limb);
666   mpz_sub_ui (mask, mask, 1);
667
668   suffix = bits_per_limb > 32 ? "ULL" : "UL";
669
670   mpz_init_set (t, x);
671
672   for (i = 0; i < size; i++)
673     {
674       if ( (i % 8) == 0)
675         printf("\n ");
676       
677       mpz_and (limb, mask, t);
678       printf (" 0x");
679       mpz_out_str (stdout, 16, limb);
680       printf ("%s,", suffix);
681       mpz_tdiv_q_2exp (t, t, bits_per_limb);
682     }
683
684   mpz_clear (t);
685   mpz_clear (mask);
686   mpz_clear (limb);
687 }
688
689 static void
690 output_bignum (const char *name, const mpz_t x,
691                unsigned size, unsigned bits_per_limb)
692 {  
693   printf ("static const mp_limb_t %s[%d] = {", name, size);
694   output_digits (x, size, bits_per_limb);
695   printf("\n};\n");
696 }
697
698 static void
699 output_point (const char *name, const struct ecc_point *p,
700               unsigned size, unsigned bits_per_limb)
701 {
702   if (name)
703     printf("static const mp_limb_t %s[%u] = {", name, 2*size);
704
705   output_digits (p->x, size, bits_per_limb);
706   output_digits (p->y, size, bits_per_limb);
707
708   if (name)
709     printf("\n};\n");
710 }
711
712 static void
713 output_point_redc (const char *name, const struct ecc_curve *ecc,
714                    const struct ecc_point *p,
715                    unsigned size, unsigned bits_per_limb)
716 {
717   mpz_t t;
718   mpz_init (t);
719
720   if (name)
721     printf("static const mp_limb_t %s[%u] = {", name, 2*size);
722     
723   mpz_mul_2exp (t, p->x, size * bits_per_limb);
724   mpz_mod (t, t, ecc->p);
725       
726   output_digits (t, size, bits_per_limb);
727
728   mpz_mul_2exp (t, p->y, size * bits_per_limb);
729   mpz_mod (t, t, ecc->p);
730       
731   output_digits (t, size, bits_per_limb);
732
733   if (name)
734     printf("\n};\n");
735
736   mpz_clear (t);
737 }
738
739 static unsigned
740 output_modulo (const char *name, const mpz_t x,
741                unsigned size, unsigned bits_per_limb)
742 {
743   mpz_t mod;
744   unsigned bits;
745
746   mpz_init (mod);
747
748   mpz_setbit (mod, bits_per_limb * size);
749   mpz_mod (mod, mod, x);
750
751   bits = mpz_sizeinbase (mod, 2);
752   assert (bits <= size * bits_per_limb - 32);
753
754   output_bignum (name, mod, size, bits_per_limb);
755   
756   mpz_clear (mod);
757   return bits;
758 }
759
760 static void
761 output_curve (const struct ecc_curve *ecc, unsigned bits_per_limb)
762 {
763   unsigned limb_size = (ecc->bit_size + bits_per_limb - 1)/bits_per_limb;
764   unsigned i;
765   unsigned bits;
766   int redc_limbs;
767   mpz_t t;
768
769   mpz_init (t);
770
771   printf ("/* For NULL. */\n#include <stddef.h>\n");
772
773   printf ("#define ECC_LIMB_SIZE %u\n", limb_size);
774   printf ("#define ECC_PIPPENGER_K %u\n", ecc->pippenger_k);
775   printf ("#define ECC_PIPPENGER_C %u\n", ecc->pippenger_c);
776
777   output_bignum ("ecc_p", ecc->p, limb_size, bits_per_limb);
778   output_bignum ("ecc_b", ecc->b, limb_size, bits_per_limb);
779   output_bignum ("ecc_q", ecc->q, limb_size, bits_per_limb);
780   output_point ("ecc_g", &ecc->g, limb_size, bits_per_limb);
781   output_point_redc ("ecc_redc_g", ecc, &ecc->g, limb_size, bits_per_limb);
782   
783   bits = output_modulo ("ecc_Bmodp", ecc->p, limb_size, bits_per_limb);
784   printf ("#define ECC_BMODP_SIZE %u\n",
785           (bits + bits_per_limb - 1) / bits_per_limb);
786   bits = output_modulo ("ecc_Bmodq", ecc->q, limb_size, bits_per_limb);
787   printf ("#define ECC_BMODQ_SIZE %u\n",
788           (bits + bits_per_limb - 1) / bits_per_limb);
789
790   if (ecc->bit_size < limb_size * bits_per_limb)
791     {
792       int shift;
793
794       mpz_set_ui (t, 0);
795       mpz_setbit (t, ecc->bit_size);
796       mpz_sub (t, t, ecc->p);      
797       output_bignum ("ecc_Bmodp_shifted", t, limb_size, bits_per_limb);
798
799       shift = limb_size * bits_per_limb - ecc->bit_size;
800       if (shift > 0)
801         {
802           /* Check condition for reducing hi limbs. If s is the
803              normalization shift and n is the bit size (so that s + n
804              = limb_size * bite_per_limb), then we need
805
806                (2^n - 1) + (2^s - 1) (2^n - p) < 2p
807
808              or equivalently,
809
810                2^s (2^n - p) <= p
811
812              To a allow a carry limb to be added in at the same time,
813              substitute s+1 for s.
814           */
815           /* FIXME: For ecdsa verify, we actually need the stricter
816              inequality < 2 q. */
817           mpz_mul_2exp (t, t, shift + 1);
818           if (mpz_cmp (t, ecc->p) > 0)
819             {
820               fprintf (stderr, "Reduction condition failed for %u-bit curve.\n",
821                        ecc->bit_size);
822               exit (EXIT_FAILURE);
823             }
824         }
825       mpz_set_ui (t, 0);
826       mpz_setbit (t, ecc->bit_size);
827       mpz_sub (t, t, ecc->q);      
828       output_bignum ("ecc_Bmodq_shifted", t, limb_size, bits_per_limb);      
829     }
830   else
831     {
832       printf ("#define ecc_Bmodp_shifted ecc_Bmodp\n");
833       printf ("#define ecc_Bmodq_shifted ecc_Bmodq\n");
834     }
835
836   mpz_add_ui (t, ecc->p, 1);
837   mpz_fdiv_q_2exp (t, t, 1);
838   output_bignum ("ecc_pp1h", t, limb_size, bits_per_limb);      
839
840   mpz_add_ui (t, ecc->q, 1);
841   mpz_fdiv_q_2exp (t, t, 1);
842   output_bignum ("ecc_qp1h", t, limb_size, bits_per_limb);  
843   
844   /* Trailing zeros in p+1 correspond to trailing ones in p. */
845   redc_limbs = mpz_scan0 (ecc->p, 0) / bits_per_limb;
846   if (redc_limbs > 0)
847     {
848       mpz_add_ui (t, ecc->p, 1);
849       mpz_fdiv_q_2exp (t, t, redc_limbs * bits_per_limb);
850       output_bignum ("ecc_redc_ppm1", t, limb_size - redc_limbs, bits_per_limb);
851     }
852   else
853     {    
854       /* Trailing zeros in p-1 correspond to zeros just above the low
855          bit of p */
856       redc_limbs = mpz_scan1 (ecc->p, 1) / bits_per_limb;
857       if (redc_limbs > 0)
858         {
859           printf ("#define ecc_redc_ppm1 (ecc_p + %d)\n",
860                   redc_limbs);
861           redc_limbs = -redc_limbs;
862         }
863       else
864         printf ("#define ecc_redc_ppm1 NULL\n");
865     }
866   printf ("#define ECC_REDC_SIZE %d\n", redc_limbs);
867
868   printf ("#if USE_REDC\n");
869   printf ("#define ecc_unit ecc_Bmodp\n");
870
871   printf ("static const mp_limb_t ecc_table[%lu] = {",
872          (unsigned long) (2*ecc->table_size * limb_size));
873   for (i = 0; i < ecc->table_size; i++)
874     output_point_redc (NULL, ecc, &ecc->table[i], limb_size, bits_per_limb);
875
876   printf("\n};\n");
877
878   printf ("#else\n");
879
880   mpz_init_set_ui (t, 1);
881   output_bignum ("ecc_unit", t, limb_size, bits_per_limb);
882   
883   printf ("static const mp_limb_t ecc_table[%lu] = {",
884          (unsigned long) (2*ecc->table_size * limb_size));
885   for (i = 0; i < ecc->table_size; i++)
886     output_point (NULL, &ecc->table[i], limb_size, bits_per_limb);
887
888   printf("\n};\n");
889   printf ("#endif\n");
890   
891   mpz_clear (t);
892 }
893
894 int
895 main (int argc, char **argv)
896 {
897   struct ecc_curve ecc;
898
899   if (argc < 4)
900     {
901       fprintf (stderr, "Usage: %s CURVE-BITS K C [BITS-PER-LIMB]\n", argv[0]);
902       return EXIT_FAILURE;
903     }
904
905   ecc_curve_init (&ecc, atoi(argv[1]));
906
907   ecc_pippenger_precompute (&ecc, atoi(argv[2]), atoi(argv[3]));
908
909   fprintf (stderr, "Table size: %lu entries\n",
910            (unsigned long) ecc.table_size);
911
912   ecc_curve_check (&ecc);
913
914   if (argc > 4)
915     output_curve (&ecc, atoi(argv[4]));
916
917   return EXIT_SUCCESS;
918 }