1 /* crypto/ec/ecp_nistp224.c */
3 * Written by Emilia Kasper (Google) for the OpenSSL project.
5 /* Copyright 2011 Google Inc.
7 * Licensed under the Apache License, Version 2.0 (the "License");
9 * you may not use this file except in compliance with the License.
10 * You may obtain a copy of the License at
12 * http://www.apache.org/licenses/LICENSE-2.0
14 * Unless required by applicable law or agreed to in writing, software
15 * distributed under the License is distributed on an "AS IS" BASIS,
16 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 * See the License for the specific language governing permissions and
18 * limitations under the License.
22 * A 64-bit implementation of the NIST P-224 elliptic curve point multiplication
24 * Inspired by Daniel J. Bernstein's public domain nistp224 implementation
25 * and Adam Langley's public domain 64-bit C implementation of curve25519
28 #include <openssl/opensslconf.h>
29 #ifndef OPENSSL_NO_EC_NISTP_64_GCC_128
33 #include <openssl/err.h>
36 #if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 1))
37 /* even with gcc, the typedef won't work for 32-bit platforms */
38 typedef __uint128_t uint128_t; /* nonstandard; implemented by gcc on 64-bit platforms */
40 #error "Need GCC 3.1 or later to define type uint128_t"
48 /******************************************************************************/
50 * INTERNAL REPRESENTATION OF FIELD ELEMENTS
52 * Field elements are represented as a_0 + 2^56*a_1 + 2^112*a_2 + 2^168*a_3
53 * using 64-bit coefficients called 'limbs',
54 * and sometimes (for multiplication results) as
55 * 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
56 * using 128-bit coefficients called 'widelimbs'.
57 * A 4-limb representation is an 'felem';
58 * a 7-widelimb representation is a 'widefelem'.
59 * Even within felems, bits of adjacent limbs overlap, and we don't always
60 * reduce the representations: we ensure that inputs to each felem
61 * multiplication satisfy a_i < 2^60, so outputs satisfy b_i < 4*2^60*2^60,
62 * and fit into a 128-bit word without overflow. The coefficients are then
63 * again partially reduced to obtain an felem satisfying a_i < 2^57.
64 * We only reduce to the unique minimal representation at the end of the
68 typedef uint64_t limb;
69 typedef uint128_t widelimb;
71 typedef limb felem[4];
72 typedef widelimb widefelem[7];
74 /* Field element represented as a byte arrary.
75 * 28*8 = 224 bits is also the group order size for the elliptic curve,
76 * and we also use this type for scalars for point multiplication.
78 typedef u8 felem_bytearray[28];
80 static const felem_bytearray nistp224_curve_params[5] = {
81 {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF, /* p */
82 0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0x00,0x00,0x00,0x00,
83 0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x01},
84 {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF, /* a */
85 0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,0xFF,0xFF,0xFF,0xFF,
86 0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFE},
87 {0xB4,0x05,0x0A,0x85,0x0C,0x04,0xB3,0xAB,0xF5,0x41, /* b */
88 0x32,0x56,0x50,0x44,0xB0,0xB7,0xD7,0xBF,0xD8,0xBA,
89 0x27,0x0B,0x39,0x43,0x23,0x55,0xFF,0xB4},
90 {0xB7,0x0E,0x0C,0xBD,0x6B,0xB4,0xBF,0x7F,0x32,0x13, /* x */
91 0x90,0xB9,0x4A,0x03,0xC1,0xD3,0x56,0xC2,0x11,0x22,
92 0x34,0x32,0x80,0xD6,0x11,0x5C,0x1D,0x21},
93 {0xbd,0x37,0x63,0x88,0xb5,0xf7,0x23,0xfb,0x4c,0x22, /* y */
94 0xdf,0xe6,0xcd,0x43,0x75,0xa0,0x5a,0x07,0x47,0x64,
95 0x44,0xd5,0x81,0x99,0x85,0x00,0x7e,0x34}
99 * Precomputed multiples of the standard generator
100 * Points are given in coordinates (X, Y, Z) where Z normally is 1
101 * (0 for the point at infinity).
102 * For each field element, slice a_0 is word 0, etc.
104 * The table has 2 * 16 elements, starting with the following:
105 * index | bits | point
106 * ------+---------+------------------------------
109 * 2 | 0 0 1 0 | 2^56G
110 * 3 | 0 0 1 1 | (2^56 + 1)G
111 * 4 | 0 1 0 0 | 2^112G
112 * 5 | 0 1 0 1 | (2^112 + 1)G
113 * 6 | 0 1 1 0 | (2^112 + 2^56)G
114 * 7 | 0 1 1 1 | (2^112 + 2^56 + 1)G
115 * 8 | 1 0 0 0 | 2^168G
116 * 9 | 1 0 0 1 | (2^168 + 1)G
117 * 10 | 1 0 1 0 | (2^168 + 2^56)G
118 * 11 | 1 0 1 1 | (2^168 + 2^56 + 1)G
119 * 12 | 1 1 0 0 | (2^168 + 2^112)G
120 * 13 | 1 1 0 1 | (2^168 + 2^112 + 1)G
121 * 14 | 1 1 1 0 | (2^168 + 2^112 + 2^56)G
122 * 15 | 1 1 1 1 | (2^168 + 2^112 + 2^56 + 1)G
123 * followed by a copy of this with each element multiplied by 2^28.
125 * The reason for this is so that we can clock bits into four different
126 * locations when doing simple scalar multiplies against the base point,
127 * and then another four locations using the second 16 elements.
129 static const felem gmul[2][16][3] =
133 {{0x3280d6115c1d21, 0xc1d356c2112234, 0x7f321390b94a03, 0xb70e0cbd6bb4bf},
134 {0xd5819985007e34, 0x75a05a07476444, 0xfb4c22dfe6cd43, 0xbd376388b5f723},
136 {{0xfd9675666ebbe9, 0xbca7664d40ce5e, 0x2242df8d8a2a43, 0x1f49bbb0f99bc5},
137 {0x29e0b892dc9c43, 0xece8608436e662, 0xdc858f185310d0, 0x9812dd4eb8d321},
139 {{0x6d3e678d5d8eb8, 0x559eed1cb362f1, 0x16e9a3bbce8a3f, 0xeedcccd8c2a748},
140 {0xf19f90ed50266d, 0xabf2b4bf65f9df, 0x313865468fafec, 0x5cb379ba910a17},
142 {{0x0641966cab26e3, 0x91fb2991fab0a0, 0xefec27a4e13a0b, 0x0499aa8a5f8ebe},
143 {0x7510407766af5d, 0x84d929610d5450, 0x81d77aae82f706, 0x6916f6d4338c5b},
145 {{0xea95ac3b1f15c6, 0x086000905e82d4, 0xdd323ae4d1c8b1, 0x932b56be7685a3},
146 {0x9ef93dea25dbbf, 0x41665960f390f0, 0xfdec76dbe2a8a7, 0x523e80f019062a},
148 {{0x822fdd26732c73, 0xa01c83531b5d0f, 0x363f37347c1ba4, 0xc391b45c84725c},
149 {0xbbd5e1b2d6ad24, 0xddfbcde19dfaec, 0xc393da7e222a7f, 0x1efb7890ede244},
151 {{0x4c9e90ca217da1, 0xd11beca79159bb, 0xff8d33c2c98b7c, 0x2610b39409f849},
152 {0x44d1352ac64da0, 0xcdbb7b2c46b4fb, 0x966c079b753c89, 0xfe67e4e820b112},
154 {{0xe28cae2df5312d, 0xc71b61d16f5c6e, 0x79b7619a3e7c4c, 0x05c73240899b47},
155 {0x9f7f6382c73e3a, 0x18615165c56bda, 0x641fab2116fd56, 0x72855882b08394},
157 {{0x0469182f161c09, 0x74a98ca8d00fb5, 0xb89da93489a3e0, 0x41c98768fb0c1d},
158 {0xe5ea05fb32da81, 0x3dce9ffbca6855, 0x1cfe2d3fbf59e6, 0x0e5e03408738a7},
160 {{0xdab22b2333e87f, 0x4430137a5dd2f6, 0xe03ab9f738beb8, 0xcb0c5d0dc34f24},
161 {0x764a7df0c8fda5, 0x185ba5c3fa2044, 0x9281d688bcbe50, 0xc40331df893881},
163 {{0xb89530796f0f60, 0xade92bd26909a3, 0x1a0c83fb4884da, 0x1765bf22a5a984},
164 {0x772a9ee75db09e, 0x23bc6c67cec16f, 0x4c1edba8b14e2f, 0xe2a215d9611369},
166 {{0x571e509fb5efb3, 0xade88696410552, 0xc8ae85fada74fe, 0x6c7e4be83bbde3},
167 {0xff9f51160f4652, 0xb47ce2495a6539, 0xa2946c53b582f4, 0x286d2db3ee9a60},
169 {{0x40bbd5081a44af, 0x0995183b13926c, 0xbcefba6f47f6d0, 0x215619e9cc0057},
170 {0x8bc94d3b0df45e, 0xf11c54a3694f6f, 0x8631b93cdfe8b5, 0xe7e3f4b0982db9},
172 {{0xb17048ab3e1c7b, 0xac38f36ff8a1d8, 0x1c29819435d2c6, 0xc813132f4c07e9},
173 {0x2891425503b11f, 0x08781030579fea, 0xf5426ba5cc9674, 0x1e28ebf18562bc},
175 {{0x9f31997cc864eb, 0x06cd91d28b5e4c, 0xff17036691a973, 0xf1aef351497c58},
176 {0xdd1f2d600564ff, 0xdead073b1402db, 0x74a684435bd693, 0xeea7471f962558},
181 {{0x9665266dddf554, 0x9613d78b60ef2d, 0xce27a34cdba417, 0xd35ab74d6afc31},
182 {0x85ccdd22deb15e, 0x2137e5783a6aab, 0xa141cffd8c93c6, 0x355a1830e90f2d},
184 {{0x1a494eadaade65, 0xd6da4da77fe53c, 0xe7992996abec86, 0x65c3553c6090e3},
185 {0xfa610b1fb09346, 0xf1c6540b8a4aaf, 0xc51a13ccd3cbab, 0x02995b1b18c28a},
187 {{0x7874568e7295ef, 0x86b419fbe38d04, 0xdc0690a7550d9a, 0xd3966a44beac33},
188 {0x2b7280ec29132f, 0xbeaa3b6a032df3, 0xdc7dd88ae41200, 0xd25e2513e3a100},
190 {{0x924857eb2efafd, 0xac2bce41223190, 0x8edaa1445553fc, 0x825800fd3562d5},
191 {0x8d79148ea96621, 0x23a01c3dd9ed8d, 0xaf8b219f9416b5, 0xd8db0cc277daea},
193 {{0x76a9c3b1a700f0, 0xe9acd29bc7e691, 0x69212d1a6b0327, 0x6322e97fe154be},
194 {0x469fc5465d62aa, 0x8d41ed18883b05, 0x1f8eae66c52b88, 0xe4fcbe9325be51},
196 {{0x825fdf583cac16, 0x020b857c7b023a, 0x683c17744b0165, 0x14ffd0a2daf2f1},
197 {0x323b36184218f9, 0x4944ec4e3b47d4, 0xc15b3080841acf, 0x0bced4b01a28bb},
199 {{0x92ac22230df5c4, 0x52f33b4063eda8, 0xcb3f19870c0c93, 0x40064f2ba65233},
200 {0xfe16f0924f8992, 0x012da25af5b517, 0x1a57bb24f723a6, 0x06f8bc76760def},
202 {{0x4a7084f7817cb9, 0xbcab0738ee9a78, 0x3ec11e11d9c326, 0xdc0fe90e0f1aae},
203 {0xcf639ea5f98390, 0x5c350aa22ffb74, 0x9afae98a4047b7, 0x956ec2d617fc45},
205 {{0x4306d648c1be6a, 0x9247cd8bc9a462, 0xf5595e377d2f2e, 0xbd1c3caff1a52e},
206 {0x045e14472409d0, 0x29f3e17078f773, 0x745a602b2d4f7d, 0x191837685cdfbb},
208 {{0x5b6ee254a8cb79, 0x4953433f5e7026, 0xe21faeb1d1def4, 0xc4c225785c09de},
209 {0x307ce7bba1e518, 0x31b125b1036db8, 0x47e91868839e8f, 0xc765866e33b9f3},
211 {{0x3bfece24f96906, 0x4794da641e5093, 0xde5df64f95db26, 0x297ecd89714b05},
212 {0x701bd3ebb2c3aa, 0x7073b4f53cb1d5, 0x13c5665658af16, 0x9895089d66fe58},
214 {{0x0fef05f78c4790, 0x2d773633b05d2e, 0x94229c3a951c94, 0xbbbd70df4911bb},
215 {0xb2c6963d2c1168, 0x105f47a72b0d73, 0x9fdf6111614080, 0x7b7e94b39e67b0},
217 {{0xad1a7d6efbe2b3, 0xf012482c0da69d, 0x6b3bdf12438345, 0x40d7558d7aa4d9},
218 {0x8a09fffb5c6d3d, 0x9a356e5d9ffd38, 0x5973f15f4f9b1c, 0xdcd5f59f63c3ea},
220 {{0xacf39f4c5ca7ab, 0x4c8071cc5fd737, 0xc64e3602cd1184, 0x0acd4644c9abba},
221 {0x6c011a36d8bf6e, 0xfecd87ba24e32a, 0x19f6f56574fad8, 0x050b204ced9405},
223 {{0xed4f1cae7d9a96, 0x5ceef7ad94c40a, 0x778e4a3bf3ef9b, 0x7405783dc3b55e},
224 {0x32477c61b6e8c6, 0xb46a97570f018b, 0x91176d0a7e95d1, 0x3df90fbc4c7d0e},
227 /* Precomputation for the group generator. */
229 felem g_pre_comp[2][16][3];
233 const EC_METHOD *EC_GFp_nistp224_method(void)
235 static const EC_METHOD ret = {
236 EC_FLAGS_DEFAULT_OCT,
237 NID_X9_62_prime_field,
238 ec_GFp_nistp224_group_init,
239 ec_GFp_simple_group_finish,
240 ec_GFp_simple_group_clear_finish,
241 ec_GFp_nist_group_copy,
242 ec_GFp_nistp224_group_set_curve,
243 ec_GFp_simple_group_get_curve,
244 ec_GFp_simple_group_get_degree,
245 ec_GFp_simple_group_check_discriminant,
246 ec_GFp_simple_point_init,
247 ec_GFp_simple_point_finish,
248 ec_GFp_simple_point_clear_finish,
249 ec_GFp_simple_point_copy,
250 ec_GFp_simple_point_set_to_infinity,
251 ec_GFp_simple_set_Jprojective_coordinates_GFp,
252 ec_GFp_simple_get_Jprojective_coordinates_GFp,
253 ec_GFp_simple_point_set_affine_coordinates,
254 ec_GFp_nistp224_point_get_affine_coordinates,
255 0 /* point_set_compressed_coordinates */,
260 ec_GFp_simple_invert,
261 ec_GFp_simple_is_at_infinity,
262 ec_GFp_simple_is_on_curve,
264 ec_GFp_simple_make_affine,
265 ec_GFp_simple_points_make_affine,
266 ec_GFp_nistp224_points_mul,
267 ec_GFp_nistp224_precompute_mult,
268 ec_GFp_nistp224_have_precompute_mult,
269 ec_GFp_nist_field_mul,
270 ec_GFp_nist_field_sqr,
272 0 /* field_encode */,
273 0 /* field_decode */,
274 0 /* field_set_to_one */ };
279 /* Helper functions to convert field elements to/from internal representation */
280 static void bin28_to_felem(felem out, const u8 in[28])
282 out[0] = *((const uint64_t *)(in)) & 0x00ffffffffffffff;
283 out[1] = (*((const uint64_t *)(in+7))) & 0x00ffffffffffffff;
284 out[2] = (*((const uint64_t *)(in+14))) & 0x00ffffffffffffff;
285 out[3] = (*((const uint64_t *)(in+21))) & 0x00ffffffffffffff;
288 static void felem_to_bin28(u8 out[28], const felem in)
291 for (i = 0; i < 7; ++i)
293 out[i] = in[0]>>(8*i);
294 out[i+7] = in[1]>>(8*i);
295 out[i+14] = in[2]>>(8*i);
296 out[i+21] = in[3]>>(8*i);
300 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */
301 static void flip_endian(u8 *out, const u8 *in, unsigned len)
304 for (i = 0; i < len; ++i)
305 out[i] = in[len-1-i];
308 /* From OpenSSL BIGNUM to internal representation */
309 static int BN_to_felem(felem out, const BIGNUM *bn)
311 felem_bytearray b_in;
312 felem_bytearray b_out;
315 /* BN_bn2bin eats leading zeroes */
316 memset(b_out, 0, sizeof b_out);
317 num_bytes = BN_num_bytes(bn);
318 if (num_bytes > sizeof b_out)
320 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
323 if (BN_is_negative(bn))
325 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
328 num_bytes = BN_bn2bin(bn, b_in);
329 flip_endian(b_out, b_in, num_bytes);
330 bin28_to_felem(out, b_out);
334 /* From internal representation to OpenSSL BIGNUM */
335 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
337 felem_bytearray b_in, b_out;
338 felem_to_bin28(b_in, in);
339 flip_endian(b_out, b_in, sizeof b_out);
340 return BN_bin2bn(b_out, sizeof b_out, out);
343 /******************************************************************************/
347 * Field operations, using the internal representation of field elements.
348 * NB! These operations are specific to our point multiplication and cannot be
349 * expected to be correct in general - e.g., multiplication with a large scalar
350 * will cause an overflow.
354 static void felem_one(felem out)
362 static void felem_assign(felem out, const felem in)
370 /* Sum two field elements: out += in */
371 static void felem_sum(felem out, const felem in)
379 /* Get negative value: out = -in */
380 /* Assumes in[i] < 2^57 */
381 static void felem_neg(felem out, const felem in)
383 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
384 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
385 static const limb two58m42m2 = (((limb) 1) << 58) -
386 (((limb) 1) << 42) - (((limb) 1) << 2);
388 /* Set to 0 mod 2^224-2^96+1 to ensure out > in */
389 out[0] = two58p2 - in[0];
390 out[1] = two58m42m2 - in[1];
391 out[2] = two58m2 - in[2];
392 out[3] = two58m2 - in[3];
395 /* Subtract field elements: out -= in */
396 /* Assumes in[i] < 2^57 */
397 static void felem_diff(felem out, const felem in)
399 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
400 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
401 static const limb two58m42m2 = (((limb) 1) << 58) -
402 (((limb) 1) << 42) - (((limb) 1) << 2);
404 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
406 out[1] += two58m42m2;
416 /* Subtract in unreduced 128-bit mode: out -= in */
417 /* Assumes in[i] < 2^119 */
418 static void widefelem_diff(widefelem out, const widefelem in)
420 static const widelimb two120 = ((widelimb) 1) << 120;
421 static const widelimb two120m64 = (((widelimb) 1) << 120) -
422 (((widelimb) 1) << 64);
423 static const widelimb two120m104m64 = (((widelimb) 1) << 120) -
424 (((widelimb) 1) << 104) - (((widelimb) 1) << 64);
426 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
431 out[4] += two120m104m64;
444 /* Subtract in mixed mode: out128 -= in64 */
446 static void felem_diff_128_64(widefelem out, const felem in)
448 static const widelimb two64p8 = (((widelimb) 1) << 64) +
449 (((widelimb) 1) << 8);
450 static const widelimb two64m8 = (((widelimb) 1) << 64) -
451 (((widelimb) 1) << 8);
452 static const widelimb two64m48m8 = (((widelimb) 1) << 64) -
453 (((widelimb) 1) << 48) - (((widelimb) 1) << 8);
455 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
457 out[1] += two64m48m8;
467 /* Multiply a field element by a scalar: out = out * scalar
468 * The scalars we actually use are small, so results fit without overflow */
469 static void felem_scalar(felem out, const limb scalar)
477 /* Multiply an unreduced field element by a scalar: out = out * scalar
478 * The scalars we actually use are small, so results fit without overflow */
479 static void widefelem_scalar(widefelem out, const widelimb scalar)
490 /* Square a field element: out = in^2 */
491 static void felem_square(widefelem out, const felem in)
493 limb tmp0, tmp1, tmp2;
494 tmp0 = 2 * in[0]; tmp1 = 2 * in[1]; tmp2 = 2 * in[2];
495 out[0] = ((widelimb) in[0]) * in[0];
496 out[1] = ((widelimb) in[0]) * tmp1;
497 out[2] = ((widelimb) in[0]) * tmp2 + ((widelimb) in[1]) * in[1];
498 out[3] = ((widelimb) in[3]) * tmp0 +
499 ((widelimb) in[1]) * tmp2;
500 out[4] = ((widelimb) in[3]) * tmp1 + ((widelimb) in[2]) * in[2];
501 out[5] = ((widelimb) in[3]) * tmp2;
502 out[6] = ((widelimb) in[3]) * in[3];
505 /* Multiply two field elements: out = in1 * in2 */
506 static void felem_mul(widefelem out, const felem in1, const felem in2)
508 out[0] = ((widelimb) in1[0]) * in2[0];
509 out[1] = ((widelimb) in1[0]) * in2[1] + ((widelimb) in1[1]) * in2[0];
510 out[2] = ((widelimb) in1[0]) * in2[2] + ((widelimb) in1[1]) * in2[1] +
511 ((widelimb) in1[2]) * in2[0];
512 out[3] = ((widelimb) in1[0]) * in2[3] + ((widelimb) in1[1]) * in2[2] +
513 ((widelimb) in1[2]) * in2[1] + ((widelimb) in1[3]) * in2[0];
514 out[4] = ((widelimb) in1[1]) * in2[3] + ((widelimb) in1[2]) * in2[2] +
515 ((widelimb) in1[3]) * in2[1];
516 out[5] = ((widelimb) in1[2]) * in2[3] + ((widelimb) in1[3]) * in2[2];
517 out[6] = ((widelimb) in1[3]) * in2[3];
521 * Reduce seven 128-bit coefficients to four 64-bit coefficients.
522 * Requires in[i] < 2^126,
523 * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] <= 2^56 + 2^16 */
524 static void felem_reduce(felem out, const widefelem in)
526 static const widelimb two127p15 = (((widelimb) 1) << 127) +
527 (((widelimb) 1) << 15);
528 static const widelimb two127m71 = (((widelimb) 1) << 127) -
529 (((widelimb) 1) << 71);
530 static const widelimb two127m71m55 = (((widelimb) 1) << 127) -
531 (((widelimb) 1) << 71) - (((widelimb) 1) << 55);
534 /* Add 0 mod 2^224-2^96+1 to ensure all differences are positive */
535 output[0] = in[0] + two127p15;
536 output[1] = in[1] + two127m71m55;
537 output[2] = in[2] + two127m71;
541 /* Eliminate in[4], in[5], in[6] */
542 output[4] += in[6] >> 16;
543 output[3] += (in[6] & 0xffff) << 40;
546 output[3] += in[5] >> 16;
547 output[2] += (in[5] & 0xffff) << 40;
550 output[2] += output[4] >> 16;
551 output[1] += (output[4] & 0xffff) << 40;
552 output[0] -= output[4];
554 /* Carry 2 -> 3 -> 4 */
555 output[3] += output[2] >> 56;
556 output[2] &= 0x00ffffffffffffff;
558 output[4] = output[3] >> 56;
559 output[3] &= 0x00ffffffffffffff;
561 /* Now output[2] < 2^56, output[3] < 2^56, output[4] < 2^72 */
563 /* Eliminate output[4] */
564 output[2] += output[4] >> 16;
565 /* output[2] < 2^56 + 2^56 = 2^57 */
566 output[1] += (output[4] & 0xffff) << 40;
567 output[0] -= output[4];
569 /* Carry 0 -> 1 -> 2 -> 3 */
570 output[1] += output[0] >> 56;
571 out[0] = output[0] & 0x00ffffffffffffff;
573 output[2] += output[1] >> 56;
574 /* output[2] < 2^57 + 2^72 */
575 out[1] = output[1] & 0x00ffffffffffffff;
576 output[3] += output[2] >> 56;
577 /* output[3] <= 2^56 + 2^16 */
578 out[2] = output[2] & 0x00ffffffffffffff;
581 * out[0] < 2^56, out[1] < 2^56, out[2] < 2^56,
582 * out[3] <= 2^56 + 2^16 (due to final carry),
588 static void felem_square_reduce(felem out, const felem in)
591 felem_square(tmp, in);
592 felem_reduce(out, tmp);
595 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
598 felem_mul(tmp, in1, in2);
599 felem_reduce(out, tmp);
602 /* Reduce to unique minimal representation.
603 * Requires 0 <= in < 2*p (always call felem_reduce first) */
604 static void felem_contract(felem out, const felem in)
606 static const int64_t two56 = ((limb) 1) << 56;
607 /* 0 <= in < 2*p, p = 2^224 - 2^96 + 1 */
608 /* if in > p , reduce in = in - 2^224 + 2^96 - 1 */
614 /* Case 1: a = 1 iff in >= 2^224 */
618 tmp[3] &= 0x00ffffffffffffff;
619 /* Case 2: a = 0 iff p <= in < 2^224, i.e.,
620 * the high 128 bits are all 1 and the lower part is non-zero */
621 a = ((in[3] & in[2] & (in[1] | 0x000000ffffffffff)) + 1) |
622 (((int64_t)(in[0] + (in[1] & 0x000000ffffffffff)) - 1) >> 63);
623 a &= 0x00ffffffffffffff;
624 /* turn a into an all-one mask (if a = 0) or an all-zero mask */
626 /* subtract 2^224 - 2^96 + 1 if a is all-one*/
627 tmp[3] &= a ^ 0xffffffffffffffff;
628 tmp[2] &= a ^ 0xffffffffffffffff;
629 tmp[1] &= (a ^ 0xffffffffffffffff) | 0x000000ffffffffff;
632 /* eliminate negative coefficients: if tmp[0] is negative, tmp[1] must
633 * be non-zero, so we only need one step */
638 /* carry 1 -> 2 -> 3 */
639 tmp[2] += tmp[1] >> 56;
640 tmp[1] &= 0x00ffffffffffffff;
642 tmp[3] += tmp[2] >> 56;
643 tmp[2] &= 0x00ffffffffffffff;
645 /* Now 0 <= out < p */
652 /* Zero-check: returns 1 if input is 0, and 0 otherwise.
653 * We know that field elements are reduced to in < 2^225,
654 * so we only need to check three cases: 0, 2^224 - 2^96 + 1,
655 * and 2^225 - 2^97 + 2 */
656 static limb felem_is_zero(const felem in)
658 limb zero, two224m96p1, two225m97p2;
660 zero = in[0] | in[1] | in[2] | in[3];
661 zero = (((int64_t)(zero) - 1) >> 63) & 1;
662 two224m96p1 = (in[0] ^ 1) | (in[1] ^ 0x00ffff0000000000)
663 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x00ffffffffffffff);
664 two224m96p1 = (((int64_t)(two224m96p1) - 1) >> 63) & 1;
665 two225m97p2 = (in[0] ^ 2) | (in[1] ^ 0x00fffe0000000000)
666 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x01ffffffffffffff);
667 two225m97p2 = (((int64_t)(two225m97p2) - 1) >> 63) & 1;
668 return (zero | two224m96p1 | two225m97p2);
671 static limb felem_is_zero_int(const felem in)
673 return (int) (felem_is_zero(in) & ((limb)1));
676 /* Invert a field element */
677 /* Computation chain copied from djb's code */
678 static void felem_inv(felem out, const felem in)
680 felem ftmp, ftmp2, ftmp3, ftmp4;
684 felem_square(tmp, in); felem_reduce(ftmp, tmp); /* 2 */
685 felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp); /* 2^2 - 1 */
686 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp); /* 2^3 - 2 */
687 felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp); /* 2^3 - 1 */
688 felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp); /* 2^4 - 2 */
689 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp); /* 2^5 - 4 */
690 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp); /* 2^6 - 8 */
691 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp); /* 2^6 - 1 */
692 felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp); /* 2^7 - 2 */
693 for (i = 0; i < 5; ++i) /* 2^12 - 2^6 */
695 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
697 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp2, tmp); /* 2^12 - 1 */
698 felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp); /* 2^13 - 2 */
699 for (i = 0; i < 11; ++i) /* 2^24 - 2^12 */
701 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
703 felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp2, tmp); /* 2^24 - 1 */
704 felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp); /* 2^25 - 2 */
705 for (i = 0; i < 23; ++i) /* 2^48 - 2^24 */
707 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
709 felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp3, tmp); /* 2^48 - 1 */
710 felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp); /* 2^49 - 2 */
711 for (i = 0; i < 47; ++i) /* 2^96 - 2^48 */
713 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
715 felem_mul(tmp, ftmp3, ftmp4); felem_reduce(ftmp3, tmp); /* 2^96 - 1 */
716 felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp); /* 2^97 - 2 */
717 for (i = 0; i < 23; ++i) /* 2^120 - 2^24 */
719 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
721 felem_mul(tmp, ftmp2, ftmp4); felem_reduce(ftmp2, tmp); /* 2^120 - 1 */
722 for (i = 0; i < 6; ++i) /* 2^126 - 2^6 */
724 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
726 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp); /* 2^126 - 1 */
727 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp); /* 2^127 - 2 */
728 felem_mul(tmp, ftmp, in); felem_reduce(ftmp, tmp); /* 2^127 - 1 */
729 for (i = 0; i < 97; ++i) /* 2^224 - 2^97 */
731 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp);
733 felem_mul(tmp, ftmp, ftmp3); felem_reduce(out, tmp); /* 2^224 - 2^96 - 1 */
736 /* Copy in constant time:
737 * if icopy == 1, copy in to out,
738 * if icopy == 0, copy out to itself. */
740 copy_conditional(felem out, const felem in, limb icopy)
743 /* icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one */
744 const limb copy = -icopy;
745 for (i = 0; i < 4; ++i)
747 const limb tmp = copy & (in[i] ^ out[i]);
752 /******************************************************************************/
754 * ELLIPTIC CURVE POINT OPERATIONS
756 * Points are represented in Jacobian projective coordinates:
757 * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
758 * or to the point at infinity if Z == 0.
763 * Double an elliptic curve point:
764 * (X', Y', Z') = 2 * (X, Y, Z), where
765 * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
766 * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^2
767 * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
768 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
769 * while x_out == y_in is not (maybe this works, but it's not tested).
772 point_double(felem x_out, felem y_out, felem z_out,
773 const felem x_in, const felem y_in, const felem z_in)
776 felem delta, gamma, beta, alpha, ftmp, ftmp2;
778 felem_assign(ftmp, x_in);
779 felem_assign(ftmp2, x_in);
782 felem_square(tmp, z_in);
783 felem_reduce(delta, tmp);
786 felem_square(tmp, y_in);
787 felem_reduce(gamma, tmp);
790 felem_mul(tmp, x_in, gamma);
791 felem_reduce(beta, tmp);
793 /* alpha = 3*(x-delta)*(x+delta) */
794 felem_diff(ftmp, delta);
795 /* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
796 felem_sum(ftmp2, delta);
797 /* ftmp2[i] < 2^57 + 2^57 = 2^58 */
798 felem_scalar(ftmp2, 3);
799 /* ftmp2[i] < 3 * 2^58 < 2^60 */
800 felem_mul(tmp, ftmp, ftmp2);
801 /* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
802 felem_reduce(alpha, tmp);
804 /* x' = alpha^2 - 8*beta */
805 felem_square(tmp, alpha);
806 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
807 felem_assign(ftmp, beta);
808 felem_scalar(ftmp, 8);
809 /* ftmp[i] < 8 * 2^57 = 2^60 */
810 felem_diff_128_64(tmp, ftmp);
811 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
812 felem_reduce(x_out, tmp);
814 /* z' = (y + z)^2 - gamma - delta */
815 felem_sum(delta, gamma);
816 /* delta[i] < 2^57 + 2^57 = 2^58 */
817 felem_assign(ftmp, y_in);
818 felem_sum(ftmp, z_in);
819 /* ftmp[i] < 2^57 + 2^57 = 2^58 */
820 felem_square(tmp, ftmp);
821 /* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
822 felem_diff_128_64(tmp, delta);
823 /* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
824 felem_reduce(z_out, tmp);
826 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
827 felem_scalar(beta, 4);
828 /* beta[i] < 4 * 2^57 = 2^59 */
829 felem_diff(beta, x_out);
830 /* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
831 felem_mul(tmp, alpha, beta);
832 /* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
833 felem_square(tmp2, gamma);
834 /* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
835 widefelem_scalar(tmp2, 8);
836 /* tmp2[i] < 8 * 2^116 = 2^119 */
837 widefelem_diff(tmp, tmp2);
838 /* tmp[i] < 2^119 + 2^120 < 2^121 */
839 felem_reduce(y_out, tmp);
843 * Add two elliptic curve points:
844 * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
845 * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
846 * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
847 * 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) -
848 * Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
849 * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2)
851 * This runs faster if 'mixed' is set, which requires Z_2 = 1 or Z_2 = 0.
854 /* This function is not entirely constant-time:
855 * it includes a branch for checking whether the two input points are equal,
856 * (while not equal to the point at infinity).
857 * This case never happens during single point multiplication,
858 * so there is no timing leak for ECDH or ECDSA signing. */
859 static void point_add(felem x3, felem y3, felem z3,
860 const felem x1, const felem y1, const felem z1,
861 const int mixed, const felem x2, const felem y2, const felem z2)
863 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, x_out, y_out, z_out;
865 limb z1_is_zero, z2_is_zero, x_equal, y_equal;
870 felem_square(tmp, z2);
871 felem_reduce(ftmp2, tmp);
874 felem_mul(tmp, ftmp2, z2);
875 felem_reduce(ftmp4, tmp);
877 /* ftmp4 = z2^3*y1 */
878 felem_mul(tmp2, ftmp4, y1);
879 felem_reduce(ftmp4, tmp2);
881 /* ftmp2 = z2^2*x1 */
882 felem_mul(tmp2, ftmp2, x1);
883 felem_reduce(ftmp2, tmp2);
887 /* We'll assume z2 = 1 (special case z2 = 0 is handled later) */
889 /* ftmp4 = z2^3*y1 */
890 felem_assign(ftmp4, y1);
892 /* ftmp2 = z2^2*x1 */
893 felem_assign(ftmp2, x1);
897 felem_square(tmp, z1);
898 felem_reduce(ftmp, tmp);
901 felem_mul(tmp, ftmp, z1);
902 felem_reduce(ftmp3, tmp);
905 felem_mul(tmp, ftmp3, y2);
906 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
908 /* ftmp3 = z1^3*y2 - z2^3*y1 */
909 felem_diff_128_64(tmp, ftmp4);
910 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
911 felem_reduce(ftmp3, tmp);
914 felem_mul(tmp, ftmp, x2);
915 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
917 /* ftmp = z1^2*x2 - z2^2*x1 */
918 felem_diff_128_64(tmp, ftmp2);
919 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
920 felem_reduce(ftmp, tmp);
922 /* the formulae are incorrect if the points are equal
923 * so we check for this and do doubling if this happens */
924 x_equal = felem_is_zero(ftmp);
925 y_equal = felem_is_zero(ftmp3);
926 z1_is_zero = felem_is_zero(z1);
927 z2_is_zero = felem_is_zero(z2);
928 /* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
929 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero)
931 point_double(x3, y3, z3, x1, y1, z1);
938 felem_mul(tmp, z1, z2);
939 felem_reduce(ftmp5, tmp);
943 /* special case z2 = 0 is handled later */
944 felem_assign(ftmp5, z1);
947 /* z_out = (z1^2*x2 - z2^2*x1)*(z1*z2) */
948 felem_mul(tmp, ftmp, ftmp5);
949 felem_reduce(z_out, tmp);
951 /* ftmp = (z1^2*x2 - z2^2*x1)^2 */
952 felem_assign(ftmp5, ftmp);
953 felem_square(tmp, ftmp);
954 felem_reduce(ftmp, tmp);
956 /* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
957 felem_mul(tmp, ftmp, ftmp5);
958 felem_reduce(ftmp5, tmp);
960 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
961 felem_mul(tmp, ftmp2, ftmp);
962 felem_reduce(ftmp2, tmp);
964 /* tmp = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
965 felem_mul(tmp, ftmp4, ftmp5);
966 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
968 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
969 felem_square(tmp2, ftmp3);
970 /* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
972 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
973 felem_diff_128_64(tmp2, ftmp5);
974 /* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
976 /* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
977 felem_assign(ftmp5, ftmp2);
978 felem_scalar(ftmp5, 2);
979 /* ftmp5[i] < 2 * 2^57 = 2^58 */
982 * x_out = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
983 * 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2
985 felem_diff_128_64(tmp2, ftmp5);
986 /* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
987 felem_reduce(x_out, tmp2);
989 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out */
990 felem_diff(ftmp2, x_out);
991 /* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
993 /* tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) */
994 felem_mul(tmp2, ftmp3, ftmp2);
995 /* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
998 * y_out = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) -
999 * z2^3*y1*(z1^2*x2 - z2^2*x1)^3
1001 widefelem_diff(tmp2, tmp);
1002 /* tmp2[i] < 2^118 + 2^120 < 2^121 */
1003 felem_reduce(y_out, tmp2);
1005 /* the result (x_out, y_out, z_out) is incorrect if one of the inputs is
1006 * the point at infinity, so we need to check for this separately */
1008 /* if point 1 is at infinity, copy point 2 to output, and vice versa */
1009 copy_conditional(x_out, x2, z1_is_zero);
1010 copy_conditional(x_out, x1, z2_is_zero);
1011 copy_conditional(y_out, y2, z1_is_zero);
1012 copy_conditional(y_out, y1, z2_is_zero);
1013 copy_conditional(z_out, z2, z1_is_zero);
1014 copy_conditional(z_out, z1, z2_is_zero);
1015 felem_assign(x3, x_out);
1016 felem_assign(y3, y_out);
1017 felem_assign(z3, z_out);
1021 * select_point selects the |idx|th point from a precomputation table and
1023 * The pre_comp array argument should be size of |size| argument
1025 static void select_point(const u64 idx, unsigned int size, const felem pre_comp[][3], felem out[3])
1028 limb *outlimbs = &out[0][0];
1029 memset(outlimbs, 0, 3 * sizeof(felem));
1031 for (i = 0; i < size; i++)
1033 const limb *inlimbs = &pre_comp[i][0][0];
1040 for (j = 0; j < 4 * 3; j++)
1041 outlimbs[j] |= inlimbs[j] & mask;
1045 /* get_bit returns the |i|th bit in |in| */
1046 static char get_bit(const felem_bytearray in, unsigned i)
1050 return (in[i >> 3] >> (i & 7)) & 1;
1053 /* Interleaved point multiplication using precomputed point multiples:
1054 * The small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[],
1055 * the scalars in scalars[]. If g_scalar is non-NULL, we also add this multiple
1056 * of the generator, using certain (large) precomputed multiples in g_pre_comp.
1057 * Output point (X, Y, Z) is stored in x_out, y_out, z_out */
1058 static void batch_mul(felem x_out, felem y_out, felem z_out,
1059 const felem_bytearray scalars[], const unsigned num_points, const u8 *g_scalar,
1060 const int mixed, const felem pre_comp[][17][3], const felem g_pre_comp[2][16][3])
1064 unsigned gen_mul = (g_scalar != NULL);
1065 felem nq[3], tmp[4];
1069 /* set nq to the point at infinity */
1070 memset(nq, 0, 3 * sizeof(felem));
1072 /* Loop over all scalars msb-to-lsb, interleaving additions
1073 * of multiples of the generator (two in each of the last 28 rounds)
1074 * and additions of other points multiples (every 5th round).
1076 skip = 1; /* save two point operations in the first round */
1077 for (i = (num_points ? 220 : 27); i >= 0; --i)
1081 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1083 /* add multiples of the generator */
1084 if (gen_mul && (i <= 27))
1086 /* first, look 28 bits upwards */
1087 bits = get_bit(g_scalar, i + 196) << 3;
1088 bits |= get_bit(g_scalar, i + 140) << 2;
1089 bits |= get_bit(g_scalar, i + 84) << 1;
1090 bits |= get_bit(g_scalar, i + 28);
1091 /* select the point to add, in constant time */
1092 select_point(bits, 16, g_pre_comp[1], tmp);
1096 /* value 1 below is argument for "mixed" */
1097 point_add(nq[0], nq[1], nq[2],
1098 nq[0], nq[1], nq[2],
1099 1, tmp[0], tmp[1], tmp[2]);
1103 memcpy(nq, tmp, 3 * sizeof(felem));
1107 /* second, look at the current position */
1108 bits = get_bit(g_scalar, i + 168) << 3;
1109 bits |= get_bit(g_scalar, i + 112) << 2;
1110 bits |= get_bit(g_scalar, i + 56) << 1;
1111 bits |= get_bit(g_scalar, i);
1112 /* select the point to add, in constant time */
1113 select_point(bits, 16, g_pre_comp[0], tmp);
1114 point_add(nq[0], nq[1], nq[2],
1115 nq[0], nq[1], nq[2],
1116 1 /* mixed */, tmp[0], tmp[1], tmp[2]);
1119 /* do other additions every 5 doublings */
1120 if (num_points && (i % 5 == 0))
1122 /* loop over all scalars */
1123 for (num = 0; num < num_points; ++num)
1125 bits = get_bit(scalars[num], i + 4) << 5;
1126 bits |= get_bit(scalars[num], i + 3) << 4;
1127 bits |= get_bit(scalars[num], i + 2) << 3;
1128 bits |= get_bit(scalars[num], i + 1) << 2;
1129 bits |= get_bit(scalars[num], i) << 1;
1130 bits |= get_bit(scalars[num], i - 1);
1131 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1133 /* select the point to add or subtract */
1134 select_point(digit, 17, pre_comp[num], tmp);
1135 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative point */
1136 copy_conditional(tmp[1], tmp[3], sign);
1140 point_add(nq[0], nq[1], nq[2],
1141 nq[0], nq[1], nq[2],
1142 mixed, tmp[0], tmp[1], tmp[2]);
1146 memcpy(nq, tmp, 3 * sizeof(felem));
1152 felem_assign(x_out, nq[0]);
1153 felem_assign(y_out, nq[1]);
1154 felem_assign(z_out, nq[2]);
1157 /******************************************************************************/
1158 /* FUNCTIONS TO MANAGE PRECOMPUTATION
1161 static NISTP224_PRE_COMP *nistp224_pre_comp_new()
1163 NISTP224_PRE_COMP *ret = NULL;
1164 ret = (NISTP224_PRE_COMP *) OPENSSL_malloc(sizeof *ret);
1167 ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1170 memset(ret->g_pre_comp, 0, sizeof(ret->g_pre_comp));
1171 ret->references = 1;
1175 static void *nistp224_pre_comp_dup(void *src_)
1177 NISTP224_PRE_COMP *src = src_;
1179 /* no need to actually copy, these objects never change! */
1180 CRYPTO_add(&src->references, 1, CRYPTO_LOCK_EC_PRE_COMP);
1185 static void nistp224_pre_comp_free(void *pre_)
1188 NISTP224_PRE_COMP *pre = pre_;
1193 i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1200 static void nistp224_pre_comp_clear_free(void *pre_)
1203 NISTP224_PRE_COMP *pre = pre_;
1208 i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1212 OPENSSL_cleanse(pre, sizeof *pre);
1216 /******************************************************************************/
1217 /* OPENSSL EC_METHOD FUNCTIONS
1220 int ec_GFp_nistp224_group_init(EC_GROUP *group)
1223 ret = ec_GFp_simple_group_init(group);
1224 group->a_is_minus3 = 1;
1228 int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1229 const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
1232 BN_CTX *new_ctx = NULL;
1233 BIGNUM *curve_p, *curve_a, *curve_b;
1236 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1238 if (((curve_p = BN_CTX_get(ctx)) == NULL) ||
1239 ((curve_a = BN_CTX_get(ctx)) == NULL) ||
1240 ((curve_b = BN_CTX_get(ctx)) == NULL)) goto err;
1241 BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
1242 BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
1243 BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
1244 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) ||
1245 (BN_cmp(curve_b, b)))
1247 ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
1248 EC_R_WRONG_CURVE_PARAMETERS);
1251 group->field_mod_func = BN_nist_mod_224;
1252 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1255 if (new_ctx != NULL)
1256 BN_CTX_free(new_ctx);
1260 /* Takes the Jacobian coordinates (X, Y, Z) of a point and returns
1261 * (X', Y') = (X/Z^2, Y/Z^3) */
1262 int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
1263 const EC_POINT *point, BIGNUM *x, BIGNUM *y, BN_CTX *ctx)
1265 felem z1, z2, x_in, y_in, x_out, y_out;
1268 if (EC_POINT_is_at_infinity(group, point))
1270 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1271 EC_R_POINT_AT_INFINITY);
1274 if ((!BN_to_felem(x_in, &point->X)) || (!BN_to_felem(y_in, &point->Y)) ||
1275 (!BN_to_felem(z1, &point->Z))) return 0;
1277 felem_square(tmp, z2); felem_reduce(z1, tmp);
1278 felem_mul(tmp, x_in, z1); felem_reduce(x_in, tmp);
1279 felem_contract(x_out, x_in);
1282 if (!felem_to_BN(x, x_out)) {
1283 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1288 felem_mul(tmp, z1, z2); felem_reduce(z1, tmp);
1289 felem_mul(tmp, y_in, z1); felem_reduce(y_in, tmp);
1290 felem_contract(y_out, y_in);
1293 if (!felem_to_BN(y, y_out)) {
1294 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1302 static void make_points_affine(size_t num, felem points[/*num*/][3], felem tmp_felems[/*num+1*/])
1304 /* Runs in constant time, unless an input is the point at infinity
1305 * (which normally shouldn't happen). */
1306 ec_GFp_nistp_points_make_affine_internal(
1311 (void (*)(void *)) felem_one,
1312 (int (*)(const void *)) felem_is_zero_int,
1313 (void (*)(void *, const void *)) felem_assign,
1314 (void (*)(void *, const void *)) felem_square_reduce,
1315 (void (*)(void *, const void *, const void *)) felem_mul_reduce,
1316 (void (*)(void *, const void *)) felem_inv,
1317 (void (*)(void *, const void *)) felem_contract);
1320 /* Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL values
1321 * Result is stored in r (r can equal one of the inputs). */
1322 int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
1323 const BIGNUM *scalar, size_t num, const EC_POINT *points[],
1324 const BIGNUM *scalars[], BN_CTX *ctx)
1330 BN_CTX *new_ctx = NULL;
1331 BIGNUM *x, *y, *z, *tmp_scalar;
1332 felem_bytearray g_secret;
1333 felem_bytearray *secrets = NULL;
1334 felem (*pre_comp)[17][3] = NULL;
1335 felem *tmp_felems = NULL;
1336 felem_bytearray tmp;
1338 int have_pre_comp = 0;
1339 size_t num_points = num;
1340 felem x_in, y_in, z_in, x_out, y_out, z_out;
1341 NISTP224_PRE_COMP *pre = NULL;
1342 const felem (*g_pre_comp)[16][3] = NULL;
1343 EC_POINT *generator = NULL;
1344 const EC_POINT *p = NULL;
1345 const BIGNUM *p_scalar = NULL;
1348 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1350 if (((x = BN_CTX_get(ctx)) == NULL) ||
1351 ((y = BN_CTX_get(ctx)) == NULL) ||
1352 ((z = BN_CTX_get(ctx)) == NULL) ||
1353 ((tmp_scalar = BN_CTX_get(ctx)) == NULL))
1358 pre = EC_EX_DATA_get_data(group->extra_data,
1359 nistp224_pre_comp_dup, nistp224_pre_comp_free,
1360 nistp224_pre_comp_clear_free);
1362 /* we have precomputation, try to use it */
1363 g_pre_comp = (const felem (*)[16][3]) pre->g_pre_comp;
1365 /* try to use the standard precomputation */
1366 g_pre_comp = &gmul[0];
1367 generator = EC_POINT_new(group);
1368 if (generator == NULL)
1370 /* get the generator from precomputation */
1371 if (!felem_to_BN(x, g_pre_comp[0][1][0]) ||
1372 !felem_to_BN(y, g_pre_comp[0][1][1]) ||
1373 !felem_to_BN(z, g_pre_comp[0][1][2]))
1375 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1378 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1379 generator, x, y, z, ctx))
1381 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1382 /* precomputation matches generator */
1385 /* we don't have valid precomputation:
1386 * treat the generator as a random point */
1387 num_points = num_points + 1;
1392 if (num_points >= 3)
1394 /* unless we precompute multiples for just one or two points,
1395 * converting those into affine form is time well spent */
1398 secrets = OPENSSL_malloc(num_points * sizeof(felem_bytearray));
1399 pre_comp = OPENSSL_malloc(num_points * 17 * 3 * sizeof(felem));
1401 tmp_felems = OPENSSL_malloc((num_points * 17 + 1) * sizeof(felem));
1402 if ((secrets == NULL) || (pre_comp == NULL) || (mixed && (tmp_felems == NULL)))
1404 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1408 /* we treat NULL scalars as 0, and NULL points as points at infinity,
1409 * i.e., they contribute nothing to the linear combination */
1410 memset(secrets, 0, num_points * sizeof(felem_bytearray));
1411 memset(pre_comp, 0, num_points * 17 * 3 * sizeof(felem));
1412 for (i = 0; i < num_points; ++i)
1417 p = EC_GROUP_get0_generator(group);
1421 /* the i^th point */
1424 p_scalar = scalars[i];
1426 if ((p_scalar != NULL) && (p != NULL))
1428 /* reduce scalar to 0 <= scalar < 2^224 */
1429 if ((BN_num_bits(p_scalar) > 224) || (BN_is_negative(p_scalar)))
1431 /* this is an unusual input, and we don't guarantee
1432 * constant-timeness */
1433 if (!BN_nnmod(tmp_scalar, p_scalar, &group->order, ctx))
1435 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1438 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1441 num_bytes = BN_bn2bin(p_scalar, tmp);
1442 flip_endian(secrets[i], tmp, num_bytes);
1443 /* precompute multiples */
1444 if ((!BN_to_felem(x_out, &p->X)) ||
1445 (!BN_to_felem(y_out, &p->Y)) ||
1446 (!BN_to_felem(z_out, &p->Z))) goto err;
1447 felem_assign(pre_comp[i][1][0], x_out);
1448 felem_assign(pre_comp[i][1][1], y_out);
1449 felem_assign(pre_comp[i][1][2], z_out);
1450 for (j = 2; j <= 16; ++j)
1455 pre_comp[i][j][0], pre_comp[i][j][1], pre_comp[i][j][2],
1456 pre_comp[i][1][0], pre_comp[i][1][1], pre_comp[i][1][2],
1457 0, pre_comp[i][j-1][0], pre_comp[i][j-1][1], pre_comp[i][j-1][2]);
1462 pre_comp[i][j][0], pre_comp[i][j][1], pre_comp[i][j][2],
1463 pre_comp[i][j/2][0], pre_comp[i][j/2][1], pre_comp[i][j/2][2]);
1469 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
1472 /* the scalar for the generator */
1473 if ((scalar != NULL) && (have_pre_comp))
1475 memset(g_secret, 0, sizeof g_secret);
1476 /* reduce scalar to 0 <= scalar < 2^224 */
1477 if ((BN_num_bits(scalar) > 224) || (BN_is_negative(scalar)))
1479 /* this is an unusual input, and we don't guarantee
1480 * constant-timeness */
1481 if (!BN_nnmod(tmp_scalar, scalar, &group->order, ctx))
1483 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1486 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1489 num_bytes = BN_bn2bin(scalar, tmp);
1490 flip_endian(g_secret, tmp, num_bytes);
1491 /* do the multiplication with generator precomputation*/
1492 batch_mul(x_out, y_out, z_out,
1493 (const felem_bytearray (*)) secrets, num_points,
1495 mixed, (const felem (*)[17][3]) pre_comp,
1499 /* do the multiplication without generator precomputation */
1500 batch_mul(x_out, y_out, z_out,
1501 (const felem_bytearray (*)) secrets, num_points,
1502 NULL, mixed, (const felem (*)[17][3]) pre_comp, NULL);
1503 /* reduce the output to its unique minimal representation */
1504 felem_contract(x_in, x_out);
1505 felem_contract(y_in, y_out);
1506 felem_contract(z_in, z_out);
1507 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
1508 (!felem_to_BN(z, z_in)))
1510 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1513 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
1517 if (generator != NULL)
1518 EC_POINT_free(generator);
1519 if (new_ctx != NULL)
1520 BN_CTX_free(new_ctx);
1521 if (secrets != NULL)
1522 OPENSSL_free(secrets);
1523 if (pre_comp != NULL)
1524 OPENSSL_free(pre_comp);
1525 if (tmp_felems != NULL)
1526 OPENSSL_free(tmp_felems);
1530 int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
1533 NISTP224_PRE_COMP *pre = NULL;
1535 BN_CTX *new_ctx = NULL;
1537 EC_POINT *generator = NULL;
1538 felem tmp_felems[32];
1540 /* throw away old precomputation */
1541 EC_EX_DATA_free_data(&group->extra_data, nistp224_pre_comp_dup,
1542 nistp224_pre_comp_free, nistp224_pre_comp_clear_free);
1544 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1546 if (((x = BN_CTX_get(ctx)) == NULL) ||
1547 ((y = BN_CTX_get(ctx)) == NULL))
1549 /* get the generator */
1550 if (group->generator == NULL) goto err;
1551 generator = EC_POINT_new(group);
1552 if (generator == NULL)
1554 BN_bin2bn(nistp224_curve_params[3], sizeof (felem_bytearray), x);
1555 BN_bin2bn(nistp224_curve_params[4], sizeof (felem_bytearray), y);
1556 if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
1558 if ((pre = nistp224_pre_comp_new()) == NULL)
1560 /* if the generator is the standard one, use built-in precomputation */
1561 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1563 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
1567 if ((!BN_to_felem(pre->g_pre_comp[0][1][0], &group->generator->X)) ||
1568 (!BN_to_felem(pre->g_pre_comp[0][1][1], &group->generator->Y)) ||
1569 (!BN_to_felem(pre->g_pre_comp[0][1][2], &group->generator->Z)))
1571 /* compute 2^56*G, 2^112*G, 2^168*G for the first table,
1572 * 2^28*G, 2^84*G, 2^140*G, 2^196*G for the second one
1574 for (i = 1; i <= 8; i <<= 1)
1577 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2],
1578 pre->g_pre_comp[0][i][0], pre->g_pre_comp[0][i][1], pre->g_pre_comp[0][i][2]);
1579 for (j = 0; j < 27; ++j)
1582 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2],
1583 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1588 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2],
1589 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1590 for (j = 0; j < 27; ++j)
1593 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2],
1594 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2]);
1597 for (i = 0; i < 2; i++)
1599 /* g_pre_comp[i][0] is the point at infinity */
1600 memset(pre->g_pre_comp[i][0], 0, sizeof(pre->g_pre_comp[i][0]));
1601 /* the remaining multiples */
1602 /* 2^56*G + 2^112*G resp. 2^84*G + 2^140*G */
1604 pre->g_pre_comp[i][6][0], pre->g_pre_comp[i][6][1],
1605 pre->g_pre_comp[i][6][2], pre->g_pre_comp[i][4][0],
1606 pre->g_pre_comp[i][4][1], pre->g_pre_comp[i][4][2],
1607 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1608 pre->g_pre_comp[i][2][2]);
1609 /* 2^56*G + 2^168*G resp. 2^84*G + 2^196*G */
1611 pre->g_pre_comp[i][10][0], pre->g_pre_comp[i][10][1],
1612 pre->g_pre_comp[i][10][2], pre->g_pre_comp[i][8][0],
1613 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1614 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1615 pre->g_pre_comp[i][2][2]);
1616 /* 2^112*G + 2^168*G resp. 2^140*G + 2^196*G */
1618 pre->g_pre_comp[i][12][0], pre->g_pre_comp[i][12][1],
1619 pre->g_pre_comp[i][12][2], pre->g_pre_comp[i][8][0],
1620 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1621 0, pre->g_pre_comp[i][4][0], pre->g_pre_comp[i][4][1],
1622 pre->g_pre_comp[i][4][2]);
1623 /* 2^56*G + 2^112*G + 2^168*G resp. 2^84*G + 2^140*G + 2^196*G */
1625 pre->g_pre_comp[i][14][0], pre->g_pre_comp[i][14][1],
1626 pre->g_pre_comp[i][14][2], pre->g_pre_comp[i][12][0],
1627 pre->g_pre_comp[i][12][1], pre->g_pre_comp[i][12][2],
1628 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1629 pre->g_pre_comp[i][2][2]);
1630 for (j = 1; j < 8; ++j)
1632 /* odd multiples: add G resp. 2^28*G */
1634 pre->g_pre_comp[i][2*j+1][0], pre->g_pre_comp[i][2*j+1][1],
1635 pre->g_pre_comp[i][2*j+1][2], pre->g_pre_comp[i][2*j][0],
1636 pre->g_pre_comp[i][2*j][1], pre->g_pre_comp[i][2*j][2],
1637 0, pre->g_pre_comp[i][1][0], pre->g_pre_comp[i][1][1],
1638 pre->g_pre_comp[i][1][2]);
1641 make_points_affine(31, &(pre->g_pre_comp[0][1]), tmp_felems);
1643 if (!EC_EX_DATA_set_data(&group->extra_data, pre, nistp224_pre_comp_dup,
1644 nistp224_pre_comp_free, nistp224_pre_comp_clear_free))
1650 if (generator != NULL)
1651 EC_POINT_free(generator);
1652 if (new_ctx != NULL)
1653 BN_CTX_free(new_ctx);
1655 nistp224_pre_comp_free(pre);
1659 int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
1661 if (EC_EX_DATA_get_data(group->extra_data, nistp224_pre_comp_dup,
1662 nistp224_pre_comp_free, nistp224_pre_comp_clear_free)
1670 static void *dummy=&dummy;