git.samba.org - gd/nettle/blob - eccdata.c

   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
 103 ecc_init (struct ecc_point *p)
 104 {
 105   mpz_init (p->x);
 106   mpz_init (p->y);
 107 }
 108
 109 static void
 110 ecc_clear (struct ecc_point *p)
 111 {
 112   mpz_clear (p->x);
 113   mpz_clear (p->y);
 114 }
 115
 116 static int
 117 ecc_zero_p (const struct ecc_point *p)
 118 {
 119   return p->is_zero;
 120 }
 121
 122 static int
 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
 130 ecc_set_zero (struct ecc_point *r)
 131 {
 132   r->is_zero = 1;
 133 }
 134
 135 static void
 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
 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
 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
 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 *
 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
 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
 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
 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 unsigned
 615 ecc_table_size(unsigned bits, unsigned k, unsigned c)
 616 {
 617   unsigned p = (bits + k-1) / k;
 618   unsigned M = (p + c-1)/c;
 619   return M;
 620 }
 621
 622 static void
 623 ecc_pippenger_precompute (struct ecc_curve *ecc, unsigned k, unsigned c)
 624 {
 625   unsigned M = ecc_table_size (ecc->bit_size, k, c);
 626   unsigned i, j;
 627
 628   if (M == ecc_table_size (ecc->bit_size, k-1, c))
 629     fprintf(stderr,
 630             "warn: Parameters k = %u, c = %d are suboptimal, could use smaller k\n",
 631             k, c);
 632
 633   ecc->pippenger_k = k;
 634   ecc->pippenger_c = c;
 635   ecc->table_size = M << c;
 636   ecc->table = ecc_alloc (ecc->table_size);
 637
 638   /* Compute the first 2^c entries */
 639   ecc_set_zero (&ecc->table[0]);
 640   ecc_set (&ecc->table[1], &ecc->g);
 641
 642   for (j = 2; j < (1U<<c); j <<= 1)
 643     {
 644       /* T[j] = 2^k T[j/2] */
 645       ecc_dup (ecc, &ecc->table[j], &ecc->table[j/2]);
 646       for (i = 1; i < k; i++)
 647         ecc_dup (ecc, &ecc->table[j], &ecc->table[j]);
 648
 649       for (i = 1; i < j; i++)
 650         ecc_add (ecc, &ecc->table[j + i], &ecc->table[j], &ecc->table[i]);
 651     }
 652   for (j = 1<<c; j < ecc->table_size; j++)
 653     {
 654       /* T[j] = 2^{kc} T[j-2^c] */
 655       ecc_dup (ecc, &ecc->table[j], &ecc->table[j - (1<<c)]);
 656       for (i = 1; i < k*c; i++)
 657         ecc_dup (ecc, &ecc->table[j], &ecc->table[j]);
 658     }
 659 }
 660
 661 static void
 662 ecc_mul_pippenger (const struct ecc_curve *ecc,
 663                    struct ecc_point *r, const mpz_t n_input)
 664 {
 665   mpz_t n;
 666   unsigned k, c;
 667   unsigned i, j;
 668   unsigned bit_rows;
 669
 670   mpz_init (n);
 671
 672   mpz_mod (n, n_input, ecc->q);
 673   ecc_set_zero (r);
 674
 675   k = ecc->pippenger_k;
 676   c = ecc->pippenger_c;
 677
 678   bit_rows = (ecc->bit_size + k - 1) / k;
 679
 680   for (i = k; i-- > 0; )
 681     {
 682       ecc_dup (ecc, r, r);
 683       for (j = 0; j * c < bit_rows; j++)
 684         {
 685           unsigned bits;
 686           mp_size_t bit_index;
 687
 688           /* Extract c bits of the exponent, stride k, starting at i + kcj, ending at
 689             i + k (cj + c - 1)*/
 690           for (bits = 0, bit_index = i + k*(c*j+c); bit_index > i + k*c*j; )
 691             {
 692               bit_index -= k;
 693               bits = (bits << 1) | mpz_tstbit (n, bit_index);
 694             }
 695
 696           ecc_add (ecc, r, r, &ecc->table[(j << c) | bits]);
 697         }
 698     }
 699   mpz_clear (n);
 700 }
 701
 702 static void
 703 ecc_point_out (FILE *f, const struct ecc_point *p)
 704 {
 705   if (p->is_zero)
 706     fprintf (f, "zero");
 707   else
 708     {
 709         fprintf (f, "(");
 710         mpz_out_str (f, 16, p->x);
 711         fprintf (f, ",\n     ");
 712         mpz_out_str (f, 16, (p)->y);
 713         fprintf (f, ")");
 714     }
 715 }
 716 #define ASSERT_EQUAL(p, q) do {                                         \
 717     if (!ecc_equal_p (p, q))                                            \
 718       {                                                                 \
 719         fprintf (stderr, "%s:%d: ASSERT_EQUAL (%s, %s) failed.\n",      \
 720                  __FILE__, __LINE__, #p, #q);                           \
 721         fprintf (stderr, "p = ");                                       \
 722         ecc_point_out (stderr, (p));                                    \
 723         fprintf (stderr, "\nq = ");                                     \
 724         ecc_point_out (stderr, (q));                                    \
 725         fprintf (stderr, "\n");                                         \
 726         abort();                                                        \
 727       }                                                                 \
 728   } while (0)
 729
 730 #define ASSERT_ZERO(p) do {                                             \
 731     if (!ecc_zero_p (p))                                                \
 732       {                                                                 \
 733         fprintf (stderr, "%s:%d: ASSERT_ZERO (%s) failed.\n",           \
 734                  __FILE__, __LINE__, #p);                               \
 735         fprintf (stderr, "p = ");                                       \
 736         ecc_point_out (stderr, (p));                                    \
 737         fprintf (stderr, "\n");                                         \
 738         abort();                                                        \
 739       }                                                                 \
 740   } while (0)
 741
 742 static void
 743 ecc_curve_check (const struct ecc_curve *ecc)
 744 {
 745   struct ecc_point p, q;
 746   mpz_t n;
 747
 748   ecc_init (&p);
 749   ecc_init (&q);
 750   mpz_init (n);
 751
 752   ecc_dup (ecc, &p, &ecc->g);
 753   if (ecc->ref)
 754     ASSERT_EQUAL (&p, &ecc->ref[0]);
 755   else
 756     {
 757       fprintf (stderr, "g2 = ");
 758       mpz_out_str (stderr, 16, p.x);
 759       fprintf (stderr, "\n     ");
 760       mpz_out_str (stderr, 16, p.y);
 761       fprintf (stderr, "\n");
 762     }
 763   ecc_add (ecc, &q, &p, &ecc->g);
 764   if (ecc->ref)
 765     ASSERT_EQUAL (&q, &ecc->ref[1]);
 766   else
 767     {
 768       fprintf (stderr, "g3 = ");
 769       mpz_out_str (stderr, 16, q.x);
 770       fprintf (stderr, "\n     ");
 771       mpz_out_str (stderr, 16, q.y);
 772       fprintf (stderr, "\n");
 773     }
 774
 775   ecc_add (ecc, &q, &q, &ecc->g);
 776   if (ecc->ref)
 777     ASSERT_EQUAL (&q, &ecc->ref[2]);
 778   else
 779     {
 780       fprintf (stderr, "g4 = ");
 781       mpz_out_str (stderr, 16, q.x);
 782       fprintf (stderr, "\n     ");
 783       mpz_out_str (stderr, 16, q.y);
 784       fprintf (stderr, "\n");
 785     }
 786
 787   ecc_dup (ecc, &q, &p);
 788   if (ecc->ref)
 789     ASSERT_EQUAL (&q, &ecc->ref[2]);
 790   else
 791     {
 792       fprintf (stderr, "g4 = ");
 793       mpz_out_str (stderr, 16, q.x);
 794       fprintf (stderr, "\n     ");
 795       mpz_out_str (stderr, 16, q.y);
 796       fprintf (stderr, "\n");
 797     }
 798
 799   ecc_mul_binary (ecc, &p, ecc->q, &ecc->g);
 800   ASSERT_ZERO (&p);
 801
 802   ecc_mul_pippenger (ecc, &q, ecc->q);
 803   ASSERT_ZERO (&q);
 804
 805   ecc_clear (&p);
 806   ecc_clear (&q);
 807   mpz_clear (n);
 808 }
 809
 810 static void
 811 output_digits (const mpz_t x,
 812                unsigned size, unsigned bits_per_limb)
 813 {
 814   mpz_t t;
 815   mpz_t mask;
 816   mpz_t limb;
 817   unsigned i;
 818   const char *suffix;
 819
 820   mpz_init (t);
 821   mpz_init (mask);
 822   mpz_init (limb);
 823
 824   mpz_setbit (mask, bits_per_limb);
 825   mpz_sub_ui (mask, mask, 1);
 826
 827   suffix = bits_per_limb > 32 ? "ULL" : "UL";
 828
 829   mpz_init_set (t, x);
 830
 831   for (i = 0; i < size; i++)
 832     {
 833       if ( (i % 8) == 0)
 834         printf("\n ");
 835
 836       mpz_and (limb, mask, t);
 837       printf (" 0x");
 838       mpz_out_str (stdout, 16, limb);
 839       printf ("%s,", suffix);
 840       mpz_tdiv_q_2exp (t, t, bits_per_limb);
 841     }
 842
 843   mpz_clear (t);
 844   mpz_clear (mask);
 845   mpz_clear (limb);
 846 }
 847
 848 static void
 849 output_bignum (const char *name, const mpz_t x,
 850                unsigned size, unsigned bits_per_limb)
 851 {
 852   printf ("static const mp_limb_t %s[%d] = {", name, size);
 853   output_digits (x, size, bits_per_limb);
 854   printf("\n};\n");
 855 }
 856
 857 static void
 858 output_point (const char *name, const struct ecc_curve *ecc,
 859               const struct ecc_point *p, int use_redc,
 860               unsigned size, unsigned bits_per_limb)
 861 {
 862   mpz_t x, y, t;
 863
 864   mpz_init (x);
 865   mpz_init (y);
 866   mpz_init (t);
 867
 868   if (name)
 869     printf("static const mp_limb_t %s[%u] = {", name, 2*size);
 870
 871   if (ecc->use_edwards)
 872     {
 873       if (ecc_zero_p (p))
 874         {
 875           mpz_set_si (x, 0);
 876           mpz_set_si (y, 1);
 877         }
 878       else if (!mpz_sgn (p->y))
 879         {
 880           assert (!mpz_sgn (p->x));
 881           mpz_set_si (x, 0);
 882           mpz_set_si (y, -1);
 883         }
 884       else
 885         {
 886           mpz_invert (x, p->y, ecc->p);
 887           mpz_mul (x, x, p->x);
 888           mpz_mul (x, x, ecc->t);
 889           mpz_mod (x, x, ecc->p);
 890
 891           mpz_sub_ui (y, p->x, 1);
 892           mpz_add_ui (t, p->x, 1);
 893           mpz_invert (t, t, ecc->p);
 894           mpz_mul (y, y, t);
 895           mpz_mod (y, y, ecc->p);
 896         }
 897     }
 898   else
 899     {
 900       mpz_set (x, p->x);
 901       mpz_set (y, p->y);
 902     }
 903   if (use_redc)
 904     {
 905       mpz_mul_2exp (x, x, size * bits_per_limb);
 906       mpz_mod (x, x, ecc->p);
 907       mpz_mul_2exp (y, y, size * bits_per_limb);
 908       mpz_mod (y, y, ecc->p);
 909     }
 910
 911   output_digits (x, size, bits_per_limb);
 912   output_digits (y, size, bits_per_limb);
 913
 914   if (name)
 915     printf("\n};\n");
 916
 917   mpz_clear (x);
 918   mpz_clear (y);
 919   mpz_clear (t);
 920 }
 921
 922 static unsigned
 923 output_modulo (const char *name, const mpz_t x,
 924                unsigned size, unsigned bits_per_limb)
 925 {
 926   mpz_t mod;
 927   unsigned bits;
 928
 929   mpz_init (mod);
 930
 931   mpz_setbit (mod, bits_per_limb * size);
 932   mpz_mod (mod, mod, x);
 933
 934   bits = mpz_sizeinbase (mod, 2);
 935   output_bignum (name, mod, size, bits_per_limb);
 936
 937   mpz_clear (mod);
 938   return bits;
 939 }
 940
 941 static void
 942 output_curve (const struct ecc_curve *ecc, unsigned bits_per_limb)
 943 {
 944   unsigned limb_size = (ecc->bit_size + bits_per_limb - 1)/bits_per_limb;
 945   unsigned i;
 946   unsigned bits, e;
 947   int redc_limbs;
 948   mpz_t t;
 949
 950   mpz_init (t);
 951
 952   printf ("/* For NULL. */\n#include <stddef.h>\n");
 953
 954   printf ("#define ECC_LIMB_SIZE %u\n", limb_size);
 955   printf ("#define ECC_PIPPENGER_K %u\n", ecc->pippenger_k);
 956   printf ("#define ECC_PIPPENGER_C %u\n", ecc->pippenger_c);
 957
 958   output_bignum ("ecc_p", ecc->p, limb_size, bits_per_limb);
 959   output_bignum ("ecc_b", ecc->b, limb_size, bits_per_limb);
 960   if (ecc->use_edwards)
 961     output_bignum ("ecc_d", ecc->d, limb_size, bits_per_limb);
 962   output_bignum ("ecc_q", ecc->q, limb_size, bits_per_limb);
 963   output_point ("ecc_g", ecc, &ecc->g, 0, limb_size, bits_per_limb);
 964
 965   bits = output_modulo ("ecc_Bmodp", ecc->p, limb_size, bits_per_limb);
 966   printf ("#define ECC_BMODP_SIZE %u\n",
 967           (bits + bits_per_limb - 1) / bits_per_limb);
 968   bits = output_modulo ("ecc_Bmodq", ecc->q, limb_size, bits_per_limb);
 969   printf ("#define ECC_BMODQ_SIZE %u\n",
 970           (bits + bits_per_limb - 1) / bits_per_limb);
 971   bits = mpz_sizeinbase (ecc->q, 2);
 972   if (bits < ecc->bit_size)
 973     {
 974       /* for curve25519, with q = 2^k + q', with a much smaller q' */
 975       unsigned mbits;
 976       unsigned shift;
 977
 978       /* Shift to align the one bit at B */
 979       shift = bits_per_limb * limb_size + 1 - bits;
 980
 981       mpz_set (t, ecc->q);
 982       mpz_clrbit (t, bits-1);
 983       mbits = mpz_sizeinbase (t, 2);
 984
 985       /* The shifted value must be a limb smaller than q. */
 986       if (mbits + shift + bits_per_limb <= bits)
 987         {
 988           /* q of the form 2^k + q', with q' a limb smaller */
 989           mpz_mul_2exp (t, t, shift);
 990           output_bignum ("ecc_mBmodq_shifted", t, limb_size, bits_per_limb);
 991         }
 992     }
 993
 994   if (ecc->bit_size < limb_size * bits_per_limb)
 995     {
 996       int shift;
 997
 998       mpz_set_ui (t, 0);
 999       mpz_setbit (t, ecc->bit_size);
1000       mpz_sub (t, t, ecc->p);
1001       output_bignum ("ecc_Bmodp_shifted", t, limb_size, bits_per_limb);
1002
1003       shift = limb_size * bits_per_limb - ecc->bit_size;
1004       if (shift > 0)
1005         {
1006           /* Check condition for reducing hi limbs. If s is the
1007              normalization shift and n is the bit size (so that s + n
1008              = limb_size * bite_per_limb), then we need
1009
1010                (2^n - 1) + (2^s - 1) (2^n - p) < 2p
1011
1012              or equivalently,
1013
1014                2^s (2^n - p) <= p
1015
1016              To a allow a carry limb to be added in at the same time,
1017              substitute s+1 for s.
1018           */
1019           /* FIXME: For ecdsa verify, we actually need the stricter
1020              inequality < 2 q. */
1021           mpz_mul_2exp (t, t, shift + 1);
1022           if (mpz_cmp (t, ecc->p) > 0)
1023             {
1024               fprintf (stderr, "Reduction condition failed for %u-bit curve.\n",
1025                        ecc->bit_size);
1026               exit (EXIT_FAILURE);
1027             }
1028         }
1029     }
1030   else
1031     printf ("#define ecc_Bmodp_shifted ecc_Bmodp\n");
1032
1033   if (bits < limb_size * bits_per_limb)
1034     {
1035       mpz_set_ui (t, 0);
1036       mpz_setbit (t, bits);
1037       mpz_sub (t, t, ecc->q);
1038       output_bignum ("ecc_Bmodq_shifted", t, limb_size, bits_per_limb);
1039     }
1040   else
1041     printf ("#define ecc_Bmodq_shifted ecc_Bmodq\n");
1042
1043   mpz_add_ui (t, ecc->p, 1);
1044   mpz_fdiv_q_2exp (t, t, 1);
1045   output_bignum ("ecc_pp1h", t, limb_size, bits_per_limb);
1046
1047   mpz_add_ui (t, ecc->q, 1);
1048   mpz_fdiv_q_2exp (t, t, 1);
1049   output_bignum ("ecc_qp1h", t, limb_size, bits_per_limb);
1050
1051   if (ecc->use_edwards)
1052     output_bignum ("ecc_edwards", ecc->t, limb_size, bits_per_limb);
1053
1054   /* Trailing zeros in p+1 correspond to trailing ones in p. */
1055   redc_limbs = mpz_scan0 (ecc->p, 0) / bits_per_limb;
1056   if (redc_limbs > 0)
1057     {
1058       mpz_add_ui (t, ecc->p, 1);
1059       mpz_fdiv_q_2exp (t, t, redc_limbs * bits_per_limb);
1060       output_bignum ("ecc_redc_ppm1", t, limb_size - redc_limbs, bits_per_limb);
1061     }
1062   else
1063     {
1064       /* Trailing zeros in p-1 correspond to zeros just above the low
1065          bit of p */
1066       redc_limbs = mpz_scan1 (ecc->p, 1) / bits_per_limb;
1067       if (redc_limbs > 0)
1068         {
1069           printf ("#define ecc_redc_ppm1 (ecc_p + %d)\n",
1070                   redc_limbs);
1071           redc_limbs = -redc_limbs;
1072         }
1073       else
1074         printf ("#define ecc_redc_ppm1 NULL\n");
1075     }
1076   printf ("#define ECC_REDC_SIZE %d\n", redc_limbs);
1077
1078   /* For mod p square root computation. */
1079   if (mpz_fdiv_ui (ecc->p, 4) == 3)
1080     {
1081       /* x = a^{(p+1)/4} gives square root of a (if it exists,
1082          otherwise the square root of -a). */
1083       e = 1;
1084       mpz_add_ui (t, ecc->p, 1);
1085       mpz_fdiv_q_2exp (t, t, 2);
1086     }
1087   else
1088     {
1089       /* p-1 = 2^e s, s odd, t = (s-1)/2*/
1090       unsigned g, i;
1091       mpz_t s;
1092       mpz_t z;
1093
1094       mpz_init (s);
1095       mpz_init (z);
1096
1097       mpz_sub_ui (s, ecc->p, 1);
1098       e = mpz_scan1 (s, 0);
1099       assert (e > 1);
1100
1101       mpz_fdiv_q_2exp (s, s, e);
1102
1103       /* Find a non-square g, g^{(p-1)/2} = -1,
1104          and z = g^{(p-1)/4 */
1105       for (g = 2; ; g++)
1106         {
1107           mpz_set_ui (z, g);
1108           mpz_powm (z, z, s, ecc->p);
1109           mpz_mul (t, z, z);
1110           mpz_mod (t, t, ecc->p);
1111
1112           for (i = 2; i < e; i++)
1113             {
1114               mpz_mul (t, t, t);
1115               mpz_mod (t, t, ecc->p);
1116             }
1117           if (mpz_cmp_ui (t, 1) != 0)
1118             break;
1119         }
1120       mpz_add_ui (t, t, 1);
1121       assert (mpz_cmp (t, ecc->p) == 0);
1122       output_bignum ("ecc_sqrt_z", z, limb_size, bits_per_limb);
1123
1124       mpz_fdiv_q_2exp (t, s, 1);
1125
1126       mpz_clear (s);
1127       mpz_clear (z);
1128     }
1129   printf ("#define ECC_SQRT_E %u\n", e);
1130   printf ("#define ECC_SQRT_T_BITS %u\n",
1131           (unsigned) mpz_sizeinbase (t, 2));
1132   output_bignum ("ecc_sqrt_t", t, limb_size, bits_per_limb);
1133
1134   printf ("#if USE_REDC\n");
1135   printf ("#define ecc_unit ecc_Bmodp\n");
1136
1137   printf ("static const mp_limb_t ecc_table[%lu] = {",
1138          (unsigned long) (2*ecc->table_size * limb_size));
1139   for (i = 0; i < ecc->table_size; i++)
1140     output_point (NULL, ecc, &ecc->table[i], 1, limb_size, bits_per_limb);
1141
1142   printf("\n};\n");
1143
1144   printf ("#else\n");
1145
1146   mpz_set_ui (t, 1);
1147   output_bignum ("ecc_unit", t, limb_size, bits_per_limb);
1148
1149   printf ("static const mp_limb_t ecc_table[%lu] = {",
1150          (unsigned long) (2*ecc->table_size * limb_size));
1151   for (i = 0; i < ecc->table_size; i++)
1152     output_point (NULL, ecc, &ecc->table[i], 0, limb_size, bits_per_limb);
1153
1154   printf("\n};\n");
1155   printf ("#endif\n");
1156
1157   mpz_clear (t);
1158 }
1159
1160 int
1161 main (int argc, char **argv)
1162 {
1163   struct ecc_curve ecc;
1164
1165   if (argc < 4)
1166     {
1167       fprintf (stderr, "Usage: %s CURVE-BITS K C [BITS-PER-LIMB]\n", argv[0]);
1168       return EXIT_FAILURE;
1169     }
1170
1171   ecc_curve_init (&ecc, atoi(argv[1]));
1172
1173   ecc_pippenger_precompute (&ecc, atoi(argv[2]), atoi(argv[3]));
1174
1175   fprintf (stderr, "Table size: %lu entries\n",
1176            (unsigned long) ecc.table_size);
1177
1178   ecc_curve_check (&ecc);
1179
1180   if (argc > 4)
1181     output_curve (&ecc, atoi(argv[4]));
1182
1183   return EXIT_SUCCESS;
1184 }