crypto/ec/ecp_nistp224.c

   1 /* crypto/ec/ecp_nistp224.c */
   2 /*
   3  * Written by Emilia Kasper (Google) for the OpenSSL project.
   4  */
   5 /* ====================================================================
   6  * Copyright (c) 2000-2010 The OpenSSL Project.  All rights reserved.
   7  *
   8  * Redistribution and use in source and binary forms, with or without
   9  * modification, are permitted provided that the following conditions
  10  * are met:
  11  *
  12  * 1. Redistributions of source code must retain the above copyright
  13  *    notice, this list of conditions and the following disclaimer.
  14  *
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in
  17  *    the documentation and/or other materials provided with the
  18  *    distribution.
  19  *
  20  * 3. All advertising materials mentioning features or use of this
  21  *    software must display the following acknowledgment:
  22  *    "This product includes software developed by the OpenSSL Project
  23  *    for use in the OpenSSL Toolkit. (http://www.OpenSSL.org/)"
  24  *
  25  * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
  26  *    endorse or promote products derived from this software without
  27  *    prior written permission. For written permission, please contact
  28  *    licensing@OpenSSL.org.
  29  *
  30  * 5. Products derived from this software may not be called "OpenSSL"
  31  *    nor may "OpenSSL" appear in their names without prior written
  32  *    permission of the OpenSSL Project.
  33  *
  34  * 6. Redistributions of any form whatsoever must retain the following
  35  *    acknowledgment:
  36  *    "This product includes software developed by the OpenSSL Project
  37  *    for use in the OpenSSL Toolkit (http://www.OpenSSL.org/)"
  38  *
  39  * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
  40  * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  41  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  42  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
  43  * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  44  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  45  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  46  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  47  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  48  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  49  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
  50  * OF THE POSSIBILITY OF SUCH DAMAGE.
  51  * ====================================================================
  52  *
  53  * This product includes cryptographic software written by Eric Young
  54  * (eay@cryptsoft.com).  This product includes software written by Tim
  55  * Hudson (tjh@cryptsoft.com).
  56  *
  57  */
  58
  59 /*
  60  * A 64-bit implementation of the NIST P-224 elliptic curve point multiplication
  61  *
  62  * Inspired by Daniel J. Bernstein's public domain nistp224 implementation
  63  * and Adam Langley's public domain 64-bit C implementation of curve25519
  64  */
  65 #include <openssl/opensslconf.h>
  66 #ifndef OPENSSL_NO_EC_NISTP224_64_GCC_128
  67 #include <stdint.h>
  68 #include <string.h>
  69 #include <openssl/err.h>
  70 #include "ec_lcl.h"
  71
  72 #if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 1))
  73   /* even with gcc, the typedef won't work for 32-bit platforms */
  74   typedef __uint128_t uint128_t; /* nonstandard; implemented by gcc on 64-bit platforms */
  75 #else
  76   #error "Need GCC 3.1 or later to define type uint128_t"
  77 #endif
  78
  79 typedef uint8_t u8;
  80
  81
  82 /******************************************************************************/
  83 /*                  INTERNAL REPRESENTATION OF FIELD ELEMENTS
  84  *
  85  * Field elements are represented as a_0 + 2^56*a_1 + 2^112*a_2 + 2^168*a_3
  86  * where each slice a_i is a 64-bit word, i.e., a field element is an fslice
  87  * array a with 4 elements, where a[i] = a_i.
  88  * Outputs from multiplications are represented as unreduced polynomials
  89  * b_0 + 2^56*b_1 + 2^112*b_2 + 2^168*b_3 + 2^224*b_4 + 2^280*b_5 + 2^336*b_6
  90  * where each b_i is a 128-bit word. We ensure that inputs to each field
  91  * multiplication satisfy a_i < 2^60, so outputs satisfy b_i < 4*2^60*2^60,
  92  * and fit into a 128-bit word without overflow. The coefficients are then
  93  * again partially reduced to a_i < 2^57. We only reduce to the unique minimal
  94  * representation at the end of the computation.
  95  *
  96  */
  97
  98 typedef uint64_t fslice;
  99
 100 /* Field element represented as a byte arrary.
 101  * 28*8 = 224 bits is also the group order size for the elliptic curve.  */
 102 typedef u8 felem_bytearray[28];
 103
 104 static const felem_bytearray nistp224_curve_params[5] = {
 105         {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,    /* p */
 106          0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0x00,0x00,0x00,0x00,
 107          0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x01},
 108         {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,    /* a */
 109          0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,0xFF,0xFF,0xFF,0xFF,
 110          0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFE},
 111         {0xB4,0x05,0x0A,0x85,0x0C,0x04,0xB3,0xAB,0xF5,0x41,    /* b */
 112          0x32,0x56,0x50,0x44,0xB0,0xB7,0xD7,0xBF,0xD8,0xBA,
 113          0x27,0x0B,0x39,0x43,0x23,0x55,0xFF,0xB4},
 114         {0xB7,0x0E,0x0C,0xBD,0x6B,0xB4,0xBF,0x7F,0x32,0x13,    /* x */
 115          0x90,0xB9,0x4A,0x03,0xC1,0xD3,0x56,0xC2,0x11,0x22,
 116          0x34,0x32,0x80,0xD6,0x11,0x5C,0x1D,0x21},
 117         {0xbd,0x37,0x63,0x88,0xb5,0xf7,0x23,0xfb,0x4c,0x22,    /* y */
 118          0xdf,0xe6,0xcd,0x43,0x75,0xa0,0x5a,0x07,0x47,0x64,
 119          0x44,0xd5,0x81,0x99,0x85,0x00,0x7e,0x34}
 120 };
 121
 122 /* Precomputed multiples of the standard generator
 123  * b_0*G + b_1*2^56*G + b_2*2^112*G + b_3*2^168*G for
 124  * (b_3, b_2, b_1, b_0) in [0,15], i.e., gmul[0] = point_at_infinity,
 125  * gmul[1] = G, gmul[2] = 2^56*G, gmul[3] = 2^56*G + G, etc.
 126  * Points are given in Jacobian projective coordinates: words 0-3 represent the
 127  * X-coordinate (slice a_0 is word 0, etc.), words 4-7 represent the
 128  * Y-coordinate and words 8-11 represent the Z-coordinate. */
 129 static const fslice gmul[16][3][4] = {
 130         {{0x00000000000000, 0x00000000000000, 0x00000000000000, 0x00000000000000},
 131          {0x00000000000000, 0x00000000000000, 0x00000000000000, 0x00000000000000},
 132          {0x00000000000000, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 133         {{0x3280d6115c1d21, 0xc1d356c2112234, 0x7f321390b94a03, 0xb70e0cbd6bb4bf},
 134          {0xd5819985007e34, 0x75a05a07476444, 0xfb4c22dfe6cd43, 0xbd376388b5f723},
 135          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 136         {{0xfd9675666ebbe9, 0xbca7664d40ce5e, 0x2242df8d8a2a43, 0x1f49bbb0f99bc5},
 137          {0x29e0b892dc9c43, 0xece8608436e662, 0xdc858f185310d0, 0x9812dd4eb8d321},
 138          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 139         {{0x6d3e678d5d8eb8, 0x559eed1cb362f1, 0x16e9a3bbce8a3f, 0xeedcccd8c2a748},
 140          {0xf19f90ed50266d, 0xabf2b4bf65f9df, 0x313865468fafec, 0x5cb379ba910a17},
 141          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 142         {{0x0641966cab26e3, 0x91fb2991fab0a0, 0xefec27a4e13a0b, 0x0499aa8a5f8ebe},
 143          {0x7510407766af5d, 0x84d929610d5450, 0x81d77aae82f706, 0x6916f6d4338c5b},
 144          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 145         {{0xea95ac3b1f15c6, 0x086000905e82d4, 0xdd323ae4d1c8b1, 0x932b56be7685a3},
 146          {0x9ef93dea25dbbf, 0x41665960f390f0, 0xfdec76dbe2a8a7, 0x523e80f019062a},
 147          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 148         {{0x822fdd26732c73, 0xa01c83531b5d0f, 0x363f37347c1ba4, 0xc391b45c84725c},
 149          {0xbbd5e1b2d6ad24, 0xddfbcde19dfaec, 0xc393da7e222a7f, 0x1efb7890ede244},
 150          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 151         {{0x4c9e90ca217da1, 0xd11beca79159bb, 0xff8d33c2c98b7c, 0x2610b39409f849},
 152          {0x44d1352ac64da0, 0xcdbb7b2c46b4fb, 0x966c079b753c89, 0xfe67e4e820b112},
 153          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 154         {{0xe28cae2df5312d, 0xc71b61d16f5c6e, 0x79b7619a3e7c4c, 0x05c73240899b47},
 155          {0x9f7f6382c73e3a, 0x18615165c56bda, 0x641fab2116fd56, 0x72855882b08394},
 156          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 157         {{0x0469182f161c09, 0x74a98ca8d00fb5, 0xb89da93489a3e0, 0x41c98768fb0c1d},
 158          {0xe5ea05fb32da81, 0x3dce9ffbca6855, 0x1cfe2d3fbf59e6, 0x0e5e03408738a7},
 159          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 160         {{0xdab22b2333e87f, 0x4430137a5dd2f6, 0xe03ab9f738beb8, 0xcb0c5d0dc34f24},
 161          {0x764a7df0c8fda5, 0x185ba5c3fa2044, 0x9281d688bcbe50, 0xc40331df893881},
 162          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 163         {{0xb89530796f0f60, 0xade92bd26909a3, 0x1a0c83fb4884da, 0x1765bf22a5a984},
 164          {0x772a9ee75db09e, 0x23bc6c67cec16f, 0x4c1edba8b14e2f, 0xe2a215d9611369},
 165          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 166         {{0x571e509fb5efb3, 0xade88696410552, 0xc8ae85fada74fe, 0x6c7e4be83bbde3},
 167          {0xff9f51160f4652, 0xb47ce2495a6539, 0xa2946c53b582f4, 0x286d2db3ee9a60},
 168          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 169         {{0x40bbd5081a44af, 0x0995183b13926c, 0xbcefba6f47f6d0, 0x215619e9cc0057},
 170          {0x8bc94d3b0df45e, 0xf11c54a3694f6f, 0x8631b93cdfe8b5, 0xe7e3f4b0982db9},
 171          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 172         {{0xb17048ab3e1c7b, 0xac38f36ff8a1d8, 0x1c29819435d2c6, 0xc813132f4c07e9},
 173          {0x2891425503b11f, 0x08781030579fea, 0xf5426ba5cc9674, 0x1e28ebf18562bc},
 174          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}},
 175         {{0x9f31997cc864eb, 0x06cd91d28b5e4c, 0xff17036691a973, 0xf1aef351497c58},
 176          {0xdd1f2d600564ff, 0xdead073b1402db, 0x74a684435bd693, 0xeea7471f962558},
 177          {0x00000000000001, 0x00000000000000, 0x00000000000000, 0x00000000000000}}
 178 };
 179
 180 /* Precomputation for the group generator. */
 181 typedef struct {
 182         fslice g_pre_comp[16][3][4];
 183         int references;
 184 } NISTP224_PRE_COMP;
 185
 186 const EC_METHOD *EC_GFp_nistp224_method(void)
 187         {
 188         static const EC_METHOD ret = {
 189                 NID_X9_62_prime_field,
 190                 ec_GFp_nistp224_group_init,
 191                 ec_GFp_simple_group_finish,
 192                 ec_GFp_simple_group_clear_finish,
 193                 ec_GFp_nist_group_copy,
 194                 ec_GFp_nistp224_group_set_curve,
 195                 ec_GFp_simple_group_get_curve,
 196                 ec_GFp_simple_group_get_degree,
 197                 ec_GFp_simple_group_check_discriminant,
 198                 ec_GFp_simple_point_init,
 199                 ec_GFp_simple_point_finish,
 200                 ec_GFp_simple_point_clear_finish,
 201                 ec_GFp_simple_point_copy,
 202                 ec_GFp_simple_point_set_to_infinity,
 203                 ec_GFp_simple_set_Jprojective_coordinates_GFp,
 204                 ec_GFp_simple_get_Jprojective_coordinates_GFp,
 205                 ec_GFp_simple_point_set_affine_coordinates,
 206                 ec_GFp_nistp224_point_get_affine_coordinates,
 207                 ec_GFp_simple_set_compressed_coordinates,
 208                 ec_GFp_simple_point2oct,
 209                 ec_GFp_simple_oct2point,
 210                 ec_GFp_simple_add,
 211                 ec_GFp_simple_dbl,
 212                 ec_GFp_simple_invert,
 213                 ec_GFp_simple_is_at_infinity,
 214                 ec_GFp_simple_is_on_curve,
 215                 ec_GFp_simple_cmp,
 216                 ec_GFp_simple_make_affine,
 217                 ec_GFp_simple_points_make_affine,
 218                 ec_GFp_nistp224_points_mul,
 219                 ec_GFp_nistp224_precompute_mult,
 220                 ec_GFp_nistp224_have_precompute_mult,
 221                 ec_GFp_nist_field_mul,
 222                 ec_GFp_nist_field_sqr,
 223                 0 /* field_div */,
 224                 0 /* field_encode */,
 225                 0 /* field_decode */,
 226                 0 /* field_set_to_one */ };
 227
 228         return &ret;
 229         }
 230
 231 /* Helper functions to convert field elements to/from internal representation */
 232 static void bin28_to_felem(fslice out[4], const u8 in[28])
 233         {
 234         out[0] = *((const uint64_t *)(in)) & 0x00ffffffffffffff;
 235         out[1] = (*((const uint64_t *)(in+7))) & 0x00ffffffffffffff;
 236         out[2] = (*((const uint64_t *)(in+14))) & 0x00ffffffffffffff;
 237         out[3] = (*((const uint64_t *)(in+21))) & 0x00ffffffffffffff;
 238         }
 239
 240 static void felem_to_bin28(u8 out[28], const fslice in[4])
 241         {
 242         unsigned i;
 243         for (i = 0; i < 7; ++i)
 244                 {
 245                 out[i]    = in[0]>>(8*i);
 246                 out[i+7]  = in[1]>>(8*i);
 247                 out[i+14] = in[2]>>(8*i);
 248                 out[i+21] = in[3]>>(8*i);
 249                 }
 250         }
 251
 252 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */
 253 static void flip_endian(u8 *out, const u8 *in, unsigned len)
 254         {
 255         unsigned i;
 256         for (i = 0; i < len; ++i)
 257                 out[i] = in[len-1-i];
 258         }
 259
 260 /* From OpenSSL BIGNUM to internal representation */
 261 static int BN_to_felem(fslice out[4], const BIGNUM *bn)
 262         {
 263         felem_bytearray b_in;
 264         felem_bytearray b_out;
 265         unsigned num_bytes;
 266
 267         /* BN_bn2bin eats leading zeroes */
 268         memset(b_out, 0, sizeof b_out);
 269         num_bytes = BN_num_bytes(bn);
 270         if (num_bytes > sizeof b_out)
 271                 {
 272                 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
 273                 return 0;
 274                 }
 275         if (BN_is_negative(bn))
 276                 {
 277                 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
 278                 return 0;
 279                 }
 280         num_bytes = BN_bn2bin(bn, b_in);
 281         flip_endian(b_out, b_in, num_bytes);
 282         bin28_to_felem(out, b_out);
 283         return 1;
 284         }
 285
 286 /* From internal representation to OpenSSL BIGNUM */
 287 static BIGNUM *felem_to_BN(BIGNUM *out, const fslice in[4])
 288         {
 289         felem_bytearray b_in, b_out;
 290         felem_to_bin28(b_in, in);
 291         flip_endian(b_out, b_in, sizeof b_out);
 292         return BN_bin2bn(b_out, sizeof b_out, out);
 293         }
 294
 295 /******************************************************************************/
 296 /*                              FIELD OPERATIONS
 297  *
 298  * Field operations, using the internal representation of field elements.
 299  * NB! These operations are specific to our point multiplication and cannot be
 300  * expected to be correct in general - e.g., multiplication with a large scalar
 301  * will cause an overflow.
 302  *
 303  */
 304
 305 /* Sum two field elements: out += in */
 306 static void felem_sum64(fslice out[4], const fslice in[4])
 307         {
 308         out[0] += in[0];
 309         out[1] += in[1];
 310         out[2] += in[2];
 311         out[3] += in[3];
 312         }
 313
 314 /* Subtract field elements: out -= in */
 315 /* Assumes in[i] < 2^57 */
 316 static void felem_diff64(fslice out[4], const fslice in[4])
 317         {
 318         static const uint64_t two58p2 = (((uint64_t) 1) << 58) + (((uint64_t) 1) << 2);
 319         static const uint64_t two58m2 = (((uint64_t) 1) << 58) - (((uint64_t) 1) << 2);
 320         static const uint64_t two58m42m2 = (((uint64_t) 1) << 58) -
 321             (((uint64_t) 1) << 42) - (((uint64_t) 1) << 2);
 322
 323         /* Add 0 mod 2^224-2^96+1 to ensure out > in */
 324         out[0] += two58p2;
 325         out[1] += two58m42m2;
 326         out[2] += two58m2;
 327         out[3] += two58m2;
 328
 329         out[0] -= in[0];
 330         out[1] -= in[1];
 331         out[2] -= in[2];
 332         out[3] -= in[3];
 333         }
 334
 335 /* Subtract in unreduced 128-bit mode: out128 -= in128 */
 336 /* Assumes in[i] < 2^119 */
 337 static void felem_diff128(uint128_t out[7], const uint128_t in[4])
 338         {
 339         static const uint128_t two120 = ((uint128_t) 1) << 120;
 340         static const uint128_t two120m64 = (((uint128_t) 1) << 120) -
 341                 (((uint128_t) 1) << 64);
 342         static const uint128_t two120m104m64 = (((uint128_t) 1) << 120) -
 343                 (((uint128_t) 1) << 104) - (((uint128_t) 1) << 64);
 344
 345         /* Add 0 mod 2^224-2^96+1 to ensure out > in */
 346         out[0] += two120;
 347         out[1] += two120m64;
 348         out[2] += two120m64;
 349         out[3] += two120;
 350         out[4] += two120m104m64;
 351         out[5] += two120m64;
 352         out[6] += two120m64;
 353
 354         out[0] -= in[0];
 355         out[1] -= in[1];
 356         out[2] -= in[2];
 357         out[3] -= in[3];
 358         out[4] -= in[4];
 359         out[5] -= in[5];
 360         out[6] -= in[6];
 361         }
 362
 363 /* Subtract in mixed mode: out128 -= in64 */
 364 /* in[i] < 2^63 */
 365 static void felem_diff_128_64(uint128_t out[7], const fslice in[4])
 366         {
 367         static const uint128_t two64p8 = (((uint128_t) 1) << 64) +
 368                 (((uint128_t) 1) << 8);
 369         static const uint128_t two64m8 = (((uint128_t) 1) << 64) -
 370                 (((uint128_t) 1) << 8);
 371         static const uint128_t two64m48m8 = (((uint128_t) 1) << 64) -
 372                 (((uint128_t) 1) << 48) - (((uint128_t) 1) << 8);
 373
 374         /* Add 0 mod 2^224-2^96+1 to ensure out > in */
 375         out[0] += two64p8;
 376         out[1] += two64m48m8;
 377         out[2] += two64m8;
 378         out[3] += two64m8;
 379
 380         out[0] -= in[0];
 381         out[1] -= in[1];
 382         out[2] -= in[2];
 383         out[3] -= in[3];
 384         }
 385
 386 /* Multiply a field element by a scalar: out64 = out64 * scalar
 387  * The scalars we actually use are small, so results fit without overflow */
 388 static void felem_scalar64(fslice out[4], const fslice scalar)
 389         {
 390         out[0] *= scalar;
 391         out[1] *= scalar;
 392         out[2] *= scalar;
 393         out[3] *= scalar;
 394         }
 395
 396 /* Multiply an unreduced field element by a scalar: out128 = out128 * scalar
 397  * The scalars we actually use are small, so results fit without overflow */
 398 static void felem_scalar128(uint128_t out[7], const uint128_t scalar)
 399         {
 400         out[0] *= scalar;
 401         out[1] *= scalar;
 402         out[2] *= scalar;
 403         out[3] *= scalar;
 404         out[4] *= scalar;
 405         out[5] *= scalar;
 406         out[6] *= scalar;
 407         }
 408
 409 /* Square a field element: out = in^2 */
 410 static void felem_square(uint128_t out[7], const fslice in[4])
 411         {
 412         out[0] = ((uint128_t) in[0]) * in[0];
 413         out[1] = ((uint128_t) in[0]) * in[1] * 2;
 414         out[2] = ((uint128_t) in[0]) * in[2] * 2 + ((uint128_t) in[1]) * in[1];
 415         out[3] = ((uint128_t) in[0]) * in[3] * 2 +
 416                 ((uint128_t) in[1]) * in[2] * 2;
 417         out[4] = ((uint128_t) in[1]) * in[3] * 2 + ((uint128_t) in[2]) * in[2];
 418         out[5] = ((uint128_t) in[2]) * in[3] * 2;
 419         out[6] = ((uint128_t) in[3]) * in[3];
 420         }
 421
 422 /* Multiply two field elements: out = in1 * in2 */
 423 static void felem_mul(uint128_t out[7], const fslice in1[4], const fslice in2[4])
 424         {
 425         out[0] = ((uint128_t) in1[0]) * in2[0];
 426         out[1] = ((uint128_t) in1[0]) * in2[1] + ((uint128_t) in1[1]) * in2[0];
 427         out[2] = ((uint128_t) in1[0]) * in2[2] + ((uint128_t) in1[1]) * in2[1] +
 428                 ((uint128_t) in1[2]) * in2[0];
 429         out[3] = ((uint128_t) in1[0]) * in2[3] + ((uint128_t) in1[1]) * in2[2] +
 430                 ((uint128_t) in1[2]) * in2[1] + ((uint128_t) in1[3]) * in2[0];
 431         out[4] = ((uint128_t) in1[1]) * in2[3] + ((uint128_t) in1[2]) * in2[2] +
 432                 ((uint128_t) in1[3]) * in2[1];
 433         out[5] = ((uint128_t) in1[2]) * in2[3] + ((uint128_t) in1[3]) * in2[2];
 434         out[6] = ((uint128_t) in1[3]) * in2[3];
 435         }
 436
 437 /* Reduce 128-bit coefficients to 64-bit coefficients. Requires in[i] < 2^126,
 438  * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] < 2^57 */
 439 static void felem_reduce(fslice out[4], const uint128_t in[7])
 440         {
 441         static const uint128_t two127p15 = (((uint128_t) 1) << 127) +
 442                 (((uint128_t) 1) << 15);
 443         static const uint128_t two127m71 = (((uint128_t) 1) << 127) -
 444                 (((uint128_t) 1) << 71);
 445         static const uint128_t two127m71m55 = (((uint128_t) 1) << 127) -
 446                 (((uint128_t) 1) << 71) - (((uint128_t) 1) << 55);
 447         uint128_t output[5];
 448
 449         /* Add 0 mod 2^224-2^96+1 to ensure all differences are positive */
 450         output[0] = in[0] + two127p15;
 451         output[1] = in[1] + two127m71m55;
 452         output[2] = in[2] + two127m71;
 453         output[3] = in[3];
 454         output[4] = in[4];
 455
 456         /* Eliminate in[4], in[5], in[6] */
 457         output[4] += in[6] >> 16;
 458         output[3] += (in[6]&0xffff) << 40;
 459         output[2] -= in[6];
 460
 461         output[3] += in[5] >> 16;
 462         output[2] += (in[5]&0xffff) << 40;
 463         output[1] -= in[5];
 464
 465         output[2] += output[4] >> 16;
 466         output[1] += (output[4]&0xffff) << 40;
 467         output[0] -= output[4];
 468         output[4] = 0;
 469
 470         /* Carry 2 -> 3 -> 4 */
 471         output[3] += output[2] >> 56;
 472         output[2] &= 0x00ffffffffffffff;
 473
 474         output[4] += output[3] >> 56;
 475         output[3] &= 0x00ffffffffffffff;
 476
 477         /* Now output[2] < 2^56, output[3] < 2^56 */
 478
 479         /* Eliminate output[4] */
 480         output[2] += output[4] >> 16;
 481         output[1] += (output[4]&0xffff) << 40;
 482         output[0] -= output[4];
 483
 484         /* Carry 0 -> 1 -> 2 -> 3 */
 485         output[1] += output[0] >> 56;
 486         out[0] = output[0] & 0x00ffffffffffffff;
 487
 488         output[2] += output[1] >> 56;
 489         out[1] = output[1] & 0x00ffffffffffffff;
 490         output[3] += output[2] >> 56;
 491         out[2] = output[2] & 0x00ffffffffffffff;
 492
 493         /* out[0] < 2^56, out[1] < 2^56, out[2] < 2^56,
 494          * out[3] < 2^57 (due to final carry) */
 495         out[3] = output[3];
 496         }
 497
 498 /* Reduce to unique minimal representation */
 499 static void felem_contract(fslice out[4], const fslice in[4])
 500         {
 501         static const int64_t two56 = ((uint64_t) 1) << 56;
 502         /* 0 <= in < 2^225 */
 503         /* if in > 2^224 , reduce in = in - 2^224 + 2^96 - 1 */
 504         int64_t tmp[4], a;
 505         tmp[0] = (int64_t) in[0] - (in[3] >> 56);
 506         tmp[1] = (int64_t) in[1] + ((in[3] >> 16) & 0x0000010000000000);
 507         tmp[2] = (int64_t) in[2];
 508         tmp[3] = (int64_t) in[3] & 0x00ffffffffffffff;
 509
 510         /* eliminate negative coefficients */
 511         a = tmp[0] >> 63;
 512         tmp[0] += two56 & a;
 513         tmp[1] -= 1 & a;
 514
 515         a = tmp[1] >> 63;
 516         tmp[1] += two56 & a;
 517         tmp[2] -= 1 & a;
 518
 519         a = tmp[2] >> 63;
 520         tmp[2] += two56 & a;
 521         tmp[3] -= 1 & a;
 522
 523         a = tmp[3] >> 63;
 524         tmp[3] += two56 & a;
 525         tmp[0] += 1 & a;
 526         tmp[1] -= (1 & a) << 40;
 527
 528         /* carry 1 -> 2 -> 3 */
 529         tmp[2] += tmp[1] >> 56;
 530         tmp[1] &= 0x00ffffffffffffff;
 531
 532         tmp[3] += tmp[2] >> 56;
 533         tmp[2] &= 0x00ffffffffffffff;
 534
 535         /* 0 <= in < 2^224 + 2^96 - 1 */
 536         /* if in > 2^224 , reduce in = in - 2^224 + 2^96 - 1 */
 537         tmp[0] -= (tmp[3] >> 56);
 538         tmp[1] += ((tmp[3] >> 16) & 0x0000010000000000);
 539         tmp[3] &= 0x00ffffffffffffff;
 540
 541         /* eliminate negative coefficients */
 542         a = tmp[0] >> 63;
 543         tmp[0] += two56 & a;
 544         tmp[1] -= 1 & a;
 545
 546         a = tmp[1] >> 63;
 547         tmp[1] += two56 & a;
 548         tmp[2] -= 1 & a;
 549
 550         a = tmp[2] >> 63;
 551         tmp[2] += two56 & a;
 552         tmp[3] -= 1 & a;
 553
 554         a = tmp[3] >> 63;
 555         tmp[3] += two56 & a;
 556         tmp[0] += 1 & a;
 557         tmp[1] -= (1 & a) << 40;
 558
 559         /* carry 1 -> 2 -> 3 */
 560         tmp[2] += tmp[1] >> 56;
 561         tmp[1] &= 0x00ffffffffffffff;
 562
 563         tmp[3] += tmp[2] >> 56;
 564         tmp[2] &= 0x00ffffffffffffff;
 565
 566         /* Now 0 <= in < 2^224 */
 567
 568         /* if in > 2^224 - 2^96, reduce */
 569         /* a = 0 iff in > 2^224 - 2^96, i.e.,
 570          * the high 128 bits are all 1 and the lower part is non-zero */
 571         a = (tmp[3] + 1) | (tmp[2] + 1) |
 572                 ((tmp[1] | 0x000000ffffffffff) + 1) |
 573                 ((((tmp[1] & 0xffff) - 1) >> 63) & ((tmp[0] - 1) >> 63));
 574         /* turn a into an all-one mask (if a = 0) or an all-zero mask */
 575         a = ((a & 0x00ffffffffffffff) - 1) >> 63;
 576         /* subtract 2^224 - 2^96 + 1 if a is all-one*/
 577         tmp[3] &= a ^ 0xffffffffffffffff;
 578         tmp[2] &= a ^ 0xffffffffffffffff;
 579         tmp[1] &= (a ^ 0xffffffffffffffff) | 0x000000ffffffffff;
 580         tmp[0] -= 1 & a;
 581         /* eliminate negative coefficients: if tmp[0] is negative, tmp[1] must be
 582          * non-zero, so we only need one step */
 583         a = tmp[0] >> 63;
 584         tmp[0] += two56 & a;
 585         tmp[1] -= 1 & a;
 586
 587         out[0] = tmp[0];
 588         out[1] = tmp[1];
 589         out[2] = tmp[2];
 590         out[3] = tmp[3];
 591         }
 592
 593 /* Zero-check: returns 1 if input is 0, and 0 otherwise.
 594  * We know that field elements are reduced to in < 2^225,
 595  * so we only need to check three cases: 0, 2^224 - 2^96 + 1,
 596  * and 2^225 - 2^97 + 2 */
 597 static fslice felem_is_zero(const fslice in[4])
 598         {
 599         fslice zero, two224m96p1, two225m97p2;
 600
 601         zero = in[0] | in[1] | in[2] | in[3];
 602         zero = (((int64_t)(zero) - 1) >> 63) & 1;
 603         two224m96p1 = (in[0] ^ 1) | (in[1] ^ 0x00ffff0000000000)
 604                 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x00ffffffffffffff);
 605         two224m96p1 = (((int64_t)(two224m96p1) - 1) >> 63) & 1;
 606         two225m97p2 = (in[0] ^ 2) | (in[1] ^ 0x00fffe0000000000)
 607                 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x01ffffffffffffff);
 608         two225m97p2 = (((int64_t)(two225m97p2) - 1) >> 63) & 1;
 609         return (zero | two224m96p1 | two225m97p2);
 610         }
 611
 612 /* Invert a field element */
 613 /* Computation chain copied from djb's code */
 614 static void felem_inv(fslice out[4], const fslice in[4])
 615         {
 616         fslice ftmp[4], ftmp2[4], ftmp3[4], ftmp4[4];
 617         uint128_t tmp[7];
 618         unsigned i;
 619
 620         felem_square(tmp, in); felem_reduce(ftmp, tmp);         /* 2 */
 621         felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp);      /* 2^2 - 1 */
 622         felem_square(tmp, ftmp); felem_reduce(ftmp, tmp);       /* 2^3 - 2 */
 623         felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp);      /* 2^3 - 1 */
 624         felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);      /* 2^4 - 2 */
 625         felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);     /* 2^5 - 4 */
 626         felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);     /* 2^6 - 8 */
 627         felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp);   /* 2^6 - 1 */
 628         felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);      /* 2^7 - 2 */
 629         for (i = 0; i < 5; ++i)                                 /* 2^12 - 2^6 */
 630                 {
 631                 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
 632                 }
 633         felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp2, tmp);  /* 2^12 - 1 */
 634         felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp);     /* 2^13 - 2 */
 635         for (i = 0; i < 11; ++i)                                /* 2^24 - 2^12 */
 636                 {
 637                 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
 638                 }
 639         felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp2, tmp); /* 2^24 - 1 */
 640         felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp);     /* 2^25 - 2 */
 641         for (i = 0; i < 23; ++i)                                /* 2^48 - 2^24 */
 642                 {
 643                 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
 644                 }
 645         felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp3, tmp); /* 2^48 - 1 */
 646         felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp);     /* 2^49 - 2 */
 647         for (i = 0; i < 47; ++i)                                /* 2^96 - 2^48 */
 648                 {
 649                 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
 650                 }
 651         felem_mul(tmp, ftmp3, ftmp4); felem_reduce(ftmp3, tmp); /* 2^96 - 1 */
 652         felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp);     /* 2^97 - 2 */
 653         for (i = 0; i < 23; ++i)                                /* 2^120 - 2^24 */
 654                 {
 655                 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
 656                 }
 657         felem_mul(tmp, ftmp2, ftmp4); felem_reduce(ftmp2, tmp); /* 2^120 - 1 */
 658         for (i = 0; i < 6; ++i)                                 /* 2^126 - 2^6 */
 659                 {
 660                 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
 661                 }
 662         felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp);   /* 2^126 - 1 */
 663         felem_square(tmp, ftmp); felem_reduce(ftmp, tmp);       /* 2^127 - 2 */
 664         felem_mul(tmp, ftmp, in); felem_reduce(ftmp, tmp);      /* 2^127 - 1 */
 665         for (i = 0; i < 97; ++i)                                /* 2^224 - 2^97 */
 666                 {
 667                 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp);
 668                 }
 669         felem_mul(tmp, ftmp, ftmp3); felem_reduce(out, tmp);    /* 2^224 - 2^96 - 1 */
 670         }
 671
 672 /* Copy in constant time:
 673  * if icopy == 1, copy in to out,
 674  * if icopy == 0, copy out to itself. */
 675 static void
 676 copy_conditional(fslice *out, const fslice *in, unsigned len, fslice icopy)
 677         {
 678         unsigned i;
 679         /* icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one */
 680         const fslice copy = -icopy;
 681         for (i = 0; i < len; ++i)
 682                 {
 683                 const fslice tmp = copy & (in[i] ^ out[i]);
 684                 out[i] ^= tmp;
 685                 }
 686         }
 687
 688 /* Copy in constant time:
 689  * if isel == 1, copy in2 to out,
 690  * if isel == 0, copy in1 to out. */
 691 static void select_conditional(fslice *out, const fslice *in1, const fslice *in2,
 692         unsigned len, fslice isel)
 693         {
 694         unsigned i;
 695         /* isel is a (64-bit) 0 or 1, so sel is either all-zero or all-one */
 696         const fslice sel = -isel;
 697         for (i = 0; i < len; ++i)
 698                 {
 699                 const fslice tmp = sel & (in1[i] ^ in2[i]);
 700                 out[i] = in1[i] ^ tmp;
 701                 }
 702 }
 703
 704 /******************************************************************************/
 705 /*                       ELLIPTIC CURVE POINT OPERATIONS
 706  *
 707  * Points are represented in Jacobian projective coordinates:
 708  * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
 709  * or to the point at infinity if Z == 0.
 710  *
 711  */
 712
 713 /* Double an elliptic curve point:
 714  * (X', Y', Z') = 2 * (X, Y, Z), where
 715  * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
 716  * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^2
 717  * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
 718  * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
 719  * while x_out == y_in is not (maybe this works, but it's not tested). */
 720 static void
 721 point_double(fslice x_out[4], fslice y_out[4], fslice z_out[4],
 722              const fslice x_in[4], const fslice y_in[4], const fslice z_in[4])
 723         {
 724         uint128_t tmp[7], tmp2[7];
 725         fslice delta[4];
 726         fslice gamma[4];
 727         fslice beta[4];
 728         fslice alpha[4];
 729         fslice ftmp[4], ftmp2[4];
 730         memcpy(ftmp, x_in, 4 * sizeof(fslice));
 731         memcpy(ftmp2, x_in, 4 * sizeof(fslice));
 732
 733         /* delta = z^2 */
 734         felem_square(tmp, z_in);
 735         felem_reduce(delta, tmp);
 736
 737         /* gamma = y^2 */
 738         felem_square(tmp, y_in);
 739         felem_reduce(gamma, tmp);
 740
 741         /* beta = x*gamma */
 742         felem_mul(tmp, x_in, gamma);
 743         felem_reduce(beta, tmp);
 744
 745         /* alpha = 3*(x-delta)*(x+delta) */
 746         felem_diff64(ftmp, delta);
 747         /* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
 748         felem_sum64(ftmp2, delta);
 749         /* ftmp2[i] < 2^57 + 2^57 = 2^58 */
 750         felem_scalar64(ftmp2, 3);
 751         /* ftmp2[i] < 3 * 2^58 < 2^60 */
 752         felem_mul(tmp, ftmp, ftmp2);
 753         /* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
 754         felem_reduce(alpha, tmp);
 755
 756         /* x' = alpha^2 - 8*beta */
 757         felem_square(tmp, alpha);
 758         /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 759         memcpy(ftmp, beta, 4 * sizeof(fslice));
 760         felem_scalar64(ftmp, 8);
 761         /* ftmp[i] < 8 * 2^57 = 2^60 */
 762         felem_diff_128_64(tmp, ftmp);
 763         /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
 764         felem_reduce(x_out, tmp);
 765
 766         /* z' = (y + z)^2 - gamma - delta */
 767         felem_sum64(delta, gamma);
 768         /* delta[i] < 2^57 + 2^57 = 2^58 */
 769         memcpy(ftmp, y_in, 4 * sizeof(fslice));
 770         felem_sum64(ftmp, z_in);
 771         /* ftmp[i] < 2^57 + 2^57 = 2^58 */
 772         felem_square(tmp, ftmp);
 773         /* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
 774         felem_diff_128_64(tmp, delta);
 775         /* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
 776         felem_reduce(z_out, tmp);
 777
 778         /* y' = alpha*(4*beta - x') - 8*gamma^2 */
 779         felem_scalar64(beta, 4);
 780         /* beta[i] < 4 * 2^57 = 2^59 */
 781         felem_diff64(beta, x_out);
 782         /* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
 783         felem_mul(tmp, alpha, beta);
 784         /* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
 785         felem_square(tmp2, gamma);
 786         /* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
 787         felem_scalar128(tmp2, 8);
 788         /* tmp2[i] < 8 * 2^116 = 2^119 */
 789         felem_diff128(tmp, tmp2);
 790         /* tmp[i] < 2^119 + 2^120 < 2^121 */
 791         felem_reduce(y_out, tmp);
 792         }
 793
 794 /* Add two elliptic curve points:
 795  * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
 796  * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
 797  * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
 798  * Y_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1) * (Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2 - X_3) -
 799  *        Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
 800  * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2) */
 801
 802 /* This function is not entirely constant-time:
 803  * it includes a branch for checking whether the two input points are equal,
 804  * (while not equal to the point at infinity).
 805  * This case never happens during single point multiplication,
 806  * so there is no timing leak for ECDH or ECDSA signing. */
 807 static void point_add(fslice x3[4], fslice y3[4], fslice z3[4],
 808         const fslice x1[4], const fslice y1[4], const fslice z1[4],
 809         const fslice x2[4], const fslice y2[4], const fslice z2[4])
 810         {
 811         fslice ftmp[4], ftmp2[4], ftmp3[4], ftmp4[4], ftmp5[4];
 812         uint128_t tmp[7], tmp2[7];
 813         fslice z1_is_zero, z2_is_zero, x_equal, y_equal;
 814
 815         /* ftmp = z1^2 */
 816         felem_square(tmp, z1);
 817         felem_reduce(ftmp, tmp);
 818
 819         /* ftmp2 = z2^2 */
 820         felem_square(tmp, z2);
 821         felem_reduce(ftmp2, tmp);
 822
 823         /* ftmp3 = z1^3 */
 824         felem_mul(tmp, ftmp, z1);
 825         felem_reduce(ftmp3, tmp);
 826
 827         /* ftmp4 = z2^3 */
 828         felem_mul(tmp, ftmp2, z2);
 829         felem_reduce(ftmp4, tmp);
 830
 831         /* ftmp3 = z1^3*y2 */
 832         felem_mul(tmp, ftmp3, y2);
 833         /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 834
 835         /* ftmp4 = z2^3*y1 */
 836         felem_mul(tmp2, ftmp4, y1);
 837         felem_reduce(ftmp4, tmp2);
 838
 839         /* ftmp3 = z1^3*y2 - z2^3*y1 */
 840         felem_diff_128_64(tmp, ftmp4);
 841         /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
 842         felem_reduce(ftmp3, tmp);
 843
 844         /* ftmp = z1^2*x2 */
 845         felem_mul(tmp, ftmp, x2);
 846         /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 847
 848         /* ftmp2 =z2^2*x1 */
 849         felem_mul(tmp2, ftmp2, x1);
 850         felem_reduce(ftmp2, tmp2);
 851
 852         /* ftmp = z1^2*x2 - z2^2*x1 */
 853         felem_diff128(tmp, tmp2);
 854         /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
 855         felem_reduce(ftmp, tmp);
 856
 857         /* the formulae are incorrect if the points are equal
 858          * so we check for this and do doubling if this happens */
 859         x_equal = felem_is_zero(ftmp);
 860         y_equal = felem_is_zero(ftmp3);
 861         z1_is_zero = felem_is_zero(z1);
 862         z2_is_zero = felem_is_zero(z2);
 863         /* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
 864         if (x_equal && y_equal && !z1_is_zero && !z2_is_zero)
 865                 {
 866                 point_double(x3, y3, z3, x1, y1, z1);
 867                 return;
 868                 }
 869
 870         /* ftmp5 = z1*z2 */
 871         felem_mul(tmp, z1, z2);
 872         felem_reduce(ftmp5, tmp);
 873
 874         /* z3 = (z1^2*x2 - z2^2*x1)*(z1*z2) */
 875         felem_mul(tmp, ftmp, ftmp5);
 876         felem_reduce(z3, tmp);
 877
 878         /* ftmp = (z1^2*x2 - z2^2*x1)^2 */
 879         memcpy(ftmp5, ftmp, 4 * sizeof(fslice));
 880         felem_square(tmp, ftmp);
 881         felem_reduce(ftmp, tmp);
 882
 883         /* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
 884         felem_mul(tmp, ftmp, ftmp5);
 885         felem_reduce(ftmp5, tmp);
 886
 887         /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
 888         felem_mul(tmp, ftmp2, ftmp);
 889         felem_reduce(ftmp2, tmp);
 890
 891         /* ftmp4 = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
 892         felem_mul(tmp, ftmp4, ftmp5);
 893         /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 894
 895         /* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
 896         felem_square(tmp2, ftmp3);
 897         /* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
 898
 899         /* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
 900         felem_diff_128_64(tmp2, ftmp5);
 901         /* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
 902
 903         /* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
 904         memcpy(ftmp5, ftmp2, 4 * sizeof(fslice));
 905         felem_scalar64(ftmp5, 2);
 906         /* ftmp5[i] < 2 * 2^57 = 2^58 */
 907
 908         /* x3 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
 909            2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
 910         felem_diff_128_64(tmp2, ftmp5);
 911         /* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
 912         felem_reduce(x3, tmp2);
 913
 914         /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x3 */
 915         felem_diff64(ftmp2, x3);
 916         /* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
 917
 918         /* tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x3) */
 919         felem_mul(tmp2, ftmp3, ftmp2);
 920         /* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
 921
 922         /* y3 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x3) -
 923            z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
 924         felem_diff128(tmp2, tmp);
 925         /* tmp2[i] < 2^118 + 2^120 < 2^121 */
 926         felem_reduce(y3, tmp2);
 927
 928         /* the result (x3, y3, z3) is incorrect if one of the inputs is the
 929          * point at infinity, so we need to check for this separately */
 930
 931         /* if point 1 is at infinity, copy point 2 to output, and vice versa */
 932         copy_conditional(x3, x2, 4, z1_is_zero);
 933         copy_conditional(x3, x1, 4, z2_is_zero);
 934         copy_conditional(y3, y2, 4, z1_is_zero);
 935         copy_conditional(y3, y1, 4, z2_is_zero);
 936         copy_conditional(z3, z2, 4, z1_is_zero);
 937         copy_conditional(z3, z1, 4, z2_is_zero);
 938         }
 939
 940 /* Select a point from an array of 16 precomputed point multiples,
 941  * in constant time: for bits = {b_0, b_1, b_2, b_3}, return the point
 942  * pre_comp[8*b_3 + 4*b_2 + 2*b_1 + b_0] */
 943 static void select_point(const fslice bits[4], const fslice pre_comp[16][3][4],
 944         fslice out[12])
 945         {
 946         fslice tmp[5][12];
 947         select_conditional(tmp[0], pre_comp[7][0], pre_comp[15][0], 12, bits[3]);
 948         select_conditional(tmp[1], pre_comp[3][0], pre_comp[11][0], 12, bits[3]);
 949         select_conditional(tmp[2], tmp[1], tmp[0], 12, bits[2]);
 950         select_conditional(tmp[0], pre_comp[5][0], pre_comp[13][0], 12, bits[3]);
 951         select_conditional(tmp[1], pre_comp[1][0], pre_comp[9][0], 12, bits[3]);
 952         select_conditional(tmp[3], tmp[1], tmp[0], 12, bits[2]);
 953         select_conditional(tmp[4], tmp[3], tmp[2], 12, bits[1]);
 954         select_conditional(tmp[0], pre_comp[6][0], pre_comp[14][0], 12, bits[3]);
 955         select_conditional(tmp[1], pre_comp[2][0], pre_comp[10][0], 12, bits[3]);
 956         select_conditional(tmp[2], tmp[1], tmp[0], 12, bits[2]);
 957         select_conditional(tmp[0], pre_comp[4][0], pre_comp[12][0], 12, bits[3]);
 958         select_conditional(tmp[1], pre_comp[0][0], pre_comp[8][0], 12, bits[3]);
 959         select_conditional(tmp[3], tmp[1], tmp[0], 12, bits[2]);
 960         select_conditional(tmp[1], tmp[3], tmp[2], 12, bits[1]);
 961         select_conditional(out, tmp[1], tmp[4], 12, bits[0]);
 962         }
 963
 964 /* Interleaved point multiplication using precomputed point multiples:
 965  * The small point multiples 0*P, 1*P, ..., 15*P are in pre_comp[],
 966  * the scalars in scalars[]. If g_scalar is non-NULL, we also add this multiple
 967  * of the generator, using certain (large) precomputed multiples in g_pre_comp.
 968  * Output point (X, Y, Z) is stored in x_out, y_out, z_out */
 969 static void batch_mul(fslice x_out[4], fslice y_out[4], fslice z_out[4],
 970         const felem_bytearray scalars[], const unsigned num_points, const u8 *g_scalar,
 971         const fslice pre_comp[][16][3][4], const fslice g_pre_comp[16][3][4])
 972         {
 973         unsigned i, j, num;
 974         unsigned gen_mul = (g_scalar != NULL);
 975         fslice nq[12], nqt[12], tmp[12];
 976         fslice bits[4];
 977         u8 byte;
 978
 979         /* set nq to the point at infinity */
 980         memset(nq, 0, 12 * sizeof(fslice));
 981
 982         /* Loop over all scalars msb-to-lsb, 4 bits at a time: for each nibble,
 983          * double 4 times, then add the precomputed point multiples.
 984          * If we are also adding multiples of the generator, then interleave
 985          * these additions with the last 56 doublings. */
 986         for (i = (num_points ? 28 : 7); i > 0; --i)
 987                 {
 988                 for (j = 0; j < 8; ++j)
 989                         {
 990                         /* double once */
 991                         point_double(nq, nq+4, nq+8, nq, nq+4, nq+8);
 992                         /* add multiples of the generator */
 993                         if ((gen_mul) && (i <= 7))
 994                                 {
 995                                 bits[3] = (g_scalar[i+20] >> (7-j)) & 1;
 996                                 bits[2] = (g_scalar[i+13] >> (7-j)) & 1;
 997                                 bits[1] = (g_scalar[i+6] >> (7-j)) & 1;
 998                                 bits[0] = (g_scalar[i-1] >> (7-j)) & 1;
 999                                 /* select the point to add, in constant time */
1000                                 select_point(bits, g_pre_comp, tmp);
1001                                 memcpy(nqt, nq, 12 * sizeof(fslice));
1002                                 point_add(nq, nq+4, nq+8, nqt, nqt+4, nqt+8,
1003                                         tmp, tmp+4, tmp+8);
1004                                 }
1005                         /* do an addition after every 4 doublings */
1006                         if (j % 4 == 3)
1007                                 {
1008                                 /* loop over all scalars */
1009                                 for (num = 0; num < num_points; ++num)
1010                                         {
1011                                         byte = scalars[num][i-1];
1012                                         bits[3] = (byte >> (10-j)) & 1;
1013                                         bits[2] = (byte >> (9-j)) & 1;
1014                                         bits[1] = (byte >> (8-j)) & 1;
1015                                         bits[0] = (byte >> (7-j)) & 1;
1016                                         /* select the point to add */
1017                                         select_point(bits,
1018                                                 pre_comp[num], tmp);
1019                                         memcpy(nqt, nq, 12 * sizeof(fslice));
1020                                         point_add(nq, nq+4, nq+8, nqt, nqt+4,
1021                                                 nqt+8, tmp, tmp+4, tmp+8);
1022                                         }
1023                                 }
1024                         }
1025                 }
1026         memcpy(x_out, nq, 4 * sizeof(fslice));
1027         memcpy(y_out, nq+4, 4 * sizeof(fslice));
1028         memcpy(z_out, nq+8, 4 * sizeof(fslice));
1029         }
1030
1031 /******************************************************************************/
1032 /*                     FUNCTIONS TO MANAGE PRECOMPUTATION
1033  */
1034
1035 static NISTP224_PRE_COMP *nistp224_pre_comp_new()
1036         {
1037         NISTP224_PRE_COMP *ret = NULL;
1038         ret = (NISTP224_PRE_COMP *)OPENSSL_malloc(sizeof(NISTP224_PRE_COMP));
1039         if (!ret)
1040                 {
1041                 ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1042                 return ret;
1043                 }
1044         memset(ret->g_pre_comp, 0, sizeof(ret->g_pre_comp));
1045         ret->references = 1;
1046         return ret;
1047         }
1048
1049 static void *nistp224_pre_comp_dup(void *src_)
1050         {
1051         NISTP224_PRE_COMP *src = src_;
1052
1053         /* no need to actually copy, these objects never change! */
1054         CRYPTO_add(&src->references, 1, CRYPTO_LOCK_EC_PRE_COMP);
1055
1056         return src_;
1057         }
1058
1059 static void nistp224_pre_comp_free(void *pre_)
1060         {
1061         int i;
1062         NISTP224_PRE_COMP *pre = pre_;
1063
1064         if (!pre)
1065                 return;
1066
1067         i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1068         if (i > 0)
1069                 return;
1070
1071         OPENSSL_free(pre);
1072         }
1073
1074 static void nistp224_pre_comp_clear_free(void *pre_)
1075         {
1076         int i;
1077         NISTP224_PRE_COMP *pre = pre_;
1078
1079         if (!pre)
1080                 return;
1081
1082         i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1083         if (i > 0)
1084                 return;
1085
1086         OPENSSL_cleanse(pre, sizeof *pre);
1087         OPENSSL_free(pre);
1088         }
1089
1090 /******************************************************************************/
1091 /*                         OPENSSL EC_METHOD FUNCTIONS
1092  */
1093
1094 int ec_GFp_nistp224_group_init(EC_GROUP *group)
1095         {
1096         int ret;
1097         ret = ec_GFp_simple_group_init(group);
1098         group->a_is_minus3 = 1;
1099         return ret;
1100         }
1101
1102 int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1103         const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
1104         {
1105         int ret = 0;
1106         BN_CTX *new_ctx = NULL;
1107         BIGNUM *curve_p, *curve_a, *curve_b;
1108
1109         if (ctx == NULL)
1110                 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1111         BN_CTX_start(ctx);
1112         if (((curve_p = BN_CTX_get(ctx)) == NULL) ||
1113                 ((curve_a = BN_CTX_get(ctx)) == NULL) ||
1114                 ((curve_b = BN_CTX_get(ctx)) == NULL)) goto err;
1115         BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
1116         BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
1117         BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
1118         if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) ||
1119                 (BN_cmp(curve_b, b)))
1120                 {
1121                 ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
1122                         EC_R_WRONG_CURVE_PARAMETERS);
1123                 goto err;
1124                 }
1125         group->field_mod_func = BN_nist_mod_224;
1126         ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1127 err:
1128         BN_CTX_end(ctx);
1129         if (new_ctx != NULL)
1130                 BN_CTX_free(new_ctx);
1131         return ret;
1132         }
1133
1134 /* Takes the Jacobian coordinates (X, Y, Z) of a point and returns
1135  * (X', Y') = (X/Z^2, Y/Z^3) */
1136 int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
1137         const EC_POINT *point, BIGNUM *x, BIGNUM *y, BN_CTX *ctx)
1138         {
1139         fslice z1[4], z2[4], x_in[4], y_in[4], x_out[4], y_out[4];
1140         uint128_t tmp[7];
1141
1142         if (EC_POINT_is_at_infinity(group, point))
1143                 {
1144                 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1145                         EC_R_POINT_AT_INFINITY);
1146                 return 0;
1147                 }
1148         if ((!BN_to_felem(x_in, &point->X)) || (!BN_to_felem(y_in, &point->Y)) ||
1149                 (!BN_to_felem(z1, &point->Z))) return 0;
1150         felem_inv(z2, z1);
1151         felem_square(tmp, z2); felem_reduce(z1, tmp);
1152         felem_mul(tmp, x_in, z1); felem_reduce(x_in, tmp);
1153         felem_contract(x_out, x_in);
1154         if (x != NULL)
1155                 {
1156                 if (!felem_to_BN(x, x_out)) {
1157                 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1158                         ERR_R_BN_LIB);
1159                 return 0;
1160                 }
1161                 }
1162         felem_mul(tmp, z1, z2); felem_reduce(z1, tmp);
1163         felem_mul(tmp, y_in, z1); felem_reduce(y_in, tmp);
1164         felem_contract(y_out, y_in);
1165         if (y != NULL)
1166                 {
1167                 if (!felem_to_BN(y, y_out)) {
1168                 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1169                         ERR_R_BN_LIB);
1170                 return 0;
1171                 }
1172                 }
1173         return 1;
1174         }
1175
1176 /* Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL values
1177  * Result is stored in r (r can equal one of the inputs). */
1178 int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
1179         const BIGNUM *scalar, size_t num, const EC_POINT *points[],
1180         const BIGNUM *scalars[], BN_CTX *ctx)
1181         {
1182         int ret = 0;
1183         int i, j;
1184         BN_CTX *new_ctx = NULL;
1185         BIGNUM *x, *y, *z, *tmp_scalar;
1186         felem_bytearray g_secret;
1187         felem_bytearray *secrets = NULL;
1188         fslice (*pre_comp)[16][3][4] = NULL;
1189         felem_bytearray tmp;
1190         unsigned num_bytes;
1191         int have_pre_comp = 0;
1192         size_t num_points = num;
1193         fslice x_in[4], y_in[4], z_in[4], x_out[4], y_out[4], z_out[4];
1194         NISTP224_PRE_COMP *pre = NULL;
1195         fslice (*g_pre_comp)[3][4] = NULL;
1196         EC_POINT *generator = NULL;
1197         const EC_POINT *p = NULL;
1198         const BIGNUM *p_scalar = NULL;
1199
1200         if (ctx == NULL)
1201                 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1202         BN_CTX_start(ctx);
1203         if (((x = BN_CTX_get(ctx)) == NULL) ||
1204                 ((y = BN_CTX_get(ctx)) == NULL) ||
1205                 ((z = BN_CTX_get(ctx)) == NULL) ||
1206                 ((tmp_scalar = BN_CTX_get(ctx)) == NULL))
1207                 goto err;
1208
1209         if (scalar != NULL)
1210                 {
1211                 pre = EC_EX_DATA_get_data(group->extra_data,
1212                         nistp224_pre_comp_dup, nistp224_pre_comp_free,
1213                         nistp224_pre_comp_clear_free);
1214                 if (pre)
1215                         /* we have precomputation, try to use it */
1216                         g_pre_comp = pre->g_pre_comp;
1217                 else
1218                         /* try to use the standard precomputation */
1219                         g_pre_comp = (fslice (*)[3][4]) gmul;
1220                 generator = EC_POINT_new(group);
1221                 if (generator == NULL)
1222                         goto err;
1223                 /* get the generator from precomputation */
1224                 if (!felem_to_BN(x, g_pre_comp[1][0]) ||
1225                         !felem_to_BN(y, g_pre_comp[1][1]) ||
1226                         !felem_to_BN(z, g_pre_comp[1][2]))
1227                         {
1228                         ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1229                         goto err;
1230                         }
1231                 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1232                                 generator, x, y, z, ctx))
1233                         goto err;
1234                 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1235                         /* precomputation matches generator */
1236                         have_pre_comp = 1;
1237                 else
1238                         /* we don't have valid precomputation:
1239                          * treat the generator as a random point */
1240                         num_points = num_points + 1;
1241                 }
1242         secrets = OPENSSL_malloc(num_points * sizeof(felem_bytearray));
1243         pre_comp = OPENSSL_malloc(num_points * 16 * 3 * 4 * sizeof(fslice));
1244
1245         if ((num_points) && ((secrets == NULL) || (pre_comp == NULL)))
1246                 {
1247                 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1248                 goto err;
1249                 }
1250
1251         /* we treat NULL scalars as 0, and NULL points as points at infinity,
1252          * i.e., they contribute nothing to the linear combination */
1253         memset(secrets, 0, num_points * sizeof(felem_bytearray));
1254         memset(pre_comp, 0, num_points * 16 * 3 * 4 * sizeof(fslice));
1255         for (i = 0; i < num_points; ++i)
1256                 {
1257                 if (i == num)
1258                         /* the generator */
1259                         {
1260                         p = EC_GROUP_get0_generator(group);
1261                         p_scalar = scalar;
1262                         }
1263                 else
1264                         /* the i^th point */
1265                         {
1266                         p = points[i];
1267                         p_scalar = scalars[i];
1268                         }
1269                 if ((p_scalar != NULL) && (p != NULL))
1270                         {
1271                         num_bytes = BN_num_bytes(p_scalar);
1272                         /* reduce scalar to 0 <= scalar < 2^224 */
1273                         if ((num_bytes > sizeof(felem_bytearray)) || (BN_is_negative(p_scalar)))
1274                                 {
1275                                 /* this is an unusual input, and we don't guarantee
1276                                  * constant-timeness */
1277                                 if (!BN_nnmod(tmp_scalar, p_scalar, &group->order, ctx))
1278                                         {
1279                                         ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1280                                         goto err;
1281                                         }
1282                                 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1283                                 }
1284                         else
1285                                 BN_bn2bin(p_scalar, tmp);
1286                         flip_endian(secrets[i], tmp, num_bytes);
1287                         /* precompute multiples */
1288                         if ((!BN_to_felem(x_out, &p->X)) ||
1289                                 (!BN_to_felem(y_out, &p->Y)) ||
1290                                 (!BN_to_felem(z_out, &p->Z))) goto err;
1291                         memcpy(pre_comp[i][1][0], x_out, 4 * sizeof(fslice));
1292                         memcpy(pre_comp[i][1][1], y_out, 4 * sizeof(fslice));
1293                         memcpy(pre_comp[i][1][2], z_out, 4 * sizeof(fslice));
1294                         for (j = 1; j < 8; ++j)
1295                                 {
1296                                 point_double(pre_comp[i][2*j][0],
1297                                         pre_comp[i][2*j][1],
1298                                         pre_comp[i][2*j][2],
1299                                         pre_comp[i][j][0],
1300                                         pre_comp[i][j][1],
1301                                         pre_comp[i][j][2]);
1302                                 point_add(pre_comp[i][2*j+1][0],
1303                                         pre_comp[i][2*j+1][1],
1304                                         pre_comp[i][2*j+1][2],
1305                                         pre_comp[i][1][0],
1306                                         pre_comp[i][1][1],
1307                                         pre_comp[i][1][2],
1308                                         pre_comp[i][2*j][0],
1309                                         pre_comp[i][2*j][1],
1310                                         pre_comp[i][2*j][2]);
1311                                 }
1312                         }
1313                 }
1314
1315         /* the scalar for the generator */
1316         if ((scalar != NULL) && (have_pre_comp))
1317                 {
1318                 memset(g_secret, 0, sizeof g_secret);
1319                 num_bytes = BN_num_bytes(scalar);
1320                 /* reduce scalar to 0 <= scalar < 2^224 */
1321                 if ((num_bytes > sizeof(felem_bytearray)) || (BN_is_negative(scalar)))
1322                         {
1323                         /* this is an unusual input, and we don't guarantee
1324                          * constant-timeness */
1325                         if (!BN_nnmod(tmp_scalar, scalar, &group->order, ctx))
1326                                 {
1327                                 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1328                                 goto err;
1329                                 }
1330                         num_bytes = BN_bn2bin(tmp_scalar, tmp);
1331                         }
1332                 else
1333                         BN_bn2bin(scalar, tmp);
1334                 flip_endian(g_secret, tmp, num_bytes);
1335                 /* do the multiplication with generator precomputation*/
1336                 batch_mul(x_out, y_out, z_out,
1337                         (const felem_bytearray (*)) secrets, num_points,
1338                         g_secret, (const fslice (*)[16][3][4]) pre_comp,
1339                         (const fslice (*)[3][4]) g_pre_comp);
1340                 }
1341         else
1342                 /* do the multiplication without generator precomputation */
1343                 batch_mul(x_out, y_out, z_out,
1344                         (const felem_bytearray (*)) secrets, num_points,
1345                         NULL, (const fslice (*)[16][3][4]) pre_comp, NULL);
1346         /* reduce the output to its unique minimal representation */
1347         felem_contract(x_in, x_out);
1348         felem_contract(y_in, y_out);
1349         felem_contract(z_in, z_out);
1350         if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
1351                 (!felem_to_BN(z, z_in)))
1352                 {
1353                 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1354                 goto err;
1355                 }
1356         ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
1357
1358 err:
1359         BN_CTX_end(ctx);
1360         if (generator != NULL)
1361                 EC_POINT_free(generator);
1362         if (new_ctx != NULL)
1363                 BN_CTX_free(new_ctx);
1364         if (secrets != NULL)
1365                 OPENSSL_free(secrets);
1366         if (pre_comp != NULL)
1367                 OPENSSL_free(pre_comp);
1368         return ret;
1369         }
1370
1371 int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
1372         {
1373         int ret = 0;
1374         NISTP224_PRE_COMP *pre = NULL;
1375         int i, j;
1376         BN_CTX *new_ctx = NULL;
1377         BIGNUM *x, *y;
1378         EC_POINT *generator = NULL;
1379
1380         /* throw away old precomputation */
1381         EC_EX_DATA_free_data(&group->extra_data, nistp224_pre_comp_dup,
1382                 nistp224_pre_comp_free, nistp224_pre_comp_clear_free);
1383         if (ctx == NULL)
1384                 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1385         BN_CTX_start(ctx);
1386         if (((x = BN_CTX_get(ctx)) == NULL) ||
1387                 ((y = BN_CTX_get(ctx)) == NULL))
1388                 goto err;
1389         /* get the generator */
1390         if (group->generator == NULL) goto err;
1391         generator = EC_POINT_new(group);
1392         if (generator == NULL)
1393                 goto err;
1394         BN_bin2bn(nistp224_curve_params[3], sizeof (felem_bytearray), x);
1395         BN_bin2bn(nistp224_curve_params[4], sizeof (felem_bytearray), y);
1396         if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
1397                 goto err;
1398         if ((pre = nistp224_pre_comp_new()) == NULL)
1399                 goto err;
1400         /* if the generator is the standard one, use built-in precomputation */
1401         if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1402                 {
1403                 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
1404                 ret = 1;
1405                 goto err;
1406                 }
1407         if ((!BN_to_felem(pre->g_pre_comp[1][0], &group->generator->X)) ||
1408                 (!BN_to_felem(pre->g_pre_comp[1][1], &group->generator->Y)) ||
1409                 (!BN_to_felem(pre->g_pre_comp[1][2], &group->generator->Z)))
1410                 goto err;
1411         /* compute 2^56*G, 2^112*G, 2^168*G */
1412         for (i = 1; i < 5; ++i)
1413                 {
1414                 point_double(pre->g_pre_comp[2*i][0], pre->g_pre_comp[2*i][1],
1415                         pre->g_pre_comp[2*i][2], pre->g_pre_comp[i][0],
1416                         pre->g_pre_comp[i][1], pre->g_pre_comp[i][2]);
1417                 for (j = 0; j < 55; ++j)
1418                         {
1419                         point_double(pre->g_pre_comp[2*i][0],
1420                                 pre->g_pre_comp[2*i][1],
1421                                 pre->g_pre_comp[2*i][2],
1422                                 pre->g_pre_comp[2*i][0],
1423                                 pre->g_pre_comp[2*i][1],
1424                                 pre->g_pre_comp[2*i][2]);
1425                         }
1426                 }
1427         /* g_pre_comp[0] is the point at infinity */
1428         memset(pre->g_pre_comp[0], 0, sizeof(pre->g_pre_comp[0]));
1429         /* the remaining multiples */
1430         /* 2^56*G + 2^112*G */
1431         point_add(pre->g_pre_comp[6][0], pre->g_pre_comp[6][1],
1432                 pre->g_pre_comp[6][2], pre->g_pre_comp[4][0],
1433                 pre->g_pre_comp[4][1], pre->g_pre_comp[4][2],
1434                 pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
1435                 pre->g_pre_comp[2][2]);
1436         /* 2^56*G + 2^168*G */
1437         point_add(pre->g_pre_comp[10][0], pre->g_pre_comp[10][1],
1438                 pre->g_pre_comp[10][2], pre->g_pre_comp[8][0],
1439                 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
1440                 pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
1441                 pre->g_pre_comp[2][2]);
1442         /* 2^112*G + 2^168*G */
1443         point_add(pre->g_pre_comp[12][0], pre->g_pre_comp[12][1],
1444                 pre->g_pre_comp[12][2], pre->g_pre_comp[8][0],
1445                 pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
1446                 pre->g_pre_comp[4][0], pre->g_pre_comp[4][1],
1447                 pre->g_pre_comp[4][2]);
1448         /* 2^56*G + 2^112*G + 2^168*G */
1449         point_add(pre->g_pre_comp[14][0], pre->g_pre_comp[14][1],
1450                 pre->g_pre_comp[14][2], pre->g_pre_comp[12][0],
1451                 pre->g_pre_comp[12][1], pre->g_pre_comp[12][2],
1452                 pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
1453                 pre->g_pre_comp[2][2]);
1454         for (i = 1; i < 8; ++i)
1455                 {
1456                 /* odd multiples: add G */
1457                 point_add(pre->g_pre_comp[2*i+1][0], pre->g_pre_comp[2*i+1][1],
1458                         pre->g_pre_comp[2*i+1][2], pre->g_pre_comp[2*i][0],
1459                         pre->g_pre_comp[2*i][1], pre->g_pre_comp[2*i][2],
1460                         pre->g_pre_comp[1][0], pre->g_pre_comp[1][1],
1461                         pre->g_pre_comp[1][2]);
1462                 }
1463
1464         if (!EC_EX_DATA_set_data(&group->extra_data, pre, nistp224_pre_comp_dup,
1465                         nistp224_pre_comp_free, nistp224_pre_comp_clear_free))
1466                 goto err;
1467         ret = 1;
1468         pre = NULL;
1469  err:
1470         BN_CTX_end(ctx);
1471         if (generator != NULL)
1472                 EC_POINT_free(generator);
1473         if (new_ctx != NULL)
1474                 BN_CTX_free(new_ctx);
1475         if (pre)
1476                 nistp224_pre_comp_free(pre);
1477         return ret;
1478         }
1479
1480 int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
1481         {
1482         if (EC_EX_DATA_get_data(group->extra_data, nistp224_pre_comp_dup,
1483                         nistp224_pre_comp_free, nistp224_pre_comp_clear_free)
1484                 != NULL)
1485                 return 1;
1486         else
1487                 return 0;
1488         }
1489
1490 #else
1491 static void *dummy=&dummy;
1492 #endif