1 /* 2 * Copyright 2010-2019 The OpenSSL Project Authors. All Rights Reserved. 3 * 4 * Licensed under the OpenSSL license (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 1717