2 * Copyright 2010-2019 The OpenSSL Project Authors. All Rights Reserved.
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
10 /* Copyright 2011 Google Inc.
12 * Licensed under the Apache License, Version 2.0 (the "License");
14 * you may not use this file except in compliance with the License.
15 * You may obtain a copy of the License at
17 * http://www.apache.org/licenses/LICENSE-2.0
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.
27 * A 64-bit implementation of the NIST P-224 elliptic curve point multiplication
29 * Inspired by Daniel J. Bernstein's public domain nistp224 implementation
30 * and Adam Langley's public domain 64-bit C implementation of curve25519
33 #include <openssl/opensslconf.h>
34 #ifdef OPENSSL_NO_EC_NISTP_64_GCC_128
35 NON_EMPTY_TRANSLATION_UNIT
40 # include <openssl/err.h>
43 # if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 1))
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
48 # error "Need GCC 3.1 or later to define type uint128_t"
54 /******************************************************************************/
56 * INTERNAL REPRESENTATION OF FIELD ELEMENTS
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
74 typedef uint64_t limb;
75 typedef uint128_t widelimb;
77 typedef limb felem[4];
78 typedef widelimb widefelem[7];
81 * Field element represented as a byte arrary. 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.
85 typedef u8 felem_bytearray[28];
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}
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.
111 * The table has 2 * 16 elements, starting with the following:
112 * index | bits | point
113 * ------+---------+------------------------------
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.
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.
136 static const felem gmul[2][16][3] = {
140 {{0x3280d6115c1d21, 0xc1d356c2112234, 0x7f321390b94a03, 0xb70e0cbd6bb4bf},
141 {0xd5819985007e34, 0x75a05a07476444, 0xfb4c22dfe6cd43, 0xbd376388b5f723},
143 {{0xfd9675666ebbe9, 0xbca7664d40ce5e, 0x2242df8d8a2a43, 0x1f49bbb0f99bc5},
144 {0x29e0b892dc9c43, 0xece8608436e662, 0xdc858f185310d0, 0x9812dd4eb8d321},
146 {{0x6d3e678d5d8eb8, 0x559eed1cb362f1, 0x16e9a3bbce8a3f, 0xeedcccd8c2a748},
147 {0xf19f90ed50266d, 0xabf2b4bf65f9df, 0x313865468fafec, 0x5cb379ba910a17},
149 {{0x0641966cab26e3, 0x91fb2991fab0a0, 0xefec27a4e13a0b, 0x0499aa8a5f8ebe},
150 {0x7510407766af5d, 0x84d929610d5450, 0x81d77aae82f706, 0x6916f6d4338c5b},
152 {{0xea95ac3b1f15c6, 0x086000905e82d4, 0xdd323ae4d1c8b1, 0x932b56be7685a3},
153 {0x9ef93dea25dbbf, 0x41665960f390f0, 0xfdec76dbe2a8a7, 0x523e80f019062a},
155 {{0x822fdd26732c73, 0xa01c83531b5d0f, 0x363f37347c1ba4, 0xc391b45c84725c},
156 {0xbbd5e1b2d6ad24, 0xddfbcde19dfaec, 0xc393da7e222a7f, 0x1efb7890ede244},
158 {{0x4c9e90ca217da1, 0xd11beca79159bb, 0xff8d33c2c98b7c, 0x2610b39409f849},
159 {0x44d1352ac64da0, 0xcdbb7b2c46b4fb, 0x966c079b753c89, 0xfe67e4e820b112},
161 {{0xe28cae2df5312d, 0xc71b61d16f5c6e, 0x79b7619a3e7c4c, 0x05c73240899b47},
162 {0x9f7f6382c73e3a, 0x18615165c56bda, 0x641fab2116fd56, 0x72855882b08394},
164 {{0x0469182f161c09, 0x74a98ca8d00fb5, 0xb89da93489a3e0, 0x41c98768fb0c1d},
165 {0xe5ea05fb32da81, 0x3dce9ffbca6855, 0x1cfe2d3fbf59e6, 0x0e5e03408738a7},
167 {{0xdab22b2333e87f, 0x4430137a5dd2f6, 0xe03ab9f738beb8, 0xcb0c5d0dc34f24},
168 {0x764a7df0c8fda5, 0x185ba5c3fa2044, 0x9281d688bcbe50, 0xc40331df893881},
170 {{0xb89530796f0f60, 0xade92bd26909a3, 0x1a0c83fb4884da, 0x1765bf22a5a984},
171 {0x772a9ee75db09e, 0x23bc6c67cec16f, 0x4c1edba8b14e2f, 0xe2a215d9611369},
173 {{0x571e509fb5efb3, 0xade88696410552, 0xc8ae85fada74fe, 0x6c7e4be83bbde3},
174 {0xff9f51160f4652, 0xb47ce2495a6539, 0xa2946c53b582f4, 0x286d2db3ee9a60},
176 {{0x40bbd5081a44af, 0x0995183b13926c, 0xbcefba6f47f6d0, 0x215619e9cc0057},
177 {0x8bc94d3b0df45e, 0xf11c54a3694f6f, 0x8631b93cdfe8b5, 0xe7e3f4b0982db9},
179 {{0xb17048ab3e1c7b, 0xac38f36ff8a1d8, 0x1c29819435d2c6, 0xc813132f4c07e9},
180 {0x2891425503b11f, 0x08781030579fea, 0xf5426ba5cc9674, 0x1e28ebf18562bc},
182 {{0x9f31997cc864eb, 0x06cd91d28b5e4c, 0xff17036691a973, 0xf1aef351497c58},
183 {0xdd1f2d600564ff, 0xdead073b1402db, 0x74a684435bd693, 0xeea7471f962558},
188 {{0x9665266dddf554, 0x9613d78b60ef2d, 0xce27a34cdba417, 0xd35ab74d6afc31},
189 {0x85ccdd22deb15e, 0x2137e5783a6aab, 0xa141cffd8c93c6, 0x355a1830e90f2d},
191 {{0x1a494eadaade65, 0xd6da4da77fe53c, 0xe7992996abec86, 0x65c3553c6090e3},
192 {0xfa610b1fb09346, 0xf1c6540b8a4aaf, 0xc51a13ccd3cbab, 0x02995b1b18c28a},
194 {{0x7874568e7295ef, 0x86b419fbe38d04, 0xdc0690a7550d9a, 0xd3966a44beac33},
195 {0x2b7280ec29132f, 0xbeaa3b6a032df3, 0xdc7dd88ae41200, 0xd25e2513e3a100},
197 {{0x924857eb2efafd, 0xac2bce41223190, 0x8edaa1445553fc, 0x825800fd3562d5},
198 {0x8d79148ea96621, 0x23a01c3dd9ed8d, 0xaf8b219f9416b5, 0xd8db0cc277daea},
200 {{0x76a9c3b1a700f0, 0xe9acd29bc7e691, 0x69212d1a6b0327, 0x6322e97fe154be},
201 {0x469fc5465d62aa, 0x8d41ed18883b05, 0x1f8eae66c52b88, 0xe4fcbe9325be51},
203 {{0x825fdf583cac16, 0x020b857c7b023a, 0x683c17744b0165, 0x14ffd0a2daf2f1},
204 {0x323b36184218f9, 0x4944ec4e3b47d4, 0xc15b3080841acf, 0x0bced4b01a28bb},
206 {{0x92ac22230df5c4, 0x52f33b4063eda8, 0xcb3f19870c0c93, 0x40064f2ba65233},
207 {0xfe16f0924f8992, 0x012da25af5b517, 0x1a57bb24f723a6, 0x06f8bc76760def},
209 {{0x4a7084f7817cb9, 0xbcab0738ee9a78, 0x3ec11e11d9c326, 0xdc0fe90e0f1aae},
210 {0xcf639ea5f98390, 0x5c350aa22ffb74, 0x9afae98a4047b7, 0x956ec2d617fc45},
212 {{0x4306d648c1be6a, 0x9247cd8bc9a462, 0xf5595e377d2f2e, 0xbd1c3caff1a52e},
213 {0x045e14472409d0, 0x29f3e17078f773, 0x745a602b2d4f7d, 0x191837685cdfbb},
215 {{0x5b6ee254a8cb79, 0x4953433f5e7026, 0xe21faeb1d1def4, 0xc4c225785c09de},
216 {0x307ce7bba1e518, 0x31b125b1036db8, 0x47e91868839e8f, 0xc765866e33b9f3},
218 {{0x3bfece24f96906, 0x4794da641e5093, 0xde5df64f95db26, 0x297ecd89714b05},
219 {0x701bd3ebb2c3aa, 0x7073b4f53cb1d5, 0x13c5665658af16, 0x9895089d66fe58},
221 {{0x0fef05f78c4790, 0x2d773633b05d2e, 0x94229c3a951c94, 0xbbbd70df4911bb},
222 {0xb2c6963d2c1168, 0x105f47a72b0d73, 0x9fdf6111614080, 0x7b7e94b39e67b0},
224 {{0xad1a7d6efbe2b3, 0xf012482c0da69d, 0x6b3bdf12438345, 0x40d7558d7aa4d9},
225 {0x8a09fffb5c6d3d, 0x9a356e5d9ffd38, 0x5973f15f4f9b1c, 0xdcd5f59f63c3ea},
227 {{0xacf39f4c5ca7ab, 0x4c8071cc5fd737, 0xc64e3602cd1184, 0x0acd4644c9abba},
228 {0x6c011a36d8bf6e, 0xfecd87ba24e32a, 0x19f6f56574fad8, 0x050b204ced9405},
230 {{0xed4f1cae7d9a96, 0x5ceef7ad94c40a, 0x778e4a3bf3ef9b, 0x7405783dc3b55e},
231 {0x32477c61b6e8c6, 0xb46a97570f018b, 0x91176d0a7e95d1, 0x3df90fbc4c7d0e},
235 /* Precomputation for the group generator. */
236 struct nistp224_pre_comp_st {
237 felem g_pre_comp[2][16][3];
242 const EC_METHOD *EC_GFp_nistp224_method(void)
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 */ ,
270 ec_GFp_simple_invert,
271 ec_GFp_simple_is_at_infinity,
272 ec_GFp_simple_is_on_curve,
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,
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,
289 ec_key_simple_generate_key,
290 ec_key_simple_check_key,
291 ec_key_simple_generate_public_key,
294 ecdh_simple_compute_key,
295 0 /* blind_coordinates */
302 * Helper functions to convert field elements to/from internal representation
304 static void bin28_to_felem(felem out, const u8 in[28])
306 out[0] = *((const uint64_t *)(in)) & 0x00ffffffffffffff;
307 out[1] = (*((const uint64_t *)(in + 7))) & 0x00ffffffffffffff;
308 out[2] = (*((const uint64_t *)(in + 14))) & 0x00ffffffffffffff;
309 out[3] = (*((const uint64_t *)(in+20))) >> 8;
312 static void felem_to_bin28(u8 out[28], const felem in)
315 for (i = 0; i < 7; ++i) {
316 out[i] = in[0] >> (8 * i);
317 out[i + 7] = in[1] >> (8 * i);
318 out[i + 14] = in[2] >> (8 * i);
319 out[i + 21] = in[3] >> (8 * i);
323 /* From OpenSSL BIGNUM to internal representation */
324 static int BN_to_felem(felem out, const BIGNUM *bn)
326 felem_bytearray b_out;
329 if (BN_is_negative(bn)) {
330 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
333 num_bytes = BN_bn2lebinpad(bn, b_out, sizeof(b_out));
335 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
338 bin28_to_felem(out, b_out);
342 /* From internal representation to OpenSSL BIGNUM */
343 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
345 felem_bytearray b_out;
346 felem_to_bin28(b_out, in);
347 return BN_lebin2bn(b_out, sizeof(b_out), out);
350 /******************************************************************************/
354 * Field operations, using the internal representation of field elements.
355 * NB! These operations are specific to our point multiplication and cannot be
356 * expected to be correct in general - e.g., multiplication with a large scalar
357 * will cause an overflow.
361 static void felem_one(felem out)
369 static void felem_assign(felem out, const felem in)
377 /* Sum two field elements: out += in */
378 static void felem_sum(felem out, const felem in)
386 /* Get negative value: out = -in */
387 /* Assumes in[i] < 2^57 */
388 static void felem_neg(felem out, const felem in)
390 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
391 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
392 static const limb two58m42m2 = (((limb) 1) << 58) -
393 (((limb) 1) << 42) - (((limb) 1) << 2);
395 /* Set to 0 mod 2^224-2^96+1 to ensure out > in */
396 out[0] = two58p2 - in[0];
397 out[1] = two58m42m2 - in[1];
398 out[2] = two58m2 - in[2];
399 out[3] = two58m2 - in[3];
402 /* Subtract field elements: out -= in */
403 /* Assumes in[i] < 2^57 */
404 static void felem_diff(felem out, const felem in)
406 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
407 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
408 static const limb two58m42m2 = (((limb) 1) << 58) -
409 (((limb) 1) << 42) - (((limb) 1) << 2);
411 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
413 out[1] += two58m42m2;
423 /* Subtract in unreduced 128-bit mode: out -= in */
424 /* Assumes in[i] < 2^119 */
425 static void widefelem_diff(widefelem out, const widefelem in)
427 static const widelimb two120 = ((widelimb) 1) << 120;
428 static const widelimb two120m64 = (((widelimb) 1) << 120) -
429 (((widelimb) 1) << 64);
430 static const widelimb two120m104m64 = (((widelimb) 1) << 120) -
431 (((widelimb) 1) << 104) - (((widelimb) 1) << 64);
433 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
438 out[4] += two120m104m64;
451 /* Subtract in mixed mode: out128 -= in64 */
453 static void felem_diff_128_64(widefelem out, const felem in)
455 static const widelimb two64p8 = (((widelimb) 1) << 64) +
456 (((widelimb) 1) << 8);
457 static const widelimb two64m8 = (((widelimb) 1) << 64) -
458 (((widelimb) 1) << 8);
459 static const widelimb two64m48m8 = (((widelimb) 1) << 64) -
460 (((widelimb) 1) << 48) - (((widelimb) 1) << 8);
462 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
464 out[1] += two64m48m8;
475 * Multiply a field element by a scalar: out = out * scalar The scalars we
476 * actually use are small, so results fit without overflow
478 static void felem_scalar(felem out, const limb scalar)
487 * Multiply an unreduced field element by a scalar: out = out * scalar The
488 * scalars we actually use are small, so results fit without overflow
490 static void widefelem_scalar(widefelem out, const widelimb scalar)
501 /* Square a field element: out = in^2 */
502 static void felem_square(widefelem out, const felem in)
504 limb tmp0, tmp1, tmp2;
508 out[0] = ((widelimb) in[0]) * in[0];
509 out[1] = ((widelimb) in[0]) * tmp1;
510 out[2] = ((widelimb) in[0]) * tmp2 + ((widelimb) in[1]) * in[1];
511 out[3] = ((widelimb) in[3]) * tmp0 + ((widelimb) in[1]) * tmp2;
512 out[4] = ((widelimb) in[3]) * tmp1 + ((widelimb) in[2]) * in[2];
513 out[5] = ((widelimb) in[3]) * tmp2;
514 out[6] = ((widelimb) in[3]) * in[3];
517 /* Multiply two field elements: out = in1 * in2 */
518 static void felem_mul(widefelem out, const felem in1, const felem in2)
520 out[0] = ((widelimb) in1[0]) * in2[0];
521 out[1] = ((widelimb) in1[0]) * in2[1] + ((widelimb) in1[1]) * in2[0];
522 out[2] = ((widelimb) in1[0]) * in2[2] + ((widelimb) in1[1]) * in2[1] +
523 ((widelimb) in1[2]) * in2[0];
524 out[3] = ((widelimb) in1[0]) * in2[3] + ((widelimb) in1[1]) * in2[2] +
525 ((widelimb) in1[2]) * in2[1] + ((widelimb) in1[3]) * in2[0];
526 out[4] = ((widelimb) in1[1]) * in2[3] + ((widelimb) in1[2]) * in2[2] +
527 ((widelimb) in1[3]) * in2[1];
528 out[5] = ((widelimb) in1[2]) * in2[3] + ((widelimb) in1[3]) * in2[2];
529 out[6] = ((widelimb) in1[3]) * in2[3];
533 * Reduce seven 128-bit coefficients to four 64-bit coefficients.
534 * Requires in[i] < 2^126,
535 * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] <= 2^56 + 2^16 */
536 static void felem_reduce(felem out, const widefelem in)
538 static const widelimb two127p15 = (((widelimb) 1) << 127) +
539 (((widelimb) 1) << 15);
540 static const widelimb two127m71 = (((widelimb) 1) << 127) -
541 (((widelimb) 1) << 71);
542 static const widelimb two127m71m55 = (((widelimb) 1) << 127) -
543 (((widelimb) 1) << 71) - (((widelimb) 1) << 55);
546 /* Add 0 mod 2^224-2^96+1 to ensure all differences are positive */
547 output[0] = in[0] + two127p15;
548 output[1] = in[1] + two127m71m55;
549 output[2] = in[2] + two127m71;
553 /* Eliminate in[4], in[5], in[6] */
554 output[4] += in[6] >> 16;
555 output[3] += (in[6] & 0xffff) << 40;
558 output[3] += in[5] >> 16;
559 output[2] += (in[5] & 0xffff) << 40;
562 output[2] += output[4] >> 16;
563 output[1] += (output[4] & 0xffff) << 40;
564 output[0] -= output[4];
566 /* Carry 2 -> 3 -> 4 */
567 output[3] += output[2] >> 56;
568 output[2] &= 0x00ffffffffffffff;
570 output[4] = output[3] >> 56;
571 output[3] &= 0x00ffffffffffffff;
573 /* Now output[2] < 2^56, output[3] < 2^56, output[4] < 2^72 */
575 /* Eliminate output[4] */
576 output[2] += output[4] >> 16;
577 /* output[2] < 2^56 + 2^56 = 2^57 */
578 output[1] += (output[4] & 0xffff) << 40;
579 output[0] -= output[4];
581 /* Carry 0 -> 1 -> 2 -> 3 */
582 output[1] += output[0] >> 56;
583 out[0] = output[0] & 0x00ffffffffffffff;
585 output[2] += output[1] >> 56;
586 /* output[2] < 2^57 + 2^72 */
587 out[1] = output[1] & 0x00ffffffffffffff;
588 output[3] += output[2] >> 56;
589 /* output[3] <= 2^56 + 2^16 */
590 out[2] = output[2] & 0x00ffffffffffffff;
593 * out[0] < 2^56, out[1] < 2^56, out[2] < 2^56,
594 * out[3] <= 2^56 + 2^16 (due to final carry),
600 static void felem_square_reduce(felem out, const felem in)
603 felem_square(tmp, in);
604 felem_reduce(out, tmp);
607 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
610 felem_mul(tmp, in1, in2);
611 felem_reduce(out, tmp);
615 * Reduce to unique minimal representation. Requires 0 <= in < 2*p (always
616 * call felem_reduce first)
618 static void felem_contract(felem out, const felem in)
620 static const int64_t two56 = ((limb) 1) << 56;
621 /* 0 <= in < 2*p, p = 2^224 - 2^96 + 1 */
622 /* if in > p , reduce in = in - 2^224 + 2^96 - 1 */
628 /* Case 1: a = 1 iff in >= 2^224 */
632 tmp[3] &= 0x00ffffffffffffff;
634 * Case 2: a = 0 iff p <= in < 2^224, i.e., the high 128 bits are all 1
635 * and the lower part is non-zero
637 a = ((in[3] & in[2] & (in[1] | 0x000000ffffffffff)) + 1) |
638 (((int64_t) (in[0] + (in[1] & 0x000000ffffffffff)) - 1) >> 63);
639 a &= 0x00ffffffffffffff;
640 /* turn a into an all-one mask (if a = 0) or an all-zero mask */
642 /* subtract 2^224 - 2^96 + 1 if a is all-one */
643 tmp[3] &= a ^ 0xffffffffffffffff;
644 tmp[2] &= a ^ 0xffffffffffffffff;
645 tmp[1] &= (a ^ 0xffffffffffffffff) | 0x000000ffffffffff;
649 * eliminate negative coefficients: if tmp[0] is negative, tmp[1] must be
650 * non-zero, so we only need one step
656 /* carry 1 -> 2 -> 3 */
657 tmp[2] += tmp[1] >> 56;
658 tmp[1] &= 0x00ffffffffffffff;
660 tmp[3] += tmp[2] >> 56;
661 tmp[2] &= 0x00ffffffffffffff;
663 /* Now 0 <= out < p */
671 * Zero-check: returns 1 if input is 0, and 0 otherwise. We know that field
672 * elements are reduced to in < 2^225, so we only need to check three cases:
673 * 0, 2^224 - 2^96 + 1, and 2^225 - 2^97 + 2
675 static limb felem_is_zero(const felem in)
677 limb zero, two224m96p1, two225m97p2;
679 zero = in[0] | in[1] | in[2] | in[3];
680 zero = (((int64_t) (zero) - 1) >> 63) & 1;
681 two224m96p1 = (in[0] ^ 1) | (in[1] ^ 0x00ffff0000000000)
682 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x00ffffffffffffff);
683 two224m96p1 = (((int64_t) (two224m96p1) - 1) >> 63) & 1;
684 two225m97p2 = (in[0] ^ 2) | (in[1] ^ 0x00fffe0000000000)
685 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x01ffffffffffffff);
686 two225m97p2 = (((int64_t) (two225m97p2) - 1) >> 63) & 1;
687 return (zero | two224m96p1 | two225m97p2);
690 static int felem_is_zero_int(const void *in)
692 return (int)(felem_is_zero(in) & ((limb) 1));
695 /* Invert a field element */
696 /* Computation chain copied from djb's code */
697 static void felem_inv(felem out, const felem in)
699 felem ftmp, ftmp2, ftmp3, ftmp4;
703 felem_square(tmp, in);
704 felem_reduce(ftmp, tmp); /* 2 */
705 felem_mul(tmp, in, ftmp);
706 felem_reduce(ftmp, tmp); /* 2^2 - 1 */
707 felem_square(tmp, ftmp);
708 felem_reduce(ftmp, tmp); /* 2^3 - 2 */
709 felem_mul(tmp, in, ftmp);
710 felem_reduce(ftmp, tmp); /* 2^3 - 1 */
711 felem_square(tmp, ftmp);
712 felem_reduce(ftmp2, tmp); /* 2^4 - 2 */
713 felem_square(tmp, ftmp2);
714 felem_reduce(ftmp2, tmp); /* 2^5 - 4 */
715 felem_square(tmp, ftmp2);
716 felem_reduce(ftmp2, tmp); /* 2^6 - 8 */
717 felem_mul(tmp, ftmp2, ftmp);
718 felem_reduce(ftmp, tmp); /* 2^6 - 1 */
719 felem_square(tmp, ftmp);
720 felem_reduce(ftmp2, tmp); /* 2^7 - 2 */
721 for (i = 0; i < 5; ++i) { /* 2^12 - 2^6 */
722 felem_square(tmp, ftmp2);
723 felem_reduce(ftmp2, tmp);
725 felem_mul(tmp, ftmp2, ftmp);
726 felem_reduce(ftmp2, tmp); /* 2^12 - 1 */
727 felem_square(tmp, ftmp2);
728 felem_reduce(ftmp3, tmp); /* 2^13 - 2 */
729 for (i = 0; i < 11; ++i) { /* 2^24 - 2^12 */
730 felem_square(tmp, ftmp3);
731 felem_reduce(ftmp3, tmp);
733 felem_mul(tmp, ftmp3, ftmp2);
734 felem_reduce(ftmp2, tmp); /* 2^24 - 1 */
735 felem_square(tmp, ftmp2);
736 felem_reduce(ftmp3, tmp); /* 2^25 - 2 */
737 for (i = 0; i < 23; ++i) { /* 2^48 - 2^24 */
738 felem_square(tmp, ftmp3);
739 felem_reduce(ftmp3, tmp);
741 felem_mul(tmp, ftmp3, ftmp2);
742 felem_reduce(ftmp3, tmp); /* 2^48 - 1 */
743 felem_square(tmp, ftmp3);
744 felem_reduce(ftmp4, tmp); /* 2^49 - 2 */
745 for (i = 0; i < 47; ++i) { /* 2^96 - 2^48 */
746 felem_square(tmp, ftmp4);
747 felem_reduce(ftmp4, tmp);
749 felem_mul(tmp, ftmp3, ftmp4);
750 felem_reduce(ftmp3, tmp); /* 2^96 - 1 */
751 felem_square(tmp, ftmp3);
752 felem_reduce(ftmp4, tmp); /* 2^97 - 2 */
753 for (i = 0; i < 23; ++i) { /* 2^120 - 2^24 */
754 felem_square(tmp, ftmp4);
755 felem_reduce(ftmp4, tmp);
757 felem_mul(tmp, ftmp2, ftmp4);
758 felem_reduce(ftmp2, tmp); /* 2^120 - 1 */
759 for (i = 0; i < 6; ++i) { /* 2^126 - 2^6 */
760 felem_square(tmp, ftmp2);
761 felem_reduce(ftmp2, tmp);
763 felem_mul(tmp, ftmp2, ftmp);
764 felem_reduce(ftmp, tmp); /* 2^126 - 1 */
765 felem_square(tmp, ftmp);
766 felem_reduce(ftmp, tmp); /* 2^127 - 2 */
767 felem_mul(tmp, ftmp, in);
768 felem_reduce(ftmp, tmp); /* 2^127 - 1 */
769 for (i = 0; i < 97; ++i) { /* 2^224 - 2^97 */
770 felem_square(tmp, ftmp);
771 felem_reduce(ftmp, tmp);
773 felem_mul(tmp, ftmp, ftmp3);
774 felem_reduce(out, tmp); /* 2^224 - 2^96 - 1 */
778 * Copy in constant time: if icopy == 1, copy in to out, if icopy == 0, copy
781 static void copy_conditional(felem out, const felem in, limb icopy)
785 * icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one
787 const limb copy = -icopy;
788 for (i = 0; i < 4; ++i) {
789 const limb tmp = copy & (in[i] ^ out[i]);
794 /******************************************************************************/
796 * ELLIPTIC CURVE POINT OPERATIONS
798 * Points are represented in Jacobian projective coordinates:
799 * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
800 * or to the point at infinity if Z == 0.
805 * Double an elliptic curve point:
806 * (X', Y', Z') = 2 * (X, Y, Z), where
807 * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
808 * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^2
809 * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
810 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
811 * while x_out == y_in is not (maybe this works, but it's not tested).
814 point_double(felem x_out, felem y_out, felem z_out,
815 const felem x_in, const felem y_in, const felem z_in)
818 felem delta, gamma, beta, alpha, ftmp, ftmp2;
820 felem_assign(ftmp, x_in);
821 felem_assign(ftmp2, x_in);
824 felem_square(tmp, z_in);
825 felem_reduce(delta, tmp);
828 felem_square(tmp, y_in);
829 felem_reduce(gamma, tmp);
832 felem_mul(tmp, x_in, gamma);
833 felem_reduce(beta, tmp);
835 /* alpha = 3*(x-delta)*(x+delta) */
836 felem_diff(ftmp, delta);
837 /* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
838 felem_sum(ftmp2, delta);
839 /* ftmp2[i] < 2^57 + 2^57 = 2^58 */
840 felem_scalar(ftmp2, 3);
841 /* ftmp2[i] < 3 * 2^58 < 2^60 */
842 felem_mul(tmp, ftmp, ftmp2);
843 /* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
844 felem_reduce(alpha, tmp);
846 /* x' = alpha^2 - 8*beta */
847 felem_square(tmp, alpha);
848 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
849 felem_assign(ftmp, beta);
850 felem_scalar(ftmp, 8);
851 /* ftmp[i] < 8 * 2^57 = 2^60 */
852 felem_diff_128_64(tmp, ftmp);
853 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
854 felem_reduce(x_out, tmp);
856 /* z' = (y + z)^2 - gamma - delta */
857 felem_sum(delta, gamma);
858 /* delta[i] < 2^57 + 2^57 = 2^58 */
859 felem_assign(ftmp, y_in);
860 felem_sum(ftmp, z_in);
861 /* ftmp[i] < 2^57 + 2^57 = 2^58 */
862 felem_square(tmp, ftmp);
863 /* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
864 felem_diff_128_64(tmp, delta);
865 /* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
866 felem_reduce(z_out, tmp);
868 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
869 felem_scalar(beta, 4);
870 /* beta[i] < 4 * 2^57 = 2^59 */
871 felem_diff(beta, x_out);
872 /* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
873 felem_mul(tmp, alpha, beta);
874 /* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
875 felem_square(tmp2, gamma);
876 /* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
877 widefelem_scalar(tmp2, 8);
878 /* tmp2[i] < 8 * 2^116 = 2^119 */
879 widefelem_diff(tmp, tmp2);
880 /* tmp[i] < 2^119 + 2^120 < 2^121 */
881 felem_reduce(y_out, tmp);
885 * Add two elliptic curve points:
886 * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
887 * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
888 * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
889 * 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) -
890 * Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
891 * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2)
893 * This runs faster if 'mixed' is set, which requires Z_2 = 1 or Z_2 = 0.
897 * This function is not entirely constant-time: it includes a branch for
898 * checking whether the two input points are equal, (while not equal to the
899 * point at infinity). This case never happens during single point
900 * multiplication, so there is no timing leak for ECDH or ECDSA signing.
902 static void point_add(felem x3, felem y3, felem z3,
903 const felem x1, const felem y1, const felem z1,
904 const int mixed, const felem x2, const felem y2,
907 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, x_out, y_out, z_out;
909 limb z1_is_zero, z2_is_zero, x_equal, y_equal;
913 felem_square(tmp, z2);
914 felem_reduce(ftmp2, tmp);
917 felem_mul(tmp, ftmp2, z2);
918 felem_reduce(ftmp4, tmp);
920 /* ftmp4 = z2^3*y1 */
921 felem_mul(tmp2, ftmp4, y1);
922 felem_reduce(ftmp4, tmp2);
924 /* ftmp2 = z2^2*x1 */
925 felem_mul(tmp2, ftmp2, x1);
926 felem_reduce(ftmp2, tmp2);
929 * We'll assume z2 = 1 (special case z2 = 0 is handled later)
932 /* ftmp4 = z2^3*y1 */
933 felem_assign(ftmp4, y1);
935 /* ftmp2 = z2^2*x1 */
936 felem_assign(ftmp2, x1);
940 felem_square(tmp, z1);
941 felem_reduce(ftmp, tmp);
944 felem_mul(tmp, ftmp, z1);
945 felem_reduce(ftmp3, tmp);
948 felem_mul(tmp, ftmp3, y2);
949 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
951 /* ftmp3 = z1^3*y2 - z2^3*y1 */
952 felem_diff_128_64(tmp, ftmp4);
953 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
954 felem_reduce(ftmp3, tmp);
957 felem_mul(tmp, ftmp, x2);
958 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
960 /* ftmp = z1^2*x2 - z2^2*x1 */
961 felem_diff_128_64(tmp, ftmp2);
962 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
963 felem_reduce(ftmp, tmp);
966 * the formulae are incorrect if the points are equal so we check for
967 * this and do doubling if this happens
969 x_equal = felem_is_zero(ftmp);
970 y_equal = felem_is_zero(ftmp3);
971 z1_is_zero = felem_is_zero(z1);
972 z2_is_zero = felem_is_zero(z2);
973 /* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
974 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero) {
975 point_double(x3, y3, z3, x1, y1, z1);
981 felem_mul(tmp, z1, z2);
982 felem_reduce(ftmp5, tmp);
984 /* special case z2 = 0 is handled later */
985 felem_assign(ftmp5, z1);
988 /* z_out = (z1^2*x2 - z2^2*x1)*(z1*z2) */
989 felem_mul(tmp, ftmp, ftmp5);
990 felem_reduce(z_out, tmp);
992 /* ftmp = (z1^2*x2 - z2^2*x1)^2 */
993 felem_assign(ftmp5, ftmp);
994 felem_square(tmp, ftmp);
995 felem_reduce(ftmp, tmp);
997 /* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
998 felem_mul(tmp, ftmp, ftmp5);
999 felem_reduce(ftmp5, tmp);
1001 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
1002 felem_mul(tmp, ftmp2, ftmp);
1003 felem_reduce(ftmp2, tmp);
1005 /* tmp = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
1006 felem_mul(tmp, ftmp4, ftmp5);
1007 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
1009 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
1010 felem_square(tmp2, ftmp3);
1011 /* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
1013 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
1014 felem_diff_128_64(tmp2, ftmp5);
1015 /* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
1017 /* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
1018 felem_assign(ftmp5, ftmp2);
1019 felem_scalar(ftmp5, 2);
1020 /* ftmp5[i] < 2 * 2^57 = 2^58 */
1023 * x_out = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
1024 * 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2
1026 felem_diff_128_64(tmp2, ftmp5);
1027 /* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
1028 felem_reduce(x_out, tmp2);
1030 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out */
1031 felem_diff(ftmp2, x_out);
1032 /* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
1035 * tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out)
1037 felem_mul(tmp2, ftmp3, ftmp2);
1038 /* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
1041 * y_out = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) -
1042 * z2^3*y1*(z1^2*x2 - z2^2*x1)^3
1044 widefelem_diff(tmp2, tmp);
1045 /* tmp2[i] < 2^118 + 2^120 < 2^121 */
1046 felem_reduce(y_out, tmp2);
1049 * the result (x_out, y_out, z_out) is incorrect if one of the inputs is
1050 * the point at infinity, so we need to check for this separately
1054 * if point 1 is at infinity, copy point 2 to output, and vice versa
1056 copy_conditional(x_out, x2, z1_is_zero);
1057 copy_conditional(x_out, x1, z2_is_zero);
1058 copy_conditional(y_out, y2, z1_is_zero);
1059 copy_conditional(y_out, y1, z2_is_zero);
1060 copy_conditional(z_out, z2, z1_is_zero);
1061 copy_conditional(z_out, z1, z2_is_zero);
1062 felem_assign(x3, x_out);
1063 felem_assign(y3, y_out);
1064 felem_assign(z3, z_out);
1068 * select_point selects the |idx|th point from a precomputation table and
1070 * The pre_comp array argument should be size of |size| argument
1072 static void select_point(const u64 idx, unsigned int size,
1073 const felem pre_comp[][3], felem out[3])
1076 limb *outlimbs = &out[0][0];
1078 memset(out, 0, sizeof(*out) * 3);
1079 for (i = 0; i < size; i++) {
1080 const limb *inlimbs = &pre_comp[i][0][0];
1087 for (j = 0; j < 4 * 3; j++)
1088 outlimbs[j] |= inlimbs[j] & mask;
1092 /* get_bit returns the |i|th bit in |in| */
1093 static char get_bit(const felem_bytearray in, unsigned i)
1097 return (in[i >> 3] >> (i & 7)) & 1;
1101 * Interleaved point multiplication using precomputed point multiples: The
1102 * small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[], the scalars
1103 * in scalars[]. If g_scalar is non-NULL, we also add this multiple of the
1104 * generator, using certain (large) precomputed multiples in g_pre_comp.
1105 * Output point (X, Y, Z) is stored in x_out, y_out, z_out
1107 static void batch_mul(felem x_out, felem y_out, felem z_out,
1108 const felem_bytearray scalars[],
1109 const unsigned num_points, const u8 *g_scalar,
1110 const int mixed, const felem pre_comp[][17][3],
1111 const felem g_pre_comp[2][16][3])
1115 unsigned gen_mul = (g_scalar != NULL);
1116 felem nq[3], tmp[4];
1120 /* set nq to the point at infinity */
1121 memset(nq, 0, sizeof(nq));
1124 * Loop over all scalars msb-to-lsb, interleaving additions of multiples
1125 * of the generator (two in each of the last 28 rounds) and additions of
1126 * other points multiples (every 5th round).
1128 skip = 1; /* save two point operations in the first
1130 for (i = (num_points ? 220 : 27); i >= 0; --i) {
1133 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1135 /* add multiples of the generator */
1136 if (gen_mul && (i <= 27)) {
1137 /* first, look 28 bits upwards */
1138 bits = get_bit(g_scalar, i + 196) << 3;
1139 bits |= get_bit(g_scalar, i + 140) << 2;
1140 bits |= get_bit(g_scalar, i + 84) << 1;
1141 bits |= get_bit(g_scalar, i + 28);
1142 /* select the point to add, in constant time */
1143 select_point(bits, 16, g_pre_comp[1], tmp);
1146 /* value 1 below is argument for "mixed" */
1147 point_add(nq[0], nq[1], nq[2],
1148 nq[0], nq[1], nq[2], 1, tmp[0], tmp[1], tmp[2]);
1150 memcpy(nq, tmp, 3 * sizeof(felem));
1154 /* second, look at the current position */
1155 bits = get_bit(g_scalar, i + 168) << 3;
1156 bits |= get_bit(g_scalar, i + 112) << 2;
1157 bits |= get_bit(g_scalar, i + 56) << 1;
1158 bits |= get_bit(g_scalar, i);
1159 /* select the point to add, in constant time */
1160 select_point(bits, 16, g_pre_comp[0], tmp);
1161 point_add(nq[0], nq[1], nq[2],
1162 nq[0], nq[1], nq[2],
1163 1 /* mixed */ , tmp[0], tmp[1], tmp[2]);
1166 /* do other additions every 5 doublings */
1167 if (num_points && (i % 5 == 0)) {
1168 /* loop over all scalars */
1169 for (num = 0; num < num_points; ++num) {
1170 bits = get_bit(scalars[num], i + 4) << 5;
1171 bits |= get_bit(scalars[num], i + 3) << 4;
1172 bits |= get_bit(scalars[num], i + 2) << 3;
1173 bits |= get_bit(scalars[num], i + 1) << 2;
1174 bits |= get_bit(scalars[num], i) << 1;
1175 bits |= get_bit(scalars[num], i - 1);
1176 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1178 /* select the point to add or subtract */
1179 select_point(digit, 17, pre_comp[num], tmp);
1180 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative
1182 copy_conditional(tmp[1], tmp[3], sign);
1185 point_add(nq[0], nq[1], nq[2],
1186 nq[0], nq[1], nq[2],
1187 mixed, tmp[0], tmp[1], tmp[2]);
1189 memcpy(nq, tmp, 3 * sizeof(felem));
1195 felem_assign(x_out, nq[0]);
1196 felem_assign(y_out, nq[1]);
1197 felem_assign(z_out, nq[2]);
1200 /******************************************************************************/
1202 * FUNCTIONS TO MANAGE PRECOMPUTATION
1205 static NISTP224_PRE_COMP *nistp224_pre_comp_new()
1207 NISTP224_PRE_COMP *ret = OPENSSL_zalloc(sizeof(*ret));
1210 ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1214 ret->references = 1;
1216 ret->lock = CRYPTO_THREAD_lock_new();
1217 if (ret->lock == NULL) {
1218 ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1225 NISTP224_PRE_COMP *EC_nistp224_pre_comp_dup(NISTP224_PRE_COMP *p)
1229 CRYPTO_atomic_add(&p->references, 1, &i, p->lock);
1233 void EC_nistp224_pre_comp_free(NISTP224_PRE_COMP *p)
1240 CRYPTO_atomic_add(&p->references, -1, &i, p->lock);
1241 REF_PRINT_COUNT("EC_nistp224", x);
1244 REF_ASSERT_ISNT(i < 0);
1246 CRYPTO_THREAD_lock_free(p->lock);
1250 /******************************************************************************/
1252 * OPENSSL EC_METHOD FUNCTIONS
1255 int ec_GFp_nistp224_group_init(EC_GROUP *group)
1258 ret = ec_GFp_simple_group_init(group);
1259 group->a_is_minus3 = 1;
1263 int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1264 const BIGNUM *a, const BIGNUM *b,
1268 BN_CTX *new_ctx = NULL;
1269 BIGNUM *curve_p, *curve_a, *curve_b;
1272 if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1275 if (((curve_p = BN_CTX_get(ctx)) == NULL) ||
1276 ((curve_a = BN_CTX_get(ctx)) == NULL) ||
1277 ((curve_b = BN_CTX_get(ctx)) == NULL))
1279 BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
1280 BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
1281 BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
1282 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) || (BN_cmp(curve_b, b))) {
1283 ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
1284 EC_R_WRONG_CURVE_PARAMETERS);
1287 group->field_mod_func = BN_nist_mod_224;
1288 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1291 BN_CTX_free(new_ctx);
1296 * Takes the Jacobian coordinates (X, Y, Z) of a point and returns (X', Y') =
1299 int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
1300 const EC_POINT *point,
1301 BIGNUM *x, BIGNUM *y,
1304 felem z1, z2, x_in, y_in, x_out, y_out;
1307 if (EC_POINT_is_at_infinity(group, point)) {
1308 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1309 EC_R_POINT_AT_INFINITY);
1312 if ((!BN_to_felem(x_in, point->X)) || (!BN_to_felem(y_in, point->Y)) ||
1313 (!BN_to_felem(z1, point->Z)))
1316 felem_square(tmp, z2);
1317 felem_reduce(z1, tmp);
1318 felem_mul(tmp, x_in, z1);
1319 felem_reduce(x_in, tmp);
1320 felem_contract(x_out, x_in);
1322 if (!felem_to_BN(x, x_out)) {
1323 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1328 felem_mul(tmp, z1, z2);
1329 felem_reduce(z1, tmp);
1330 felem_mul(tmp, y_in, z1);
1331 felem_reduce(y_in, tmp);
1332 felem_contract(y_out, y_in);
1334 if (!felem_to_BN(y, y_out)) {
1335 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1343 static void make_points_affine(size_t num, felem points[ /* num */ ][3],
1344 felem tmp_felems[ /* num+1 */ ])
1347 * Runs in constant time, unless an input is the point at infinity (which
1348 * normally shouldn't happen).
1350 ec_GFp_nistp_points_make_affine_internal(num,
1354 (void (*)(void *))felem_one,
1356 (void (*)(void *, const void *))
1358 (void (*)(void *, const void *))
1359 felem_square_reduce, (void (*)
1366 (void (*)(void *, const void *))
1368 (void (*)(void *, const void *))
1373 * Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL
1374 * values Result is stored in r (r can equal one of the inputs).
1376 int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
1377 const BIGNUM *scalar, size_t num,
1378 const EC_POINT *points[],
1379 const BIGNUM *scalars[], BN_CTX *ctx)
1385 BN_CTX *new_ctx = NULL;
1386 BIGNUM *x, *y, *z, *tmp_scalar;
1387 felem_bytearray g_secret;
1388 felem_bytearray *secrets = NULL;
1389 felem (*pre_comp)[17][3] = NULL;
1390 felem *tmp_felems = NULL;
1392 int have_pre_comp = 0;
1393 size_t num_points = num;
1394 felem x_in, y_in, z_in, x_out, y_out, z_out;
1395 NISTP224_PRE_COMP *pre = NULL;
1396 const felem(*g_pre_comp)[16][3] = NULL;
1397 EC_POINT *generator = NULL;
1398 const EC_POINT *p = NULL;
1399 const BIGNUM *p_scalar = NULL;
1402 if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1405 if (((x = BN_CTX_get(ctx)) == NULL) ||
1406 ((y = BN_CTX_get(ctx)) == NULL) ||
1407 ((z = BN_CTX_get(ctx)) == NULL) ||
1408 ((tmp_scalar = BN_CTX_get(ctx)) == NULL))
1411 if (scalar != NULL) {
1412 pre = group->pre_comp.nistp224;
1414 /* we have precomputation, try to use it */
1415 g_pre_comp = (const felem(*)[16][3])pre->g_pre_comp;
1417 /* try to use the standard precomputation */
1418 g_pre_comp = &gmul[0];
1419 generator = EC_POINT_new(group);
1420 if (generator == NULL)
1422 /* get the generator from precomputation */
1423 if (!felem_to_BN(x, g_pre_comp[0][1][0]) ||
1424 !felem_to_BN(y, g_pre_comp[0][1][1]) ||
1425 !felem_to_BN(z, g_pre_comp[0][1][2])) {
1426 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1429 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1433 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1434 /* precomputation matches generator */
1438 * we don't have valid precomputation: treat the generator as a
1441 num_points = num_points + 1;
1444 if (num_points > 0) {
1445 if (num_points >= 3) {
1447 * unless we precompute multiples for just one or two points,
1448 * converting those into affine form is time well spent
1452 secrets = OPENSSL_zalloc(sizeof(*secrets) * num_points);
1453 pre_comp = OPENSSL_zalloc(sizeof(*pre_comp) * num_points);
1456 OPENSSL_malloc(sizeof(felem) * (num_points * 17 + 1));
1457 if ((secrets == NULL) || (pre_comp == NULL)
1458 || (mixed && (tmp_felems == NULL))) {
1459 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1464 * we treat NULL scalars as 0, and NULL points as points at infinity,
1465 * i.e., they contribute nothing to the linear combination
1467 for (i = 0; i < num_points; ++i) {
1470 p = EC_GROUP_get0_generator(group);
1473 /* the i^th point */
1475 p_scalar = scalars[i];
1477 if ((p_scalar != NULL) && (p != NULL)) {
1478 /* reduce scalar to 0 <= scalar < 2^224 */
1479 if ((BN_num_bits(p_scalar) > 224)
1480 || (BN_is_negative(p_scalar))) {
1482 * this is an unusual input, and we don't guarantee
1485 if (!BN_nnmod(tmp_scalar, p_scalar, group->order, ctx)) {
1486 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1489 num_bytes = BN_bn2lebinpad(tmp_scalar,
1490 secrets[i], sizeof(secrets[i]));
1492 num_bytes = BN_bn2lebinpad(p_scalar,
1493 secrets[i], sizeof(secrets[i]));
1495 if (num_bytes < 0) {
1496 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1499 /* precompute multiples */
1500 if ((!BN_to_felem(x_out, p->X)) ||
1501 (!BN_to_felem(y_out, p->Y)) ||
1502 (!BN_to_felem(z_out, p->Z)))
1504 felem_assign(pre_comp[i][1][0], x_out);
1505 felem_assign(pre_comp[i][1][1], y_out);
1506 felem_assign(pre_comp[i][1][2], z_out);
1507 for (j = 2; j <= 16; ++j) {
1509 point_add(pre_comp[i][j][0], pre_comp[i][j][1],
1510 pre_comp[i][j][2], pre_comp[i][1][0],
1511 pre_comp[i][1][1], pre_comp[i][1][2], 0,
1512 pre_comp[i][j - 1][0],
1513 pre_comp[i][j - 1][1],
1514 pre_comp[i][j - 1][2]);
1516 point_double(pre_comp[i][j][0], pre_comp[i][j][1],
1517 pre_comp[i][j][2], pre_comp[i][j / 2][0],
1518 pre_comp[i][j / 2][1],
1519 pre_comp[i][j / 2][2]);
1525 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
1528 /* the scalar for the generator */
1529 if ((scalar != NULL) && (have_pre_comp)) {
1530 memset(g_secret, 0, sizeof(g_secret));
1531 /* reduce scalar to 0 <= scalar < 2^224 */
1532 if ((BN_num_bits(scalar) > 224) || (BN_is_negative(scalar))) {
1534 * this is an unusual input, and we don't guarantee
1537 if (!BN_nnmod(tmp_scalar, scalar, group->order, ctx)) {
1538 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1541 num_bytes = BN_bn2lebinpad(tmp_scalar, g_secret, sizeof(g_secret));
1543 num_bytes = BN_bn2lebinpad(scalar, g_secret, sizeof(g_secret));
1545 /* do the multiplication with generator precomputation */
1546 batch_mul(x_out, y_out, z_out,
1547 (const felem_bytearray(*))secrets, num_points,
1549 mixed, (const felem(*)[17][3])pre_comp, g_pre_comp);
1551 /* do the multiplication without generator precomputation */
1552 batch_mul(x_out, y_out, z_out,
1553 (const felem_bytearray(*))secrets, num_points,
1554 NULL, mixed, (const felem(*)[17][3])pre_comp, NULL);
1556 /* reduce the output to its unique minimal representation */
1557 felem_contract(x_in, x_out);
1558 felem_contract(y_in, y_out);
1559 felem_contract(z_in, z_out);
1560 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
1561 (!felem_to_BN(z, z_in))) {
1562 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1565 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
1569 EC_POINT_free(generator);
1570 BN_CTX_free(new_ctx);
1571 OPENSSL_free(secrets);
1572 OPENSSL_free(pre_comp);
1573 OPENSSL_free(tmp_felems);
1577 int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
1580 NISTP224_PRE_COMP *pre = NULL;
1582 BN_CTX *new_ctx = NULL;
1584 EC_POINT *generator = NULL;
1585 felem tmp_felems[32];
1587 /* throw away old precomputation */
1588 EC_pre_comp_free(group);
1590 if ((ctx = new_ctx = BN_CTX_new()) == NULL)
1593 if (((x = BN_CTX_get(ctx)) == NULL) || ((y = BN_CTX_get(ctx)) == NULL))
1595 /* get the generator */
1596 if (group->generator == NULL)
1598 generator = EC_POINT_new(group);
1599 if (generator == NULL)
1601 BN_bin2bn(nistp224_curve_params[3], sizeof(felem_bytearray), x);
1602 BN_bin2bn(nistp224_curve_params[4], sizeof(felem_bytearray), y);
1603 if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
1605 if ((pre = nistp224_pre_comp_new()) == NULL)
1608 * if the generator is the standard one, use built-in precomputation
1610 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx)) {
1611 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
1614 if ((!BN_to_felem(pre->g_pre_comp[0][1][0], group->generator->X)) ||
1615 (!BN_to_felem(pre->g_pre_comp[0][1][1], group->generator->Y)) ||
1616 (!BN_to_felem(pre->g_pre_comp[0][1][2], group->generator->Z)))
1619 * compute 2^56*G, 2^112*G, 2^168*G for the first table, 2^28*G, 2^84*G,
1620 * 2^140*G, 2^196*G for the second one
1622 for (i = 1; i <= 8; i <<= 1) {
1623 point_double(pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1],
1624 pre->g_pre_comp[1][i][2], pre->g_pre_comp[0][i][0],
1625 pre->g_pre_comp[0][i][1], pre->g_pre_comp[0][i][2]);
1626 for (j = 0; j < 27; ++j) {
1627 point_double(pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1],
1628 pre->g_pre_comp[1][i][2], pre->g_pre_comp[1][i][0],
1629 pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1633 point_double(pre->g_pre_comp[0][2 * i][0],
1634 pre->g_pre_comp[0][2 * i][1],
1635 pre->g_pre_comp[0][2 * i][2], pre->g_pre_comp[1][i][0],
1636 pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1637 for (j = 0; j < 27; ++j) {
1638 point_double(pre->g_pre_comp[0][2 * i][0],
1639 pre->g_pre_comp[0][2 * i][1],
1640 pre->g_pre_comp[0][2 * i][2],
1641 pre->g_pre_comp[0][2 * i][0],
1642 pre->g_pre_comp[0][2 * i][1],
1643 pre->g_pre_comp[0][2 * i][2]);
1646 for (i = 0; i < 2; i++) {
1647 /* g_pre_comp[i][0] is the point at infinity */
1648 memset(pre->g_pre_comp[i][0], 0, sizeof(pre->g_pre_comp[i][0]));
1649 /* the remaining multiples */
1650 /* 2^56*G + 2^112*G resp. 2^84*G + 2^140*G */
1651 point_add(pre->g_pre_comp[i][6][0], pre->g_pre_comp[i][6][1],
1652 pre->g_pre_comp[i][6][2], pre->g_pre_comp[i][4][0],
1653 pre->g_pre_comp[i][4][1], pre->g_pre_comp[i][4][2],
1654 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1655 pre->g_pre_comp[i][2][2]);
1656 /* 2^56*G + 2^168*G resp. 2^84*G + 2^196*G */
1657 point_add(pre->g_pre_comp[i][10][0], pre->g_pre_comp[i][10][1],
1658 pre->g_pre_comp[i][10][2], pre->g_pre_comp[i][8][0],
1659 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1660 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1661 pre->g_pre_comp[i][2][2]);
1662 /* 2^112*G + 2^168*G resp. 2^140*G + 2^196*G */
1663 point_add(pre->g_pre_comp[i][12][0], pre->g_pre_comp[i][12][1],
1664 pre->g_pre_comp[i][12][2], pre->g_pre_comp[i][8][0],
1665 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1666 0, pre->g_pre_comp[i][4][0], pre->g_pre_comp[i][4][1],
1667 pre->g_pre_comp[i][4][2]);
1669 * 2^56*G + 2^112*G + 2^168*G resp. 2^84*G + 2^140*G + 2^196*G
1671 point_add(pre->g_pre_comp[i][14][0], pre->g_pre_comp[i][14][1],
1672 pre->g_pre_comp[i][14][2], pre->g_pre_comp[i][12][0],
1673 pre->g_pre_comp[i][12][1], pre->g_pre_comp[i][12][2],
1674 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1675 pre->g_pre_comp[i][2][2]);
1676 for (j = 1; j < 8; ++j) {
1677 /* odd multiples: add G resp. 2^28*G */
1678 point_add(pre->g_pre_comp[i][2 * j + 1][0],
1679 pre->g_pre_comp[i][2 * j + 1][1],
1680 pre->g_pre_comp[i][2 * j + 1][2],
1681 pre->g_pre_comp[i][2 * j][0],
1682 pre->g_pre_comp[i][2 * j][1],
1683 pre->g_pre_comp[i][2 * j][2], 0,
1684 pre->g_pre_comp[i][1][0], pre->g_pre_comp[i][1][1],
1685 pre->g_pre_comp[i][1][2]);
1688 make_points_affine(31, &(pre->g_pre_comp[0][1]), tmp_felems);
1691 SETPRECOMP(group, nistp224, pre);
1696 EC_POINT_free(generator);
1697 BN_CTX_free(new_ctx);
1698 EC_nistp224_pre_comp_free(pre);
1702 int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
1704 return HAVEPRECOMP(group, nistp224);