crypto/openssl/crypto/ec/ecp_nistp521.c

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