crypto/ec/ecp_nistp224.c

   1 /*
   2  * Copyright 2010-2018 The OpenSSL Project Authors. All Rights Reserved.
   3  *
   4  * Licensed under the Apache License 2.0 (the "License").  You may not use
   5  * this file except in compliance with the License.  You can obtain a copy
   6  * in the file LICENSE in the source distribution or at
   7  * https://www.openssl.org/source/license.html
   8  */
   9
  10 /* Copyright 2011 Google Inc.
  11  *
  12  * Licensed under the Apache License, Version 2.0 (the "License");
  13  *
  14  * you may not use this file except in compliance with the License.
  15  * You may obtain a copy of the License at
  16  *
  17  *     http://www.apache.org/licenses/LICENSE-2.0
  18  *
  19  *  Unless required by applicable law or agreed to in writing, software
  20  *  distributed under the License is distributed on an "AS IS" BASIS,
  21  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  22  *  See the License for the specific language governing permissions and
  23  *  limitations under the License.
  24  */
  25
  26 /*
  27  * A 64-bit implementation of the NIST P-224 elliptic curve point multiplication
  28  *
  29  * Inspired by Daniel J. Bernstein's public domain nistp224 implementation
  30  * and Adam Langley's public domain 64-bit C implementation of curve25519
  31  */
  32
  33 #include <openssl/opensslconf.h>
  34 #ifdef OPENSSL_NO_EC_NISTP_64_GCC_128
  35 NON_EMPTY_TRANSLATION_UNIT
  36 #else
  37
  38 # include <stdint.h>
  39 # include <string.h>
  40 # include <openssl/err.h>
  41 # include "ec_lcl.h"
  42
  43 # if defined(__SIZEOF_INT128__) && __SIZEOF_INT128__==16
  44   /* even with gcc, the typedef won't work for 32-bit platforms */
  45 typedef __uint128_t uint128_t;  /* nonstandard; implemented by gcc on 64-bit
  46                                  * platforms */
  47 # else
  48 #  error "Your compiler doesn't appear to support 128-bit integer types"
  49 # endif
  50
  51 typedef uint8_t u8;
  52 typedef uint64_t u64;
  53
  54 /******************************************************************************/
  55 /*-
  56  * INTERNAL REPRESENTATION OF FIELD ELEMENTS
  57  *
  58  * Field elements are represented as a_0 + 2^56*a_1 + 2^112*a_2 + 2^168*a_3
  59  * using 64-bit coefficients called 'limbs',
  60  * and sometimes (for multiplication results) as
  61  * 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
  62  * using 128-bit coefficients called 'widelimbs'.
  63  * A 4-limb representation is an 'felem';
  64  * a 7-widelimb representation is a 'widefelem'.
  65  * Even within felems, bits of adjacent limbs overlap, and we don't always
  66  * reduce the representations: we ensure that inputs to each felem
  67  * multiplication satisfy a_i < 2^60, so outputs satisfy b_i < 4*2^60*2^60,
  68  * and fit into a 128-bit word without overflow. The coefficients are then
  69  * again partially reduced to obtain an felem satisfying a_i < 2^57.
  70  * We only reduce to the unique minimal representation at the end of the
  71  * computation.
  72  */
  73
  74 typedef uint64_t limb;
  75 typedef uint128_t widelimb;
  76
  77 typedef limb felem[4];
  78 typedef widelimb widefelem[7];
  79
  80 /*
  81  * Field element represented as a byte array. 28*8 = 224 bits is also the
  82  * group order size for the elliptic curve, and we also use this type for
  83  * scalars for point multiplication.
  84  */
  85 typedef u8 felem_bytearray[28];
  86
  87 static const felem_bytearray nistp224_curve_params[5] = {
  88     {0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, /* p */
  89      0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x00, 0x00, 0x00, 0x00,
  90      0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01},
  91     {0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, /* a */
  92      0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFE, 0xFF, 0xFF, 0xFF, 0xFF,
  93      0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFE},
  94     {0xB4, 0x05, 0x0A, 0x85, 0x0C, 0x04, 0xB3, 0xAB, 0xF5, 0x41, /* b */
  95      0x32, 0x56, 0x50, 0x44, 0xB0, 0xB7, 0xD7, 0xBF, 0xD8, 0xBA,
  96      0x27, 0x0B, 0x39, 0x43, 0x23, 0x55, 0xFF, 0xB4},
  97     {0xB7, 0x0E, 0x0C, 0xBD, 0x6B, 0xB4, 0xBF, 0x7F, 0x32, 0x13, /* x */
  98      0x90, 0xB9, 0x4A, 0x03, 0xC1, 0xD3, 0x56, 0xC2, 0x11, 0x22,
  99      0x34, 0x32, 0x80, 0xD6, 0x11, 0x5C, 0x1D, 0x21},
 100     {0xbd, 0x37, 0x63, 0x88, 0xb5, 0xf7, 0x23, 0xfb, 0x4c, 0x22, /* y */
 101      0xdf, 0xe6, 0xcd, 0x43, 0x75, 0xa0, 0x5a, 0x07, 0x47, 0x64,
 102      0x44, 0xd5, 0x81, 0x99, 0x85, 0x00, 0x7e, 0x34}
 103 };
 104
 105 /*-
 106  * Precomputed multiples of the standard generator
 107  * Points are given in coordinates (X, Y, Z) where Z normally is 1
 108  * (0 for the point at infinity).
 109  * For each field element, slice a_0 is word 0, etc.
 110  *
 111  * The table has 2 * 16 elements, starting with the following:
 112  * index | bits    | point
 113  * ------+---------+------------------------------
 114  *     0 | 0 0 0 0 | 0G
 115  *     1 | 0 0 0 1 | 1G
 116  *     2 | 0 0 1 0 | 2^56G
 117  *     3 | 0 0 1 1 | (2^56 + 1)G
 118  *     4 | 0 1 0 0 | 2^112G
 119  *     5 | 0 1 0 1 | (2^112 + 1)G
 120  *     6 | 0 1 1 0 | (2^112 + 2^56)G
 121  *     7 | 0 1 1 1 | (2^112 + 2^56 + 1)G
 122  *     8 | 1 0 0 0 | 2^168G
 123  *     9 | 1 0 0 1 | (2^168 + 1)G
 124  *    10 | 1 0 1 0 | (2^168 + 2^56)G
 125  *    11 | 1 0 1 1 | (2^168 + 2^56 + 1)G
 126  *    12 | 1 1 0 0 | (2^168 + 2^112)G
 127  *    13 | 1 1 0 1 | (2^168 + 2^112 + 1)G
 128  *    14 | 1 1 1 0 | (2^168 + 2^112 + 2^56)G
 129  *    15 | 1 1 1 1 | (2^168 + 2^112 + 2^56 + 1)G
 130  * followed by a copy of this with each element multiplied by 2^28.
 131  *
 132  * The reason for this is so that we can clock bits into four different
 133  * locations when doing simple scalar multiplies against the base point,
 134  * and then another four locations using the second 16 elements.
 135  */
 136 static const felem gmul[2][16][3] = {
 137 {{{0, 0, 0, 0},
 138   {0, 0, 0, 0},
 139   {0, 0, 0, 0}},
 140  {{0x3280d6115c1d21, 0xc1d356c2112234, 0x7f321390b94a03, 0xb70e0cbd6bb4bf},
 141   {0xd5819985007e34, 0x75a05a07476444, 0xfb4c22dfe6cd43, 0xbd376388b5f723},
 142   {1, 0, 0, 0}},
 143  {{0xfd9675666ebbe9, 0xbca7664d40ce5e, 0x2242df8d8a2a43, 0x1f49bbb0f99bc5},
 144   {0x29e0b892dc9c43, 0xece8608436e662, 0xdc858f185310d0, 0x9812dd4eb8d321},
 145   {1, 0, 0, 0}},
 146  {{0x6d3e678d5d8eb8, 0x559eed1cb362f1, 0x16e9a3bbce8a3f, 0xeedcccd8c2a748},
 147   {0xf19f90ed50266d, 0xabf2b4bf65f9df, 0x313865468fafec, 0x5cb379ba910a17},
 148   {1, 0, 0, 0}},
 149  {{0x0641966cab26e3, 0x91fb2991fab0a0, 0xefec27a4e13a0b, 0x0499aa8a5f8ebe},
 150   {0x7510407766af5d, 0x84d929610d5450, 0x81d77aae82f706, 0x6916f6d4338c5b},
 151   {1, 0, 0, 0}},
 152  {{0xea95ac3b1f15c6, 0x086000905e82d4, 0xdd323ae4d1c8b1, 0x932b56be7685a3},
 153   {0x9ef93dea25dbbf, 0x41665960f390f0, 0xfdec76dbe2a8a7, 0x523e80f019062a},
 154   {1, 0, 0, 0}},
 155  {{0x822fdd26732c73, 0xa01c83531b5d0f, 0x363f37347c1ba4, 0xc391b45c84725c},
 156   {0xbbd5e1b2d6ad24, 0xddfbcde19dfaec, 0xc393da7e222a7f, 0x1efb7890ede244},
 157   {1, 0, 0, 0}},
 158  {{0x4c9e90ca217da1, 0xd11beca79159bb, 0xff8d33c2c98b7c, 0x2610b39409f849},
 159   {0x44d1352ac64da0, 0xcdbb7b2c46b4fb, 0x966c079b753c89, 0xfe67e4e820b112},
 160   {1, 0, 0, 0}},
 161  {{0xe28cae2df5312d, 0xc71b61d16f5c6e, 0x79b7619a3e7c4c, 0x05c73240899b47},
 162   {0x9f7f6382c73e3a, 0x18615165c56bda, 0x641fab2116fd56, 0x72855882b08394},
 163   {1, 0, 0, 0}},
 164  {{0x0469182f161c09, 0x74a98ca8d00fb5, 0xb89da93489a3e0, 0x41c98768fb0c1d},
 165   {0xe5ea05fb32da81, 0x3dce9ffbca6855, 0x1cfe2d3fbf59e6, 0x0e5e03408738a7},
 166   {1, 0, 0, 0}},
 167  {{0xdab22b2333e87f, 0x4430137a5dd2f6, 0xe03ab9f738beb8, 0xcb0c5d0dc34f24},
 168   {0x764a7df0c8fda5, 0x185ba5c3fa2044, 0x9281d688bcbe50, 0xc40331df893881},
 169   {1, 0, 0, 0}},
 170  {{0xb89530796f0f60, 0xade92bd26909a3, 0x1a0c83fb4884da, 0x1765bf22a5a984},
 171   {0x772a9ee75db09e, 0x23bc6c67cec16f, 0x4c1edba8b14e2f, 0xe2a215d9611369},
 172   {1, 0, 0, 0}},
 173  {{0x571e509fb5efb3, 0xade88696410552, 0xc8ae85fada74fe, 0x6c7e4be83bbde3},
 174   {0xff9f51160f4652, 0xb47ce2495a6539, 0xa2946c53b582f4, 0x286d2db3ee9a60},
 175   {1, 0, 0, 0}},
 176  {{0x40bbd5081a44af, 0x0995183b13926c, 0xbcefba6f47f6d0, 0x215619e9cc0057},
 177   {0x8bc94d3b0df45e, 0xf11c54a3694f6f, 0x8631b93cdfe8b5, 0xe7e3f4b0982db9},
 178   {1, 0, 0, 0}},
 179  {{0xb17048ab3e1c7b, 0xac38f36ff8a1d8, 0x1c29819435d2c6, 0xc813132f4c07e9},
 180   {0x2891425503b11f, 0x08781030579fea, 0xf5426ba5cc9674, 0x1e28ebf18562bc},
 181   {1, 0, 0, 0}},
 182  {{0x9f31997cc864eb, 0x06cd91d28b5e4c, 0xff17036691a973, 0xf1aef351497c58},
 183   {0xdd1f2d600564ff, 0xdead073b1402db, 0x74a684435bd693, 0xeea7471f962558},
 184   {1, 0, 0, 0}}},
 185 {{{0, 0, 0, 0},
 186   {0, 0, 0, 0},
 187   {0, 0, 0, 0}},
 188  {{0x9665266dddf554, 0x9613d78b60ef2d, 0xce27a34cdba417, 0xd35ab74d6afc31},
 189   {0x85ccdd22deb15e, 0x2137e5783a6aab, 0xa141cffd8c93c6, 0x355a1830e90f2d},
 190   {1, 0, 0, 0}},
 191  {{0x1a494eadaade65, 0xd6da4da77fe53c, 0xe7992996abec86, 0x65c3553c6090e3},
 192   {0xfa610b1fb09346, 0xf1c6540b8a4aaf, 0xc51a13ccd3cbab, 0x02995b1b18c28a},
 193   {1, 0, 0, 0}},
 194  {{0x7874568e7295ef, 0x86b419fbe38d04, 0xdc0690a7550d9a, 0xd3966a44beac33},
 195   {0x2b7280ec29132f, 0xbeaa3b6a032df3, 0xdc7dd88ae41200, 0xd25e2513e3a100},
 196   {1, 0, 0, 0}},
 197  {{0x924857eb2efafd, 0xac2bce41223190, 0x8edaa1445553fc, 0x825800fd3562d5},
 198   {0x8d79148ea96621, 0x23a01c3dd9ed8d, 0xaf8b219f9416b5, 0xd8db0cc277daea},
 199   {1, 0, 0, 0}},
 200  {{0x76a9c3b1a700f0, 0xe9acd29bc7e691, 0x69212d1a6b0327, 0x6322e97fe154be},
 201   {0x469fc5465d62aa, 0x8d41ed18883b05, 0x1f8eae66c52b88, 0xe4fcbe9325be51},
 202   {1, 0, 0, 0}},
 203  {{0x825fdf583cac16, 0x020b857c7b023a, 0x683c17744b0165, 0x14ffd0a2daf2f1},
 204   {0x323b36184218f9, 0x4944ec4e3b47d4, 0xc15b3080841acf, 0x0bced4b01a28bb},
 205   {1, 0, 0, 0}},
 206  {{0x92ac22230df5c4, 0x52f33b4063eda8, 0xcb3f19870c0c93, 0x40064f2ba65233},
 207   {0xfe16f0924f8992, 0x012da25af5b517, 0x1a57bb24f723a6, 0x06f8bc76760def},
 208   {1, 0, 0, 0}},
 209  {{0x4a7084f7817cb9, 0xbcab0738ee9a78, 0x3ec11e11d9c326, 0xdc0fe90e0f1aae},
 210   {0xcf639ea5f98390, 0x5c350aa22ffb74, 0x9afae98a4047b7, 0x956ec2d617fc45},
 211   {1, 0, 0, 0}},
 212  {{0x4306d648c1be6a, 0x9247cd8bc9a462, 0xf5595e377d2f2e, 0xbd1c3caff1a52e},
 213   {0x045e14472409d0, 0x29f3e17078f773, 0x745a602b2d4f7d, 0x191837685cdfbb},
 214   {1, 0, 0, 0}},
 215  {{0x5b6ee254a8cb79, 0x4953433f5e7026, 0xe21faeb1d1def4, 0xc4c225785c09de},
 216   {0x307ce7bba1e518, 0x31b125b1036db8, 0x47e91868839e8f, 0xc765866e33b9f3},
 217   {1, 0, 0, 0}},
 218  {{0x3bfece24f96906, 0x4794da641e5093, 0xde5df64f95db26, 0x297ecd89714b05},
 219   {0x701bd3ebb2c3aa, 0x7073b4f53cb1d5, 0x13c5665658af16, 0x9895089d66fe58},
 220   {1, 0, 0, 0}},
 221  {{0x0fef05f78c4790, 0x2d773633b05d2e, 0x94229c3a951c94, 0xbbbd70df4911bb},
 222   {0xb2c6963d2c1168, 0x105f47a72b0d73, 0x9fdf6111614080, 0x7b7e94b39e67b0},
 223   {1, 0, 0, 0}},
 224  {{0xad1a7d6efbe2b3, 0xf012482c0da69d, 0x6b3bdf12438345, 0x40d7558d7aa4d9},
 225   {0x8a09fffb5c6d3d, 0x9a356e5d9ffd38, 0x5973f15f4f9b1c, 0xdcd5f59f63c3ea},
 226   {1, 0, 0, 0}},
 227  {{0xacf39f4c5ca7ab, 0x4c8071cc5fd737, 0xc64e3602cd1184, 0x0acd4644c9abba},
 228   {0x6c011a36d8bf6e, 0xfecd87ba24e32a, 0x19f6f56574fad8, 0x050b204ced9405},
 229   {1, 0, 0, 0}},
 230  {{0xed4f1cae7d9a96, 0x5ceef7ad94c40a, 0x778e4a3bf3ef9b, 0x7405783dc3b55e},
 231   {0x32477c61b6e8c6, 0xb46a97570f018b, 0x91176d0a7e95d1, 0x3df90fbc4c7d0e},
 232   {1, 0, 0, 0}}}
 233 };
 234
 235 /* Precomputation for the group generator. */
 236 struct nistp224_pre_comp_st {
 237     felem g_pre_comp[2][16][3];
 238     CRYPTO_REF_COUNT references;
 239     CRYPTO_RWLOCK *lock;
 240 };
 241
 242 const EC_METHOD *EC_GFp_nistp224_method(void)
 243 {
 244     static const EC_METHOD ret = {
 245         EC_FLAGS_DEFAULT_OCT,
 246         NID_X9_62_prime_field,
 247         ec_GFp_nistp224_group_init,
 248         ec_GFp_simple_group_finish,
 249         ec_GFp_simple_group_clear_finish,
 250         ec_GFp_nist_group_copy,
 251         ec_GFp_nistp224_group_set_curve,
 252         ec_GFp_simple_group_get_curve,
 253         ec_GFp_simple_group_get_degree,
 254         ec_group_simple_order_bits,
 255         ec_GFp_simple_group_check_discriminant,
 256         ec_GFp_simple_point_init,
 257         ec_GFp_simple_point_finish,
 258         ec_GFp_simple_point_clear_finish,
 259         ec_GFp_simple_point_copy,
 260         ec_GFp_simple_point_set_to_infinity,
 261         ec_GFp_simple_set_Jprojective_coordinates_GFp,
 262         ec_GFp_simple_get_Jprojective_coordinates_GFp,
 263         ec_GFp_simple_point_set_affine_coordinates,
 264         ec_GFp_nistp224_point_get_affine_coordinates,
 265         0 /* point_set_compressed_coordinates */ ,
 266         0 /* point2oct */ ,
 267         0 /* oct2point */ ,
 268         ec_GFp_simple_add,
 269         ec_GFp_simple_dbl,
 270         ec_GFp_simple_invert,
 271         ec_GFp_simple_is_at_infinity,
 272         ec_GFp_simple_is_on_curve,
 273         ec_GFp_simple_cmp,
 274         ec_GFp_simple_make_affine,
 275         ec_GFp_simple_points_make_affine,
 276         ec_GFp_nistp224_points_mul,
 277         ec_GFp_nistp224_precompute_mult,
 278         ec_GFp_nistp224_have_precompute_mult,
 279         ec_GFp_nist_field_mul,
 280         ec_GFp_nist_field_sqr,
 281         0 /* field_div */ ,
 282         ec_GFp_simple_field_inv,
 283         0 /* field_encode */ ,
 284         0 /* field_decode */ ,
 285         0,                      /* field_set_to_one */
 286         ec_key_simple_priv2oct,
 287         ec_key_simple_oct2priv,
 288         0, /* set private */
 289         ec_key_simple_generate_key,
 290         ec_key_simple_check_key,
 291         ec_key_simple_generate_public_key,
 292         0, /* keycopy */
 293         0, /* keyfinish */
 294         ecdh_simple_compute_key,
 295         0, /* field_inverse_mod_ord */
 296         0, /* blind_coordinates */
 297         0, /* ladder_pre */
 298         0, /* ladder_step */
 299         0  /* ladder_post */
 300     };
 301
 302     return &ret;
 303 }
 304
 305 /*
 306  * Helper functions to convert field elements to/from internal representation
 307  */
 308 static void bin28_to_felem(felem out, const u8 in[28])
 309 {
 310     out[0] = *((const uint64_t *)(in)) & 0x00ffffffffffffff;
 311     out[1] = (*((const uint64_t *)(in + 7))) & 0x00ffffffffffffff;
 312     out[2] = (*((const uint64_t *)(in + 14))) & 0x00ffffffffffffff;
 313     out[3] = (*((const uint64_t *)(in+20))) >> 8;
 314 }
 315
 316 static void felem_to_bin28(u8 out[28], const felem in)
 317 {
 318     unsigned i;
 319     for (i = 0; i < 7; ++i) {
 320         out[i] = in[0] >> (8 * i);
 321         out[i + 7] = in[1] >> (8 * i);
 322         out[i + 14] = in[2] >> (8 * i);
 323         out[i + 21] = in[3] >> (8 * i);
 324     }
 325 }
 326
 327 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */
 328 static void flip_endian(u8 *out, const u8 *in, unsigned len)
 329 {
 330     unsigned i;
 331     for (i = 0; i < len; ++i)
 332         out[i] = in[len - 1 - i];
 333 }
 334
 335 /* From OpenSSL BIGNUM to internal representation */
 336 static int BN_to_felem(felem out, const BIGNUM *bn)
 337 {
 338     felem_bytearray b_in;
 339     felem_bytearray b_out;
 340     unsigned num_bytes;
 341
 342     /* BN_bn2bin eats leading zeroes */
 343     memset(b_out, 0, sizeof(b_out));
 344     num_bytes = BN_num_bytes(bn);
 345     if (num_bytes > sizeof(b_out)) {
 346         ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
 347         return 0;
 348     }
 349     if (BN_is_negative(bn)) {
 350         ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
 351         return 0;
 352     }
 353     num_bytes = BN_bn2bin(bn, b_in);
 354     flip_endian(b_out, b_in, num_bytes);
 355     bin28_to_felem(out, b_out);
 356     return 1;
 357 }
 358
 359 /* From internal representation to OpenSSL BIGNUM */
 360 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
 361 {
 362     felem_bytearray b_in, b_out;
 363     felem_to_bin28(b_in, in);
 364     flip_endian(b_out, b_in, sizeof(b_out));
 365     return BN_bin2bn(b_out, sizeof(b_out), out);
 366 }
 367
 368 /******************************************************************************/
 369 /*-
 370  *                              FIELD OPERATIONS
 371  *
 372  * Field operations, using the internal representation of field elements.
 373  * NB! These operations are specific to our point multiplication and cannot be
 374  * expected to be correct in general - e.g., multiplication with a large scalar
 375  * will cause an overflow.
 376  *
 377  */
 378
 379 static void felem_one(felem out)
 380 {
 381     out[0] = 1;
 382     out[1] = 0;
 383     out[2] = 0;
 384     out[3] = 0;
 385 }
 386
 387 static void felem_assign(felem out, const felem in)
 388 {
 389     out[0] = in[0];
 390     out[1] = in[1];
 391     out[2] = in[2];
 392     out[3] = in[3];
 393 }
 394
 395 /* Sum two field elements: out += in */
 396 static void felem_sum(felem out, const felem in)
 397 {
 398     out[0] += in[0];
 399     out[1] += in[1];
 400     out[2] += in[2];
 401     out[3] += in[3];
 402 }
 403
 404 /* Subtract field elements: out -= in */
 405 /* Assumes in[i] < 2^57 */
 406 static void felem_diff(felem out, const felem in)
 407 {
 408     static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
 409     static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
 410     static const limb two58m42m2 = (((limb) 1) << 58) -
 411         (((limb) 1) << 42) - (((limb) 1) << 2);
 412
 413     /* Add 0 mod 2^224-2^96+1 to ensure out > in */
 414     out[0] += two58p2;
 415     out[1] += two58m42m2;
 416     out[2] += two58m2;
 417     out[3] += two58m2;
 418
 419     out[0] -= in[0];
 420     out[1] -= in[1];
 421     out[2] -= in[2];
 422     out[3] -= in[3];
 423 }
 424
 425 /* Subtract in unreduced 128-bit mode: out -= in */
 426 /* Assumes in[i] < 2^119 */
 427 static void widefelem_diff(widefelem out, const widefelem in)
 428 {
 429     static const widelimb two120 = ((widelimb) 1) << 120;
 430     static const widelimb two120m64 = (((widelimb) 1) << 120) -
 431         (((widelimb) 1) << 64);
 432     static const widelimb two120m104m64 = (((widelimb) 1) << 120) -
 433         (((widelimb) 1) << 104) - (((widelimb) 1) << 64);
 434
 435     /* Add 0 mod 2^224-2^96+1 to ensure out > in */
 436     out[0] += two120;
 437     out[1] += two120m64;
 438     out[2] += two120m64;
 439     out[3] += two120;
 440     out[4] += two120m104m64;
 441     out[5] += two120m64;
 442     out[6] += two120m64;
 443
 444     out[0] -= in[0];
 445     out[1] -= in[1];
 446     out[2] -= in[2];
 447     out[3] -= in[3];
 448     out[4] -= in[4];
 449     out[5] -= in[5];
 450     out[6] -= in[6];
 451 }
 452
 453 /* Subtract in mixed mode: out128 -= in64 */
 454 /* in[i] < 2^63 */
 455 static void felem_diff_128_64(widefelem out, const felem in)
 456 {
 457     static const widelimb two64p8 = (((widelimb) 1) << 64) +
 458         (((widelimb) 1) << 8);
 459     static const widelimb two64m8 = (((widelimb) 1) << 64) -
 460         (((widelimb) 1) << 8);
 461     static const widelimb two64m48m8 = (((widelimb) 1) << 64) -
 462         (((widelimb) 1) << 48) - (((widelimb) 1) << 8);
 463
 464     /* Add 0 mod 2^224-2^96+1 to ensure out > in */
 465     out[0] += two64p8;
 466     out[1] += two64m48m8;
 467     out[2] += two64m8;
 468     out[3] += two64m8;
 469
 470     out[0] -= in[0];
 471     out[1] -= in[1];
 472     out[2] -= in[2];
 473     out[3] -= in[3];
 474 }
 475
 476 /*
 477  * Multiply a field element by a scalar: out = out * scalar The scalars we
 478  * actually use are small, so results fit without overflow
 479  */
 480 static void felem_scalar(felem out, const limb scalar)
 481 {
 482     out[0] *= scalar;
 483     out[1] *= scalar;
 484     out[2] *= scalar;
 485     out[3] *= scalar;
 486 }
 487
 488 /*
 489  * Multiply an unreduced field element by a scalar: out = out * scalar The
 490  * scalars we actually use are small, so results fit without overflow
 491  */
 492 static void widefelem_scalar(widefelem out, const widelimb scalar)
 493 {
 494     out[0] *= scalar;
 495     out[1] *= scalar;
 496     out[2] *= scalar;
 497     out[3] *= scalar;
 498     out[4] *= scalar;
 499     out[5] *= scalar;
 500     out[6] *= scalar;
 501 }
 502
 503 /* Square a field element: out = in^2 */
 504 static void felem_square(widefelem out, const felem in)
 505 {
 506     limb tmp0, tmp1, tmp2;
 507     tmp0 = 2 * in[0];
 508     tmp1 = 2 * in[1];
 509     tmp2 = 2 * in[2];
 510     out[0] = ((widelimb) in[0]) * in[0];
 511     out[1] = ((widelimb) in[0]) * tmp1;
 512     out[2] = ((widelimb) in[0]) * tmp2 + ((widelimb) in[1]) * in[1];
 513     out[3] = ((widelimb) in[3]) * tmp0 + ((widelimb) in[1]) * tmp2;
 514     out[4] = ((widelimb) in[3]) * tmp1 + ((widelimb) in[2]) * in[2];
 515     out[5] = ((widelimb) in[3]) * tmp2;
 516     out[6] = ((widelimb) in[3]) * in[3];
 517 }
 518
 519 /* Multiply two field elements: out = in1 * in2 */
 520 static void felem_mul(widefelem out, const felem in1, const felem in2)
 521 {
 522     out[0] = ((widelimb) in1[0]) * in2[0];
 523     out[1] = ((widelimb) in1[0]) * in2[1] + ((widelimb) in1[1]) * in2[0];
 524     out[2] = ((widelimb) in1[0]) * in2[2] + ((widelimb) in1[1]) * in2[1] +
 525              ((widelimb) in1[2]) * in2[0];
 526     out[3] = ((widelimb) in1[0]) * in2[3] + ((widelimb) in1[1]) * in2[2] +
 527              ((widelimb) in1[2]) * in2[1] + ((widelimb) in1[3]) * in2[0];
 528     out[4] = ((widelimb) in1[1]) * in2[3] + ((widelimb) in1[2]) * in2[2] +
 529              ((widelimb) in1[3]) * in2[1];
 530     out[5] = ((widelimb) in1[2]) * in2[3] + ((widelimb) in1[3]) * in2[2];
 531     out[6] = ((widelimb) in1[3]) * in2[3];
 532 }
 533
 534 /*-
 535  * Reduce seven 128-bit coefficients to four 64-bit coefficients.
 536  * Requires in[i] < 2^126,
 537  * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] <= 2^56 + 2^16 */
 538 static void felem_reduce(felem out, const widefelem in)
 539 {
 540     static const widelimb two127p15 = (((widelimb) 1) << 127) +
 541         (((widelimb) 1) << 15);
 542     static const widelimb two127m71 = (((widelimb) 1) << 127) -
 543         (((widelimb) 1) << 71);
 544     static const widelimb two127m71m55 = (((widelimb) 1) << 127) -
 545         (((widelimb) 1) << 71) - (((widelimb) 1) << 55);
 546     widelimb output[5];
 547
 548     /* Add 0 mod 2^224-2^96+1 to ensure all differences are positive */
 549     output[0] = in[0] + two127p15;
 550     output[1] = in[1] + two127m71m55;
 551     output[2] = in[2] + two127m71;
 552     output[3] = in[3];
 553     output[4] = in[4];
 554
 555     /* Eliminate in[4], in[5], in[6] */
 556     output[4] += in[6] >> 16;
 557     output[3] += (in[6] & 0xffff) << 40;
 558     output[2] -= in[6];
 559
 560     output[3] += in[5] >> 16;
 561     output[2] += (in[5] & 0xffff) << 40;
 562     output[1] -= in[5];
 563
 564     output[2] += output[4] >> 16;
 565     output[1] += (output[4] & 0xffff) << 40;
 566     output[0] -= output[4];
 567
 568     /* Carry 2 -> 3 -> 4 */
 569     output[3] += output[2] >> 56;
 570     output[2] &= 0x00ffffffffffffff;
 571
 572     output[4] = output[3] >> 56;
 573     output[3] &= 0x00ffffffffffffff;
 574
 575     /* Now output[2] < 2^56, output[3] < 2^56, output[4] < 2^72 */
 576
 577     /* Eliminate output[4] */
 578     output[2] += output[4] >> 16;
 579     /* output[2] < 2^56 + 2^56 = 2^57 */
 580     output[1] += (output[4] & 0xffff) << 40;
 581     output[0] -= output[4];
 582
 583     /* Carry 0 -> 1 -> 2 -> 3 */
 584     output[1] += output[0] >> 56;
 585     out[0] = output[0] & 0x00ffffffffffffff;
 586
 587     output[2] += output[1] >> 56;
 588     /* output[2] < 2^57 + 2^72 */
 589     out[1] = output[1] & 0x00ffffffffffffff;
 590     output[3] += output[2] >> 56;
 591     /* output[3] <= 2^56 + 2^16 */
 592     out[2] = output[2] & 0x00ffffffffffffff;
 593
 594     /*-
 595      * out[0] < 2^56, out[1] < 2^56, out[2] < 2^56,
 596      * out[3] <= 2^56 + 2^16 (due to final carry),
 597      * so out < 2*p
 598      */
 599     out[3] = output[3];
 600 }
 601
 602 static void felem_square_reduce(felem out, const felem in)
 603 {
 604     widefelem tmp;
 605     felem_square(tmp, in);
 606     felem_reduce(out, tmp);
 607 }
 608
 609 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
 610 {
 611     widefelem tmp;
 612     felem_mul(tmp, in1, in2);
 613     felem_reduce(out, tmp);
 614 }
 615
 616 /*
 617  * Reduce to unique minimal representation. Requires 0 <= in < 2*p (always
 618  * call felem_reduce first)
 619  */
 620 static void felem_contract(felem out, const felem in)
 621 {
 622     static const int64_t two56 = ((limb) 1) << 56;
 623     /* 0 <= in < 2*p, p = 2^224 - 2^96 + 1 */
 624     /* if in > p , reduce in = in - 2^224 + 2^96 - 1 */
 625     int64_t tmp[4], a;
 626     tmp[0] = in[0];
 627     tmp[1] = in[1];
 628     tmp[2] = in[2];
 629     tmp[3] = in[3];
 630     /* Case 1: a = 1 iff in >= 2^224 */
 631     a = (in[3] >> 56);
 632     tmp[0] -= a;
 633     tmp[1] += a << 40;
 634     tmp[3] &= 0x00ffffffffffffff;
 635     /*
 636      * Case 2: a = 0 iff p <= in < 2^224, i.e., the high 128 bits are all 1
 637      * and the lower part is non-zero
 638      */
 639     a = ((in[3] & in[2] & (in[1] | 0x000000ffffffffff)) + 1) |
 640         (((int64_t) (in[0] + (in[1] & 0x000000ffffffffff)) - 1) >> 63);
 641     a &= 0x00ffffffffffffff;
 642     /* turn a into an all-one mask (if a = 0) or an all-zero mask */
 643     a = (a - 1) >> 63;
 644     /* subtract 2^224 - 2^96 + 1 if a is all-one */
 645     tmp[3] &= a ^ 0xffffffffffffffff;
 646     tmp[2] &= a ^ 0xffffffffffffffff;
 647     tmp[1] &= (a ^ 0xffffffffffffffff) | 0x000000ffffffffff;
 648     tmp[0] -= 1 & a;
 649
 650     /*
 651      * eliminate negative coefficients: if tmp[0] is negative, tmp[1] must be
 652      * non-zero, so we only need one step
 653      */
 654     a = tmp[0] >> 63;
 655     tmp[0] += two56 & a;
 656     tmp[1] -= 1 & a;
 657
 658     /* carry 1 -> 2 -> 3 */
 659     tmp[2] += tmp[1] >> 56;
 660     tmp[1] &= 0x00ffffffffffffff;
 661
 662     tmp[3] += tmp[2] >> 56;
 663     tmp[2] &= 0x00ffffffffffffff;
 664
 665     /* Now 0 <= out < p */
 666     out[0] = tmp[0];
 667     out[1] = tmp[1];
 668     out[2] = tmp[2];
 669     out[3] = tmp[3];
 670 }
 671
 672 /*
 673  * Get negative value: out = -in
 674  * Requires in[i] < 2^63,
 675  * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] <= 2^56 + 2^16
 676  */
 677 static void felem_neg(felem out, const felem in)
 678 {
 679     widefelem tmp = {0};
 680     felem_diff_128_64(tmp, in);
 681     felem_reduce(out, tmp);
 682 }
 683
 684 /*
 685  * Zero-check: returns 1 if input is 0, and 0 otherwise. We know that field
 686  * elements are reduced to in < 2^225, so we only need to check three cases:
 687  * 0, 2^224 - 2^96 + 1, and 2^225 - 2^97 + 2
 688  */
 689 static limb felem_is_zero(const felem in)
 690 {
 691     limb zero, two224m96p1, two225m97p2;
 692
 693     zero = in[0] | in[1] | in[2] | in[3];
 694     zero = (((int64_t) (zero) - 1) >> 63) & 1;
 695     two224m96p1 = (in[0] ^ 1) | (in[1] ^ 0x00ffff0000000000)
 696         | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x00ffffffffffffff);
 697     two224m96p1 = (((int64_t) (two224m96p1) - 1) >> 63) & 1;
 698     two225m97p2 = (in[0] ^ 2) | (in[1] ^ 0x00fffe0000000000)
 699         | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x01ffffffffffffff);
 700     two225m97p2 = (((int64_t) (two225m97p2) - 1) >> 63) & 1;
 701     return (zero | two224m96p1 | two225m97p2);
 702 }
 703
 704 static int felem_is_zero_int(const void *in)
 705 {
 706     return (int)(felem_is_zero(in) & ((limb) 1));
 707 }
 708
 709 /* Invert a field element */
 710 /* Computation chain copied from djb's code */
 711 static void felem_inv(felem out, const felem in)
 712 {
 713     felem ftmp, ftmp2, ftmp3, ftmp4;
 714     widefelem tmp;
 715     unsigned i;
 716
 717     felem_square(tmp, in);
 718     felem_reduce(ftmp, tmp);    /* 2 */
 719     felem_mul(tmp, in, ftmp);
 720     felem_reduce(ftmp, tmp);    /* 2^2 - 1 */
 721     felem_square(tmp, ftmp);
 722     felem_reduce(ftmp, tmp);    /* 2^3 - 2 */
 723     felem_mul(tmp, in, ftmp);
 724     felem_reduce(ftmp, tmp);    /* 2^3 - 1 */
 725     felem_square(tmp, ftmp);
 726     felem_reduce(ftmp2, tmp);   /* 2^4 - 2 */
 727     felem_square(tmp, ftmp2);
 728     felem_reduce(ftmp2, tmp);   /* 2^5 - 4 */
 729     felem_square(tmp, ftmp2);
 730     felem_reduce(ftmp2, tmp);   /* 2^6 - 8 */
 731     felem_mul(tmp, ftmp2, ftmp);
 732     felem_reduce(ftmp, tmp);    /* 2^6 - 1 */
 733     felem_square(tmp, ftmp);
 734     felem_reduce(ftmp2, tmp);   /* 2^7 - 2 */
 735     for (i = 0; i < 5; ++i) {   /* 2^12 - 2^6 */
 736         felem_square(tmp, ftmp2);
 737         felem_reduce(ftmp2, tmp);
 738     }
 739     felem_mul(tmp, ftmp2, ftmp);
 740     felem_reduce(ftmp2, tmp);   /* 2^12 - 1 */
 741     felem_square(tmp, ftmp2);
 742     felem_reduce(ftmp3, tmp);   /* 2^13 - 2 */
 743     for (i = 0; i < 11; ++i) {  /* 2^24 - 2^12 */
 744         felem_square(tmp, ftmp3);
 745         felem_reduce(ftmp3, tmp);
 746     }
 747     felem_mul(tmp, ftmp3, ftmp2);
 748     felem_reduce(ftmp2, tmp);   /* 2^24 - 1 */
 749     felem_square(tmp, ftmp2);
 750     felem_reduce(ftmp3, tmp);   /* 2^25 - 2 */
 751     for (i = 0; i < 23; ++i) {  /* 2^48 - 2^24 */
 752         felem_square(tmp, ftmp3);
 753         felem_reduce(ftmp3, tmp);
 754     }
 755     felem_mul(tmp, ftmp3, ftmp2);
 756     felem_reduce(ftmp3, tmp);   /* 2^48 - 1 */
 757     felem_square(tmp, ftmp3);
 758     felem_reduce(ftmp4, tmp);   /* 2^49 - 2 */
 759     for (i = 0; i < 47; ++i) {  /* 2^96 - 2^48 */
 760         felem_square(tmp, ftmp4);
 761         felem_reduce(ftmp4, tmp);
 762     }
 763     felem_mul(tmp, ftmp3, ftmp4);
 764     felem_reduce(ftmp3, tmp);   /* 2^96 - 1 */
 765     felem_square(tmp, ftmp3);
 766     felem_reduce(ftmp4, tmp);   /* 2^97 - 2 */
 767     for (i = 0; i < 23; ++i) {  /* 2^120 - 2^24 */
 768         felem_square(tmp, ftmp4);
 769         felem_reduce(ftmp4, tmp);
 770     }
 771     felem_mul(tmp, ftmp2, ftmp4);
 772     felem_reduce(ftmp2, tmp);   /* 2^120 - 1 */
 773     for (i = 0; i < 6; ++i) {   /* 2^126 - 2^6 */
 774         felem_square(tmp, ftmp2);
 775         felem_reduce(ftmp2, tmp);
 776     }
 777     felem_mul(tmp, ftmp2, ftmp);
 778     felem_reduce(ftmp, tmp);    /* 2^126 - 1 */
 779     felem_square(tmp, ftmp);
 780     felem_reduce(ftmp, tmp);    /* 2^127 - 2 */
 781     felem_mul(tmp, ftmp, in);
 782     felem_reduce(ftmp, tmp);    /* 2^127 - 1 */
 783     for (i = 0; i < 97; ++i) {  /* 2^224 - 2^97 */
 784         felem_square(tmp, ftmp);
 785         felem_reduce(ftmp, tmp);
 786     }
 787     felem_mul(tmp, ftmp, ftmp3);
 788     felem_reduce(out, tmp);     /* 2^224 - 2^96 - 1 */
 789 }
 790
 791 /*
 792  * Copy in constant time: if icopy == 1, copy in to out, if icopy == 0, copy
 793  * out to itself.
 794  */
 795 static void copy_conditional(felem out, const felem in, limb icopy)
 796 {
 797     unsigned i;
 798     /*
 799      * icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one
 800      */
 801     const limb copy = -icopy;
 802     for (i = 0; i < 4; ++i) {
 803         const limb tmp = copy & (in[i] ^ out[i]);
 804         out[i] ^= tmp;
 805     }
 806 }
 807
 808 /******************************************************************************/
 809 /*-
 810  *                       ELLIPTIC CURVE POINT OPERATIONS
 811  *
 812  * Points are represented in Jacobian projective coordinates:
 813  * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
 814  * or to the point at infinity if Z == 0.
 815  *
 816  */
 817
 818 /*-
 819  * Double an elliptic curve point:
 820  * (X', Y', Z') = 2 * (X, Y, Z), where
 821  * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
 822  * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^4
 823  * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
 824  * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
 825  * while x_out == y_in is not (maybe this works, but it's not tested).
 826  */
 827 static void
 828 point_double(felem x_out, felem y_out, felem z_out,
 829              const felem x_in, const felem y_in, const felem z_in)
 830 {
 831     widefelem tmp, tmp2;
 832     felem delta, gamma, beta, alpha, ftmp, ftmp2;
 833
 834     felem_assign(ftmp, x_in);
 835     felem_assign(ftmp2, x_in);
 836
 837     /* delta = z^2 */
 838     felem_square(tmp, z_in);
 839     felem_reduce(delta, tmp);
 840
 841     /* gamma = y^2 */
 842     felem_square(tmp, y_in);
 843     felem_reduce(gamma, tmp);
 844
 845     /* beta = x*gamma */
 846     felem_mul(tmp, x_in, gamma);
 847     felem_reduce(beta, tmp);
 848
 849     /* alpha = 3*(x-delta)*(x+delta) */
 850     felem_diff(ftmp, delta);
 851     /* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
 852     felem_sum(ftmp2, delta);
 853     /* ftmp2[i] < 2^57 + 2^57 = 2^58 */
 854     felem_scalar(ftmp2, 3);
 855     /* ftmp2[i] < 3 * 2^58 < 2^60 */
 856     felem_mul(tmp, ftmp, ftmp2);
 857     /* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
 858     felem_reduce(alpha, tmp);
 859
 860     /* x' = alpha^2 - 8*beta */
 861     felem_square(tmp, alpha);
 862     /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 863     felem_assign(ftmp, beta);
 864     felem_scalar(ftmp, 8);
 865     /* ftmp[i] < 8 * 2^57 = 2^60 */
 866     felem_diff_128_64(tmp, ftmp);
 867     /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
 868     felem_reduce(x_out, tmp);
 869
 870     /* z' = (y + z)^2 - gamma - delta */
 871     felem_sum(delta, gamma);
 872     /* delta[i] < 2^57 + 2^57 = 2^58 */
 873     felem_assign(ftmp, y_in);
 874     felem_sum(ftmp, z_in);
 875     /* ftmp[i] < 2^57 + 2^57 = 2^58 */
 876     felem_square(tmp, ftmp);
 877     /* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
 878     felem_diff_128_64(tmp, delta);
 879     /* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
 880     felem_reduce(z_out, tmp);
 881
 882     /* y' = alpha*(4*beta - x') - 8*gamma^2 */
 883     felem_scalar(beta, 4);
 884     /* beta[i] < 4 * 2^57 = 2^59 */
 885     felem_diff(beta, x_out);
 886     /* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
 887     felem_mul(tmp, alpha, beta);
 888     /* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
 889     felem_square(tmp2, gamma);
 890     /* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
 891     widefelem_scalar(tmp2, 8);
 892     /* tmp2[i] < 8 * 2^116 = 2^119 */
 893     widefelem_diff(tmp, tmp2);
 894     /* tmp[i] < 2^119 + 2^120 < 2^121 */
 895     felem_reduce(y_out, tmp);
 896 }
 897
 898 /*-
 899  * Add two elliptic curve points:
 900  * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
 901  * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
 902  * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
 903  * 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) -
 904  *        Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
 905  * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2)
 906  *
 907  * This runs faster if 'mixed' is set, which requires Z_2 = 1 or Z_2 = 0.
 908  */
 909
 910 /*
 911  * This function is not entirely constant-time: it includes a branch for
 912  * checking whether the two input points are equal, (while not equal to the
 913  * point at infinity). This case never happens during single point
 914  * multiplication, so there is no timing leak for ECDH or ECDSA signing.
 915  */
 916 static void point_add(felem x3, felem y3, felem z3,
 917                       const felem x1, const felem y1, const felem z1,
 918                       const int mixed, const felem x2, const felem y2,
 919                       const felem z2)
 920 {
 921     felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, x_out, y_out, z_out;
 922     widefelem tmp, tmp2;
 923     limb z1_is_zero, z2_is_zero, x_equal, y_equal;
 924
 925     if (!mixed) {
 926         /* ftmp2 = z2^2 */
 927         felem_square(tmp, z2);
 928         felem_reduce(ftmp2, tmp);
 929
 930         /* ftmp4 = z2^3 */
 931         felem_mul(tmp, ftmp2, z2);
 932         felem_reduce(ftmp4, tmp);
 933
 934         /* ftmp4 = z2^3*y1 */
 935         felem_mul(tmp2, ftmp4, y1);
 936         felem_reduce(ftmp4, tmp2);
 937
 938         /* ftmp2 = z2^2*x1 */
 939         felem_mul(tmp2, ftmp2, x1);
 940         felem_reduce(ftmp2, tmp2);
 941     } else {
 942         /*
 943          * We'll assume z2 = 1 (special case z2 = 0 is handled later)
 944          */
 945
 946         /* ftmp4 = z2^3*y1 */
 947         felem_assign(ftmp4, y1);
 948
 949         /* ftmp2 = z2^2*x1 */
 950         felem_assign(ftmp2, x1);
 951     }
 952
 953     /* ftmp = z1^2 */
 954     felem_square(tmp, z1);
 955     felem_reduce(ftmp, tmp);
 956
 957     /* ftmp3 = z1^3 */
 958     felem_mul(tmp, ftmp, z1);
 959     felem_reduce(ftmp3, tmp);
 960
 961     /* tmp = z1^3*y2 */
 962     felem_mul(tmp, ftmp3, y2);
 963     /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 964
 965     /* ftmp3 = z1^3*y2 - z2^3*y1 */
 966     felem_diff_128_64(tmp, ftmp4);
 967     /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
 968     felem_reduce(ftmp3, tmp);
 969
 970     /* tmp = z1^2*x2 */
 971     felem_mul(tmp, ftmp, x2);
 972     /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
 973
 974     /* ftmp = z1^2*x2 - z2^2*x1 */
 975     felem_diff_128_64(tmp, ftmp2);
 976     /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
 977     felem_reduce(ftmp, tmp);
 978
 979     /*
 980      * the formulae are incorrect if the points are equal so we check for
 981      * this and do doubling if this happens
 982      */
 983     x_equal = felem_is_zero(ftmp);
 984     y_equal = felem_is_zero(ftmp3);
 985     z1_is_zero = felem_is_zero(z1);
 986     z2_is_zero = felem_is_zero(z2);
 987     /* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
 988     if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
 989         point_double(x3, y3, z3, x1, y1, z1);
 990         return;
 991     }
 992
 993     /* ftmp5 = z1*z2 */
 994     if (!mixed) {
 995         felem_mul(tmp, z1, z2);
 996         felem_reduce(ftmp5, tmp);
 997     } else {
 998         /* special case z2 = 0 is handled later */
 999         felem_assign(ftmp5, z1);
1000     }
1001
1002     /* z_out = (z1^2*x2 - z2^2*x1)*(z1*z2) */
1003     felem_mul(tmp, ftmp, ftmp5);
1004     felem_reduce(z_out, tmp);
1005
1006     /* ftmp = (z1^2*x2 - z2^2*x1)^2 */
1007     felem_assign(ftmp5, ftmp);
1008     felem_square(tmp, ftmp);
1009     felem_reduce(ftmp, tmp);
1010
1011     /* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
1012     felem_mul(tmp, ftmp, ftmp5);
1013     felem_reduce(ftmp5, tmp);
1014
1015     /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
1016     felem_mul(tmp, ftmp2, ftmp);
1017     felem_reduce(ftmp2, tmp);
1018
1019     /* tmp = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
1020     felem_mul(tmp, ftmp4, ftmp5);
1021     /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
1022
1023     /* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
1024     felem_square(tmp2, ftmp3);
1025     /* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
1026
1027     /* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
1028     felem_diff_128_64(tmp2, ftmp5);
1029     /* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
1030
1031     /* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
1032     felem_assign(ftmp5, ftmp2);
1033     felem_scalar(ftmp5, 2);
1034     /* ftmp5[i] < 2 * 2^57 = 2^58 */
1035
1036     /*-
1037      * x_out = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
1038      *  2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2
1039      */
1040     felem_diff_128_64(tmp2, ftmp5);
1041     /* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
1042     felem_reduce(x_out, tmp2);
1043
1044     /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out */
1045     felem_diff(ftmp2, x_out);
1046     /* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
1047
1048     /*
1049      * tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out)
1050      */
1051     felem_mul(tmp2, ftmp3, ftmp2);
1052     /* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
1053
1054     /*-
1055      * y_out = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) -
1056      *  z2^3*y1*(z1^2*x2 - z2^2*x1)^3
1057      */
1058     widefelem_diff(tmp2, tmp);
1059     /* tmp2[i] < 2^118 + 2^120 < 2^121 */
1060     felem_reduce(y_out, tmp2);
1061
1062     /*
1063      * the result (x_out, y_out, z_out) is incorrect if one of the inputs is
1064      * the point at infinity, so we need to check for this separately
1065      */
1066
1067     /*
1068      * if point 1 is at infinity, copy point 2 to output, and vice versa
1069      */
1070     copy_conditional(x_out, x2, z1_is_zero);
1071     copy_conditional(x_out, x1, z2_is_zero);
1072     copy_conditional(y_out, y2, z1_is_zero);
1073     copy_conditional(y_out, y1, z2_is_zero);
1074     copy_conditional(z_out, z2, z1_is_zero);
1075     copy_conditional(z_out, z1, z2_is_zero);
1076     felem_assign(x3, x_out);
1077     felem_assign(y3, y_out);
1078     felem_assign(z3, z_out);
1079 }
1080
1081 /*
1082  * select_point selects the |idx|th point from a precomputation table and
1083  * copies it to out.
1084  * The pre_comp array argument should be size of |size| argument
1085  */
1086 static void select_point(const u64 idx, unsigned int size,
1087                          const felem pre_comp[][3], felem out[3])
1088 {
1089     unsigned i, j;
1090     limb *outlimbs = &out[0][0];
1091
1092     memset(out, 0, sizeof(*out) * 3);
1093     for (i = 0; i < size; i++) {
1094         const limb *inlimbs = &pre_comp[i][0][0];
1095         u64 mask = i ^ idx;
1096         mask |= mask >> 4;
1097         mask |= mask >> 2;
1098         mask |= mask >> 1;
1099         mask &= 1;
1100         mask--;
1101         for (j = 0; j < 4 * 3; j++)
1102             outlimbs[j] |= inlimbs[j] & mask;
1103     }
1104 }
1105
1106 /* get_bit returns the |i|th bit in |in| */
1107 static char get_bit(const felem_bytearray in, unsigned i)
1108 {
1109     if (i >= 224)
1110         return 0;
1111     return (in[i >> 3] >> (i & 7)) & 1;
1112 }
1113
1114 /*
1115  * Interleaved point multiplication using precomputed point multiples: The
1116  * small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[], the scalars
1117  * in scalars[]. If g_scalar is non-NULL, we also add this multiple of the
1118  * generator, using certain (large) precomputed multiples in g_pre_comp.
1119  * Output point (X, Y, Z) is stored in x_out, y_out, z_out
1120  */
1121 static void batch_mul(felem x_out, felem y_out, felem z_out,
1122                       const felem_bytearray scalars[],
1123                       const unsigned num_points, const u8 *g_scalar,
1124                       const int mixed, const felem pre_comp[][17][3],
1125                       const felem g_pre_comp[2][16][3])
1126 {
1127     int i, skip;
1128     unsigned num;
1129     unsigned gen_mul = (g_scalar != NULL);
1130     felem nq[3], tmp[4];
1131     u64 bits;
1132     u8 sign, digit;
1133
1134     /* set nq to the point at infinity */
1135     memset(nq, 0, sizeof(nq));
1136
1137     /*
1138      * Loop over all scalars msb-to-lsb, interleaving additions of multiples
1139      * of the generator (two in each of the last 28 rounds) and additions of
1140      * other points multiples (every 5th round).
1141      */
1142     skip = 1;                   /* save two point operations in the first
1143                                  * round */
1144     for (i = (num_points ? 220 : 27); i >= 0; --i) {
1145         /* double */
1146         if (!skip)
1147             point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1148
1149         /* add multiples of the generator */
1150         if (gen_mul && (i <= 27)) {
1151             /* first, look 28 bits upwards */
1152             bits = get_bit(g_scalar, i + 196) << 3;
1153             bits |= get_bit(g_scalar, i + 140) << 2;
1154             bits |= get_bit(g_scalar, i + 84) << 1;
1155             bits |= get_bit(g_scalar, i + 28);
1156             /* select the point to add, in constant time */
1157             select_point(bits, 16, g_pre_comp[1], tmp);
1158
1159             if (!skip) {
1160                 /* value 1 below is argument for "mixed" */
1161                 point_add(nq[0], nq[1], nq[2],
1162                           nq[0], nq[1], nq[2], 1, tmp[0], tmp[1], tmp[2]);
1163             } else {
1164                 memcpy(nq, tmp, 3 * sizeof(felem));
1165                 skip = 0;
1166             }
1167
1168             /* second, look at the current position */
1169             bits = get_bit(g_scalar, i + 168) << 3;
1170             bits |= get_bit(g_scalar, i + 112) << 2;
1171             bits |= get_bit(g_scalar, i + 56) << 1;
1172             bits |= get_bit(g_scalar, i);
1173             /* select the point to add, in constant time */
1174             select_point(bits, 16, g_pre_comp[0], tmp);
1175             point_add(nq[0], nq[1], nq[2],
1176                       nq[0], nq[1], nq[2],
1177                       1 /* mixed */ , tmp[0], tmp[1], tmp[2]);
1178         }
1179
1180         /* do other additions every 5 doublings */
1181         if (num_points && (i % 5 == 0)) {
1182             /* loop over all scalars */
1183             for (num = 0; num < num_points; ++num) {
1184                 bits = get_bit(scalars[num], i + 4) << 5;
1185                 bits |= get_bit(scalars[num], i + 3) << 4;
1186                 bits |= get_bit(scalars[num], i + 2) << 3;
1187                 bits |= get_bit(scalars[num], i + 1) << 2;
1188                 bits |= get_bit(scalars[num], i) << 1;
1189                 bits |= get_bit(scalars[num], i - 1);
1190                 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1191
1192                 /* select the point to add or subtract */
1193                 select_point(digit, 17, pre_comp[num], tmp);
1194                 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative
1195                                             * point */
1196                 copy_conditional(tmp[1], tmp[3], sign);
1197
1198                 if (!skip) {
1199                     point_add(nq[0], nq[1], nq[2],
1200                               nq[0], nq[1], nq[2],
1201                               mixed, tmp[0], tmp[1], tmp[2]);
1202                 } else {
1203                     memcpy(nq, tmp, 3 * sizeof(felem));
1204                     skip = 0;
1205                 }
1206             }
1207         }
1208     }
1209     felem_assign(x_out, nq[0]);
1210     felem_assign(y_out, nq[1]);
1211     felem_assign(z_out, nq[2]);
1212 }
1213
1214 /******************************************************************************/
1215 /*
1216  * FUNCTIONS TO MANAGE PRECOMPUTATION
1217  */
1218
1219 static NISTP224_PRE_COMP *nistp224_pre_comp_new(void)
1220 {
1221     NISTP224_PRE_COMP *ret = OPENSSL_zalloc(sizeof(*ret));
1222
1223     if (!ret) {
1224         ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1225         return ret;
1226     }
1227
1228     ret->references = 1;
1229
1230     ret->lock = CRYPTO_THREAD_lock_new();
1231     if (ret->lock == NULL) {
1232         ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1233         OPENSSL_free(ret);
1234         return NULL;
1235     }
1236     return ret;
1237 }
1238
1239 NISTP224_PRE_COMP *EC_nistp224_pre_comp_dup(NISTP224_PRE_COMP *p)
1240 {
1241     int i;
1242     if (p != NULL)
1243         CRYPTO_UP_REF(&p->references, &i, p->lock);
1244     return p;
1245 }
1246
1247 void EC_nistp224_pre_comp_free(NISTP224_PRE_COMP *p)
1248 {
1249     int i;
1250
1251     if (p == NULL)
1252         return;
1253
1254     CRYPTO_DOWN_REF(&p->references, &i, p->lock);
1255     REF_PRINT_COUNT("EC_nistp224", x);
1256     if (i > 0)
1257         return;
1258     REF_ASSERT_ISNT(i < 0);
1259
1260     CRYPTO_THREAD_lock_free(p->lock);
1261     OPENSSL_free(p);
1262 }
1263
1264 /******************************************************************************/
1265 /*
1266  * OPENSSL EC_METHOD FUNCTIONS
1267  */
1268
1269 int ec_GFp_nistp224_group_init(EC_GROUP *group)
1270 {
1271     int ret;
1272     ret = ec_GFp_simple_group_init(group);
1273     group->a_is_minus3 = 1;
1274     return ret;
1275 }
1276
1277 int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1278                                     const BIGNUM *a, const BIGNUM *b,
1279                                     BN_CTX *ctx)
1280 {
1281     int ret = 0;
1282     BN_CTX *new_ctx = NULL;
1283     BIGNUM *curve_p, *curve_a, *curve_b;
1284
1285     if (ctx == NULL)
1286         if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1287             return 0;
1288     BN_CTX_start(ctx);
1289     curve_p = BN_CTX_get(ctx);
1290     curve_a = BN_CTX_get(ctx);
1291     curve_b = BN_CTX_get(ctx);
1292     if (curve_b == NULL)
1293         goto err;
1294     BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
1295     BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
1296     BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
1297     if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) || (BN_cmp(curve_b, b))) {
1298         ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
1299               EC_R_WRONG_CURVE_PARAMETERS);
1300         goto err;
1301     }
1302     group->field_mod_func = BN_nist_mod_224;
1303     ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1304  err:
1305     BN_CTX_end(ctx);
1306     BN_CTX_free(new_ctx);
1307     return ret;
1308 }
1309
1310 /*
1311  * Takes the Jacobian coordinates (X, Y, Z) of a point and returns (X', Y') =
1312  * (X/Z^2, Y/Z^3)
1313  */
1314 int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
1315                                                  const EC_POINT *point,
1316                                                  BIGNUM *x, BIGNUM *y,
1317                                                  BN_CTX *ctx)
1318 {
1319     felem z1, z2, x_in, y_in, x_out, y_out;
1320     widefelem tmp;
1321
1322     if (EC_POINT_is_at_infinity(group, point)) {
1323         ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1324               EC_R_POINT_AT_INFINITY);
1325         return 0;
1326     }
1327     if ((!BN_to_felem(x_in, point->X)) || (!BN_to_felem(y_in, point->Y)) ||
1328         (!BN_to_felem(z1, point->Z)))
1329         return 0;
1330     felem_inv(z2, z1);
1331     felem_square(tmp, z2);
1332     felem_reduce(z1, tmp);
1333     felem_mul(tmp, x_in, z1);
1334     felem_reduce(x_in, tmp);
1335     felem_contract(x_out, x_in);
1336     if (x != NULL) {
1337         if (!felem_to_BN(x, x_out)) {
1338             ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1339                   ERR_R_BN_LIB);
1340             return 0;
1341         }
1342     }
1343     felem_mul(tmp, z1, z2);
1344     felem_reduce(z1, tmp);
1345     felem_mul(tmp, y_in, z1);
1346     felem_reduce(y_in, tmp);
1347     felem_contract(y_out, y_in);
1348     if (y != NULL) {
1349         if (!felem_to_BN(y, y_out)) {
1350             ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1351                   ERR_R_BN_LIB);
1352             return 0;
1353         }
1354     }
1355     return 1;
1356 }
1357
1358 static void make_points_affine(size_t num, felem points[ /* num */ ][3],
1359                                felem tmp_felems[ /* num+1 */ ])
1360 {
1361     /*
1362      * Runs in constant time, unless an input is the point at infinity (which
1363      * normally shouldn't happen).
1364      */
1365     ec_GFp_nistp_points_make_affine_internal(num,
1366                                              points,
1367                                              sizeof(felem),
1368                                              tmp_felems,
1369                                              (void (*)(void *))felem_one,
1370                                              felem_is_zero_int,
1371                                              (void (*)(void *, const void *))
1372                                              felem_assign,
1373                                              (void (*)(void *, const void *))
1374                                              felem_square_reduce, (void (*)
1375                                                                    (void *,
1376                                                                     const void
1377                                                                     *,
1378                                                                     const void
1379                                                                     *))
1380                                              felem_mul_reduce,
1381                                              (void (*)(void *, const void *))
1382                                              felem_inv,
1383                                              (void (*)(void *, const void *))
1384                                              felem_contract);
1385 }
1386
1387 /*
1388  * Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL
1389  * values Result is stored in r (r can equal one of the inputs).
1390  */
1391 int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
1392                                const BIGNUM *scalar, size_t num,
1393                                const EC_POINT *points[],
1394                                const BIGNUM *scalars[], BN_CTX *ctx)
1395 {
1396     int ret = 0;
1397     int j;
1398     unsigned i;
1399     int mixed = 0;
1400     BIGNUM *x, *y, *z, *tmp_scalar;
1401     felem_bytearray g_secret;
1402     felem_bytearray *secrets = NULL;
1403     felem (*pre_comp)[17][3] = NULL;
1404     felem *tmp_felems = NULL;
1405     felem_bytearray tmp;
1406     unsigned num_bytes;
1407     int have_pre_comp = 0;
1408     size_t num_points = num;
1409     felem x_in, y_in, z_in, x_out, y_out, z_out;
1410     NISTP224_PRE_COMP *pre = NULL;
1411     const felem(*g_pre_comp)[16][3] = NULL;
1412     EC_POINT *generator = NULL;
1413     const EC_POINT *p = NULL;
1414     const BIGNUM *p_scalar = NULL;
1415
1416     BN_CTX_start(ctx);
1417     x = BN_CTX_get(ctx);
1418     y = BN_CTX_get(ctx);
1419     z = BN_CTX_get(ctx);
1420     tmp_scalar = BN_CTX_get(ctx);
1421     if (tmp_scalar == NULL)
1422         goto err;
1423
1424     if (scalar != NULL) {
1425         pre = group->pre_comp.nistp224;
1426         if (pre)
1427             /* we have precomputation, try to use it */
1428             g_pre_comp = (const felem(*)[16][3])pre->g_pre_comp;
1429         else
1430             /* try to use the standard precomputation */
1431             g_pre_comp = &gmul[0];
1432         generator = EC_POINT_new(group);
1433         if (generator == NULL)
1434             goto err;
1435         /* get the generator from precomputation */
1436         if (!felem_to_BN(x, g_pre_comp[0][1][0]) ||
1437             !felem_to_BN(y, g_pre_comp[0][1][1]) ||
1438             !felem_to_BN(z, g_pre_comp[0][1][2])) {
1439             ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1440             goto err;
1441         }
1442         if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1443                                                       generator, x, y, z,
1444                                                       ctx))
1445             goto err;
1446         if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1447             /* precomputation matches generator */
1448             have_pre_comp = 1;
1449         else
1450             /*
1451              * we don't have valid precomputation: treat the generator as a
1452              * random point
1453              */
1454             num_points = num_points + 1;
1455     }
1456
1457     if (num_points > 0) {
1458         if (num_points >= 3) {
1459             /*
1460              * unless we precompute multiples for just one or two points,
1461              * converting those into affine form is time well spent
1462              */
1463             mixed = 1;
1464         }
1465         secrets = OPENSSL_zalloc(sizeof(*secrets) * num_points);
1466         pre_comp = OPENSSL_zalloc(sizeof(*pre_comp) * num_points);
1467         if (mixed)
1468             tmp_felems =
1469                 OPENSSL_malloc(sizeof(felem) * (num_points * 17 + 1));
1470         if ((secrets == NULL) || (pre_comp == NULL)
1471             || (mixed && (tmp_felems == NULL))) {
1472             ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1473             goto err;
1474         }
1475
1476         /*
1477          * we treat NULL scalars as 0, and NULL points as points at infinity,
1478          * i.e., they contribute nothing to the linear combination
1479          */
1480         for (i = 0; i < num_points; ++i) {
1481             if (i == num)
1482                 /* the generator */
1483             {
1484                 p = EC_GROUP_get0_generator(group);
1485                 p_scalar = scalar;
1486             } else
1487                 /* the i^th point */
1488             {
1489                 p = points[i];
1490                 p_scalar = scalars[i];
1491             }
1492             if ((p_scalar != NULL) && (p != NULL)) {
1493                 /* reduce scalar to 0 <= scalar < 2^224 */
1494                 if ((BN_num_bits(p_scalar) > 224)
1495                     || (BN_is_negative(p_scalar))) {
1496                     /*
1497                      * this is an unusual input, and we don't guarantee
1498                      * constant-timeness
1499                      */
1500                     if (!BN_nnmod(tmp_scalar, p_scalar, group->order, ctx)) {
1501                         ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1502                         goto err;
1503                     }
1504                     num_bytes = BN_bn2bin(tmp_scalar, tmp);
1505                 } else
1506                     num_bytes = BN_bn2bin(p_scalar, tmp);
1507                 flip_endian(secrets[i], tmp, num_bytes);
1508                 /* precompute multiples */
1509                 if ((!BN_to_felem(x_out, p->X)) ||
1510                     (!BN_to_felem(y_out, p->Y)) ||
1511                     (!BN_to_felem(z_out, p->Z)))
1512                     goto err;
1513                 felem_assign(pre_comp[i][1][0], x_out);
1514                 felem_assign(pre_comp[i][1][1], y_out);
1515                 felem_assign(pre_comp[i][1][2], z_out);
1516                 for (j = 2; j <= 16; ++j) {
1517                     if (j & 1) {
1518                         point_add(pre_comp[i][j][0], pre_comp[i][j][1],
1519                                   pre_comp[i][j][2], pre_comp[i][1][0],
1520                                   pre_comp[i][1][1], pre_comp[i][1][2], 0,
1521                                   pre_comp[i][j - 1][0],
1522                                   pre_comp[i][j - 1][1],
1523                                   pre_comp[i][j - 1][2]);
1524                     } else {
1525                         point_double(pre_comp[i][j][0], pre_comp[i][j][1],
1526                                      pre_comp[i][j][2], pre_comp[i][j / 2][0],
1527                                      pre_comp[i][j / 2][1],
1528                                      pre_comp[i][j / 2][2]);
1529                     }
1530                 }
1531             }
1532         }
1533         if (mixed)
1534             make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
1535     }
1536
1537     /* the scalar for the generator */
1538     if ((scalar != NULL) && (have_pre_comp)) {
1539         memset(g_secret, 0, sizeof(g_secret));
1540         /* reduce scalar to 0 <= scalar < 2^224 */
1541         if ((BN_num_bits(scalar) > 224) || (BN_is_negative(scalar))) {
1542             /*
1543              * this is an unusual input, and we don't guarantee
1544              * constant-timeness
1545              */
1546             if (!BN_nnmod(tmp_scalar, scalar, group->order, ctx)) {
1547                 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1548                 goto err;
1549             }
1550             num_bytes = BN_bn2bin(tmp_scalar, tmp);
1551         } else
1552             num_bytes = BN_bn2bin(scalar, tmp);
1553         flip_endian(g_secret, tmp, num_bytes);
1554         /* do the multiplication with generator precomputation */
1555         batch_mul(x_out, y_out, z_out,
1556                   (const felem_bytearray(*))secrets, num_points,
1557                   g_secret,
1558                   mixed, (const felem(*)[17][3])pre_comp, g_pre_comp);
1559     } else
1560         /* do the multiplication without generator precomputation */
1561         batch_mul(x_out, y_out, z_out,
1562                   (const felem_bytearray(*))secrets, num_points,
1563                   NULL, mixed, (const felem(*)[17][3])pre_comp, NULL);
1564     /* reduce the output to its unique minimal representation */
1565     felem_contract(x_in, x_out);
1566     felem_contract(y_in, y_out);
1567     felem_contract(z_in, z_out);
1568     if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
1569         (!felem_to_BN(z, z_in))) {
1570         ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1571         goto err;
1572     }
1573     ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
1574
1575  err:
1576     BN_CTX_end(ctx);
1577     EC_POINT_free(generator);
1578     OPENSSL_free(secrets);
1579     OPENSSL_free(pre_comp);
1580     OPENSSL_free(tmp_felems);
1581     return ret;
1582 }
1583
1584 int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
1585 {
1586     int ret = 0;
1587     NISTP224_PRE_COMP *pre = NULL;
1588     int i, j;
1589     BN_CTX *new_ctx = NULL;
1590     BIGNUM *x, *y;
1591     EC_POINT *generator = NULL;
1592     felem tmp_felems[32];
1593
1594     /* throw away old precomputation */
1595     EC_pre_comp_free(group);
1596     if (ctx == NULL)
1597         if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1598             return 0;
1599     BN_CTX_start(ctx);
1600     x = BN_CTX_get(ctx);
1601     y = BN_CTX_get(ctx);
1602     if (y == NULL)
1603         goto err;
1604     /* get the generator */
1605     if (group->generator == NULL)
1606         goto err;
1607     generator = EC_POINT_new(group);
1608     if (generator == NULL)
1609         goto err;
1610     BN_bin2bn(nistp224_curve_params[3], sizeof(felem_bytearray), x);
1611     BN_bin2bn(nistp224_curve_params[4], sizeof(felem_bytearray), y);
1612     if (!EC_POINT_set_affine_coordinates(group, generator, x, y, ctx))
1613         goto err;
1614     if ((pre = nistp224_pre_comp_new()) == NULL)
1615         goto err;
1616     /*
1617      * if the generator is the standard one, use built-in precomputation
1618      */
1619     if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) {
1620         memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
1621         goto done;
1622     }
1623     if ((!BN_to_felem(pre->g_pre_comp[0][1][0], group->generator->X)) ||
1624         (!BN_to_felem(pre->g_pre_comp[0][1][1], group->generator->Y)) ||
1625         (!BN_to_felem(pre->g_pre_comp[0][1][2], group->generator->Z)))
1626         goto err;
1627     /*
1628      * compute 2^56*G, 2^112*G, 2^168*G for the first table, 2^28*G, 2^84*G,
1629      * 2^140*G, 2^196*G for the second one
1630      */
1631     for (i = 1; i <= 8; i <<= 1) {
1632         point_double(pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1],
1633                      pre->g_pre_comp[1][i][2], pre->g_pre_comp[0][i][0],
1634                      pre->g_pre_comp[0][i][1], pre->g_pre_comp[0][i][2]);
1635         for (j = 0; j < 27; ++j) {
1636             point_double(pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1],
1637                          pre->g_pre_comp[1][i][2], pre->g_pre_comp[1][i][0],
1638                          pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1639         }
1640         if (i == 8)
1641             break;
1642         point_double(pre->g_pre_comp[0][2 * i][0],
1643                      pre->g_pre_comp[0][2 * i][1],
1644                      pre->g_pre_comp[0][2 * i][2], pre->g_pre_comp[1][i][0],
1645                      pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1646         for (j = 0; j < 27; ++j) {
1647             point_double(pre->g_pre_comp[0][2 * i][0],
1648                          pre->g_pre_comp[0][2 * i][1],
1649                          pre->g_pre_comp[0][2 * i][2],
1650                          pre->g_pre_comp[0][2 * i][0],
1651                          pre->g_pre_comp[0][2 * i][1],
1652                          pre->g_pre_comp[0][2 * i][2]);
1653         }
1654     }
1655     for (i = 0; i < 2; i++) {
1656         /* g_pre_comp[i][0] is the point at infinity */
1657         memset(pre->g_pre_comp[i][0], 0, sizeof(pre->g_pre_comp[i][0]));
1658         /* the remaining multiples */
1659         /* 2^56*G + 2^112*G resp. 2^84*G + 2^140*G */
1660         point_add(pre->g_pre_comp[i][6][0], pre->g_pre_comp[i][6][1],
1661                   pre->g_pre_comp[i][6][2], pre->g_pre_comp[i][4][0],
1662                   pre->g_pre_comp[i][4][1], pre->g_pre_comp[i][4][2],
1663                   0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1664                   pre->g_pre_comp[i][2][2]);
1665         /* 2^56*G + 2^168*G resp. 2^84*G + 2^196*G */
1666         point_add(pre->g_pre_comp[i][10][0], pre->g_pre_comp[i][10][1],
1667                   pre->g_pre_comp[i][10][2], pre->g_pre_comp[i][8][0],
1668                   pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1669                   0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1670                   pre->g_pre_comp[i][2][2]);
1671         /* 2^112*G + 2^168*G resp. 2^140*G + 2^196*G */
1672         point_add(pre->g_pre_comp[i][12][0], pre->g_pre_comp[i][12][1],
1673                   pre->g_pre_comp[i][12][2], pre->g_pre_comp[i][8][0],
1674                   pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1675                   0, pre->g_pre_comp[i][4][0], pre->g_pre_comp[i][4][1],
1676                   pre->g_pre_comp[i][4][2]);
1677         /*
1678          * 2^56*G + 2^112*G + 2^168*G resp. 2^84*G + 2^140*G + 2^196*G
1679          */
1680         point_add(pre->g_pre_comp[i][14][0], pre->g_pre_comp[i][14][1],
1681                   pre->g_pre_comp[i][14][2], pre->g_pre_comp[i][12][0],
1682                   pre->g_pre_comp[i][12][1], pre->g_pre_comp[i][12][2],
1683                   0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1684                   pre->g_pre_comp[i][2][2]);
1685         for (j = 1; j < 8; ++j) {
1686             /* odd multiples: add G resp. 2^28*G */
1687             point_add(pre->g_pre_comp[i][2 * j + 1][0],
1688                       pre->g_pre_comp[i][2 * j + 1][1],
1689                       pre->g_pre_comp[i][2 * j + 1][2],
1690                       pre->g_pre_comp[i][2 * j][0],
1691                       pre->g_pre_comp[i][2 * j][1],
1692                       pre->g_pre_comp[i][2 * j][2], 0,
1693                       pre->g_pre_comp[i][1][0], pre->g_pre_comp[i][1][1],
1694                       pre->g_pre_comp[i][1][2]);
1695         }
1696     }
1697     make_points_affine(31, &(pre->g_pre_comp[0][1]), tmp_felems);
1698
1699  done:
1700     SETPRECOMP(group, nistp224, pre);
1701     pre = NULL;
1702     ret = 1;
1703  err:
1704     BN_CTX_end(ctx);
1705     EC_POINT_free(generator);
1706     BN_CTX_free(new_ctx);
1707     EC_nistp224_pre_comp_free(pre);
1708     return ret;
1709 }
1710
1711 int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
1712 {
1713     return HAVEPRECOMP(group, nistp224);
1714 }
1715
1716 #endif