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 void
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
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 }