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 /******************************************************************************/
49 /* INTERNAL REPRESENTATION OF FIELD ELEMENTS
51 * Field elements are represented as a_0 + 2^56*a_1 + 2^112*a_2 + 2^168*a_3
52 * using 64-bit coefficients called 'limbs',
53 * and sometimes (for multiplication results) as
54 * 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
55 * using 128-bit coefficients called 'widelimbs'.
56 * A 4-limb representation is an 'felem';
57 * a 7-widelimb representation is a 'widefelem'.
58 * Even within felems, bits of adjacent limbs overlap, and we don't always
59 * reduce the representations: we ensure that inputs to each felem
60 * multiplication satisfy a_i < 2^60, so outputs satisfy b_i < 4*2^60*2^60,
61 * and fit into a 128-bit word without overflow. The coefficients are then
62 * again partially reduced to obtain an felem satisfying a_i < 2^57.
63 * We only reduce to the unique minimal representation at the end of the
67 typedef uint64_t limb;
68 typedef uint128_t widelimb;
70 typedef limb felem[4];
71 typedef widelimb widefelem[7];
73 /* Field element represented as a byte arrary.
74 * 28*8 = 224 bits is also the group order size for the elliptic curve,
75 * and we also use this type for scalars for point multiplication.
77 typedef u8 felem_bytearray[28];
79 static const felem_bytearray nistp224_curve_params[5] = {
80 {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF, /* p */
81 0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0x00,0x00,0x00,0x00,
82 0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x01},
83 {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF, /* a */
84 0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,0xFF,0xFF,0xFF,0xFF,
85 0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFE},
86 {0xB4,0x05,0x0A,0x85,0x0C,0x04,0xB3,0xAB,0xF5,0x41, /* b */
87 0x32,0x56,0x50,0x44,0xB0,0xB7,0xD7,0xBF,0xD8,0xBA,
88 0x27,0x0B,0x39,0x43,0x23,0x55,0xFF,0xB4},
89 {0xB7,0x0E,0x0C,0xBD,0x6B,0xB4,0xBF,0x7F,0x32,0x13, /* x */
90 0x90,0xB9,0x4A,0x03,0xC1,0xD3,0x56,0xC2,0x11,0x22,
91 0x34,0x32,0x80,0xD6,0x11,0x5C,0x1D,0x21},
92 {0xbd,0x37,0x63,0x88,0xb5,0xf7,0x23,0xfb,0x4c,0x22, /* y */
93 0xdf,0xe6,0xcd,0x43,0x75,0xa0,0x5a,0x07,0x47,0x64,
94 0x44,0xd5,0x81,0x99,0x85,0x00,0x7e,0x34}
97 /* Precomputed multiples of the standard generator
98 * Points are given in coordinates (X, Y, Z) where Z normally is 1
99 * (0 for the point at infinity).
100 * For each field element, slice a_0 is word 0, etc.
102 * The table has 2 * 16 elements, starting with the following:
103 * index | bits | point
104 * ------+---------+------------------------------
107 * 2 | 0 0 1 0 | 2^56G
108 * 3 | 0 0 1 1 | (2^56 + 1)G
109 * 4 | 0 1 0 0 | 2^112G
110 * 5 | 0 1 0 1 | (2^112 + 1)G
111 * 6 | 0 1 1 0 | (2^112 + 2^56)G
112 * 7 | 0 1 1 1 | (2^112 + 2^56 + 1)G
113 * 8 | 1 0 0 0 | 2^168G
114 * 9 | 1 0 0 1 | (2^168 + 1)G
115 * 10 | 1 0 1 0 | (2^168 + 2^56)G
116 * 11 | 1 0 1 1 | (2^168 + 2^56 + 1)G
117 * 12 | 1 1 0 0 | (2^168 + 2^112)G
118 * 13 | 1 1 0 1 | (2^168 + 2^112 + 1)G
119 * 14 | 1 1 1 0 | (2^168 + 2^112 + 2^56)G
120 * 15 | 1 1 1 1 | (2^168 + 2^112 + 2^56 + 1)G
121 * followed by a copy of this with each element multiplied by 2^28.
123 * The reason for this is so that we can clock bits into four different
124 * locations when doing simple scalar multiplies against the base point,
125 * and then another four locations using the second 16 elements.
127 static const felem gmul[2][16][3] =
131 {{0x3280d6115c1d21, 0xc1d356c2112234, 0x7f321390b94a03, 0xb70e0cbd6bb4bf},
132 {0xd5819985007e34, 0x75a05a07476444, 0xfb4c22dfe6cd43, 0xbd376388b5f723},
134 {{0xfd9675666ebbe9, 0xbca7664d40ce5e, 0x2242df8d8a2a43, 0x1f49bbb0f99bc5},
135 {0x29e0b892dc9c43, 0xece8608436e662, 0xdc858f185310d0, 0x9812dd4eb8d321},
137 {{0x6d3e678d5d8eb8, 0x559eed1cb362f1, 0x16e9a3bbce8a3f, 0xeedcccd8c2a748},
138 {0xf19f90ed50266d, 0xabf2b4bf65f9df, 0x313865468fafec, 0x5cb379ba910a17},
140 {{0x0641966cab26e3, 0x91fb2991fab0a0, 0xefec27a4e13a0b, 0x0499aa8a5f8ebe},
141 {0x7510407766af5d, 0x84d929610d5450, 0x81d77aae82f706, 0x6916f6d4338c5b},
143 {{0xea95ac3b1f15c6, 0x086000905e82d4, 0xdd323ae4d1c8b1, 0x932b56be7685a3},
144 {0x9ef93dea25dbbf, 0x41665960f390f0, 0xfdec76dbe2a8a7, 0x523e80f019062a},
146 {{0x822fdd26732c73, 0xa01c83531b5d0f, 0x363f37347c1ba4, 0xc391b45c84725c},
147 {0xbbd5e1b2d6ad24, 0xddfbcde19dfaec, 0xc393da7e222a7f, 0x1efb7890ede244},
149 {{0x4c9e90ca217da1, 0xd11beca79159bb, 0xff8d33c2c98b7c, 0x2610b39409f849},
150 {0x44d1352ac64da0, 0xcdbb7b2c46b4fb, 0x966c079b753c89, 0xfe67e4e820b112},
152 {{0xe28cae2df5312d, 0xc71b61d16f5c6e, 0x79b7619a3e7c4c, 0x05c73240899b47},
153 {0x9f7f6382c73e3a, 0x18615165c56bda, 0x641fab2116fd56, 0x72855882b08394},
155 {{0x0469182f161c09, 0x74a98ca8d00fb5, 0xb89da93489a3e0, 0x41c98768fb0c1d},
156 {0xe5ea05fb32da81, 0x3dce9ffbca6855, 0x1cfe2d3fbf59e6, 0x0e5e03408738a7},
158 {{0xdab22b2333e87f, 0x4430137a5dd2f6, 0xe03ab9f738beb8, 0xcb0c5d0dc34f24},
159 {0x764a7df0c8fda5, 0x185ba5c3fa2044, 0x9281d688bcbe50, 0xc40331df893881},
161 {{0xb89530796f0f60, 0xade92bd26909a3, 0x1a0c83fb4884da, 0x1765bf22a5a984},
162 {0x772a9ee75db09e, 0x23bc6c67cec16f, 0x4c1edba8b14e2f, 0xe2a215d9611369},
164 {{0x571e509fb5efb3, 0xade88696410552, 0xc8ae85fada74fe, 0x6c7e4be83bbde3},
165 {0xff9f51160f4652, 0xb47ce2495a6539, 0xa2946c53b582f4, 0x286d2db3ee9a60},
167 {{0x40bbd5081a44af, 0x0995183b13926c, 0xbcefba6f47f6d0, 0x215619e9cc0057},
168 {0x8bc94d3b0df45e, 0xf11c54a3694f6f, 0x8631b93cdfe8b5, 0xe7e3f4b0982db9},
170 {{0xb17048ab3e1c7b, 0xac38f36ff8a1d8, 0x1c29819435d2c6, 0xc813132f4c07e9},
171 {0x2891425503b11f, 0x08781030579fea, 0xf5426ba5cc9674, 0x1e28ebf18562bc},
173 {{0x9f31997cc864eb, 0x06cd91d28b5e4c, 0xff17036691a973, 0xf1aef351497c58},
174 {0xdd1f2d600564ff, 0xdead073b1402db, 0x74a684435bd693, 0xeea7471f962558},
179 {{0x9665266dddf554, 0x9613d78b60ef2d, 0xce27a34cdba417, 0xd35ab74d6afc31},
180 {0x85ccdd22deb15e, 0x2137e5783a6aab, 0xa141cffd8c93c6, 0x355a1830e90f2d},
182 {{0x1a494eadaade65, 0xd6da4da77fe53c, 0xe7992996abec86, 0x65c3553c6090e3},
183 {0xfa610b1fb09346, 0xf1c6540b8a4aaf, 0xc51a13ccd3cbab, 0x02995b1b18c28a},
185 {{0x7874568e7295ef, 0x86b419fbe38d04, 0xdc0690a7550d9a, 0xd3966a44beac33},
186 {0x2b7280ec29132f, 0xbeaa3b6a032df3, 0xdc7dd88ae41200, 0xd25e2513e3a100},
188 {{0x924857eb2efafd, 0xac2bce41223190, 0x8edaa1445553fc, 0x825800fd3562d5},
189 {0x8d79148ea96621, 0x23a01c3dd9ed8d, 0xaf8b219f9416b5, 0xd8db0cc277daea},
191 {{0x76a9c3b1a700f0, 0xe9acd29bc7e691, 0x69212d1a6b0327, 0x6322e97fe154be},
192 {0x469fc5465d62aa, 0x8d41ed18883b05, 0x1f8eae66c52b88, 0xe4fcbe9325be51},
194 {{0x825fdf583cac16, 0x020b857c7b023a, 0x683c17744b0165, 0x14ffd0a2daf2f1},
195 {0x323b36184218f9, 0x4944ec4e3b47d4, 0xc15b3080841acf, 0x0bced4b01a28bb},
197 {{0x92ac22230df5c4, 0x52f33b4063eda8, 0xcb3f19870c0c93, 0x40064f2ba65233},
198 {0xfe16f0924f8992, 0x012da25af5b517, 0x1a57bb24f723a6, 0x06f8bc76760def},
200 {{0x4a7084f7817cb9, 0xbcab0738ee9a78, 0x3ec11e11d9c326, 0xdc0fe90e0f1aae},
201 {0xcf639ea5f98390, 0x5c350aa22ffb74, 0x9afae98a4047b7, 0x956ec2d617fc45},
203 {{0x4306d648c1be6a, 0x9247cd8bc9a462, 0xf5595e377d2f2e, 0xbd1c3caff1a52e},
204 {0x045e14472409d0, 0x29f3e17078f773, 0x745a602b2d4f7d, 0x191837685cdfbb},
206 {{0x5b6ee254a8cb79, 0x4953433f5e7026, 0xe21faeb1d1def4, 0xc4c225785c09de},
207 {0x307ce7bba1e518, 0x31b125b1036db8, 0x47e91868839e8f, 0xc765866e33b9f3},
209 {{0x3bfece24f96906, 0x4794da641e5093, 0xde5df64f95db26, 0x297ecd89714b05},
210 {0x701bd3ebb2c3aa, 0x7073b4f53cb1d5, 0x13c5665658af16, 0x9895089d66fe58},
212 {{0x0fef05f78c4790, 0x2d773633b05d2e, 0x94229c3a951c94, 0xbbbd70df4911bb},
213 {0xb2c6963d2c1168, 0x105f47a72b0d73, 0x9fdf6111614080, 0x7b7e94b39e67b0},
215 {{0xad1a7d6efbe2b3, 0xf012482c0da69d, 0x6b3bdf12438345, 0x40d7558d7aa4d9},
216 {0x8a09fffb5c6d3d, 0x9a356e5d9ffd38, 0x5973f15f4f9b1c, 0xdcd5f59f63c3ea},
218 {{0xacf39f4c5ca7ab, 0x4c8071cc5fd737, 0xc64e3602cd1184, 0x0acd4644c9abba},
219 {0x6c011a36d8bf6e, 0xfecd87ba24e32a, 0x19f6f56574fad8, 0x050b204ced9405},
221 {{0xed4f1cae7d9a96, 0x5ceef7ad94c40a, 0x778e4a3bf3ef9b, 0x7405783dc3b55e},
222 {0x32477c61b6e8c6, 0xb46a97570f018b, 0x91176d0a7e95d1, 0x3df90fbc4c7d0e},
225 /* Precomputation for the group generator. */
227 felem g_pre_comp[2][16][3];
231 const EC_METHOD *EC_GFp_nistp224_method(void)
233 static const EC_METHOD ret = {
234 EC_FLAGS_DEFAULT_OCT,
235 NID_X9_62_prime_field,
236 ec_GFp_nistp224_group_init,
237 ec_GFp_simple_group_finish,
238 ec_GFp_simple_group_clear_finish,
239 ec_GFp_nist_group_copy,
240 ec_GFp_nistp224_group_set_curve,
241 ec_GFp_simple_group_get_curve,
242 ec_GFp_simple_group_get_degree,
243 ec_GFp_simple_group_check_discriminant,
244 ec_GFp_simple_point_init,
245 ec_GFp_simple_point_finish,
246 ec_GFp_simple_point_clear_finish,
247 ec_GFp_simple_point_copy,
248 ec_GFp_simple_point_set_to_infinity,
249 ec_GFp_simple_set_Jprojective_coordinates_GFp,
250 ec_GFp_simple_get_Jprojective_coordinates_GFp,
251 ec_GFp_simple_point_set_affine_coordinates,
252 ec_GFp_nistp224_point_get_affine_coordinates,
253 0 /* point_set_compressed_coordinates */,
258 ec_GFp_simple_invert,
259 ec_GFp_simple_is_at_infinity,
260 ec_GFp_simple_is_on_curve,
262 ec_GFp_simple_make_affine,
263 ec_GFp_simple_points_make_affine,
264 ec_GFp_nistp224_points_mul,
265 ec_GFp_nistp224_precompute_mult,
266 ec_GFp_nistp224_have_precompute_mult,
267 ec_GFp_nist_field_mul,
268 ec_GFp_nist_field_sqr,
270 0 /* field_encode */,
271 0 /* field_decode */,
272 0 /* field_set_to_one */ };
277 /* Helper functions to convert field elements to/from internal representation */
278 static void bin28_to_felem(felem out, const u8 in[28])
280 out[0] = *((const uint64_t *)(in)) & 0x00ffffffffffffff;
281 out[1] = (*((const uint64_t *)(in+7))) & 0x00ffffffffffffff;
282 out[2] = (*((const uint64_t *)(in+14))) & 0x00ffffffffffffff;
283 out[3] = (*((const uint64_t *)(in+21))) & 0x00ffffffffffffff;
286 static void felem_to_bin28(u8 out[28], const felem in)
289 for (i = 0; i < 7; ++i)
291 out[i] = in[0]>>(8*i);
292 out[i+7] = in[1]>>(8*i);
293 out[i+14] = in[2]>>(8*i);
294 out[i+21] = in[3]>>(8*i);
298 /* To preserve endianness when using BN_bn2bin and BN_bin2bn */
299 static void flip_endian(u8 *out, const u8 *in, unsigned len)
302 for (i = 0; i < len; ++i)
303 out[i] = in[len-1-i];
306 /* From OpenSSL BIGNUM to internal representation */
307 static int BN_to_felem(felem out, const BIGNUM *bn)
309 felem_bytearray b_in;
310 felem_bytearray b_out;
313 /* BN_bn2bin eats leading zeroes */
314 memset(b_out, 0, sizeof b_out);
315 num_bytes = BN_num_bytes(bn);
316 if (num_bytes > sizeof b_out)
318 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
321 if (BN_is_negative(bn))
323 ECerr(EC_F_BN_TO_FELEM, EC_R_BIGNUM_OUT_OF_RANGE);
326 num_bytes = BN_bn2bin(bn, b_in);
327 flip_endian(b_out, b_in, num_bytes);
328 bin28_to_felem(out, b_out);
332 /* From internal representation to OpenSSL BIGNUM */
333 static BIGNUM *felem_to_BN(BIGNUM *out, const felem in)
335 felem_bytearray b_in, b_out;
336 felem_to_bin28(b_in, in);
337 flip_endian(b_out, b_in, sizeof b_out);
338 return BN_bin2bn(b_out, sizeof b_out, out);
341 /******************************************************************************/
344 * Field operations, using the internal representation of field elements.
345 * NB! These operations are specific to our point multiplication and cannot be
346 * expected to be correct in general - e.g., multiplication with a large scalar
347 * will cause an overflow.
351 static void felem_one(felem out)
359 static void felem_assign(felem out, const felem in)
367 /* Sum two field elements: out += in */
368 static void felem_sum(felem out, const felem in)
376 /* Get negative value: out = -in */
377 /* Assumes in[i] < 2^57 */
378 static void felem_neg(felem out, const felem in)
380 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
381 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
382 static const limb two58m42m2 = (((limb) 1) << 58) -
383 (((limb) 1) << 42) - (((limb) 1) << 2);
385 /* Set to 0 mod 2^224-2^96+1 to ensure out > in */
386 out[0] = two58p2 - in[0];
387 out[1] = two58m42m2 - in[1];
388 out[2] = two58m2 - in[2];
389 out[3] = two58m2 - in[3];
392 /* Subtract field elements: out -= in */
393 /* Assumes in[i] < 2^57 */
394 static void felem_diff(felem out, const felem in)
396 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
397 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
398 static const limb two58m42m2 = (((limb) 1) << 58) -
399 (((limb) 1) << 42) - (((limb) 1) << 2);
401 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
403 out[1] += two58m42m2;
413 /* Subtract in unreduced 128-bit mode: out -= in */
414 /* Assumes in[i] < 2^119 */
415 static void widefelem_diff(widefelem out, const widefelem in)
417 static const widelimb two120 = ((widelimb) 1) << 120;
418 static const widelimb two120m64 = (((widelimb) 1) << 120) -
419 (((widelimb) 1) << 64);
420 static const widelimb two120m104m64 = (((widelimb) 1) << 120) -
421 (((widelimb) 1) << 104) - (((widelimb) 1) << 64);
423 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
428 out[4] += two120m104m64;
441 /* Subtract in mixed mode: out128 -= in64 */
443 static void felem_diff_128_64(widefelem out, const felem in)
445 static const widelimb two64p8 = (((widelimb) 1) << 64) +
446 (((widelimb) 1) << 8);
447 static const widelimb two64m8 = (((widelimb) 1) << 64) -
448 (((widelimb) 1) << 8);
449 static const widelimb two64m48m8 = (((widelimb) 1) << 64) -
450 (((widelimb) 1) << 48) - (((widelimb) 1) << 8);
452 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
454 out[1] += two64m48m8;
464 /* Multiply a field element by a scalar: out = out * scalar
465 * The scalars we actually use are small, so results fit without overflow */
466 static void felem_scalar(felem out, const limb scalar)
474 /* Multiply an unreduced field element by a scalar: out = out * scalar
475 * The scalars we actually use are small, so results fit without overflow */
476 static void widefelem_scalar(widefelem out, const widelimb scalar)
487 /* Square a field element: out = in^2 */
488 static void felem_square(widefelem out, const felem in)
490 limb tmp0, tmp1, tmp2;
491 tmp0 = 2 * in[0]; tmp1 = 2 * in[1]; tmp2 = 2 * in[2];
492 out[0] = ((widelimb) in[0]) * in[0];
493 out[1] = ((widelimb) in[0]) * tmp1;
494 out[2] = ((widelimb) in[0]) * tmp2 + ((widelimb) in[1]) * in[1];
495 out[3] = ((widelimb) in[3]) * tmp0 +
496 ((widelimb) in[1]) * tmp2;
497 out[4] = ((widelimb) in[3]) * tmp1 + ((widelimb) in[2]) * in[2];
498 out[5] = ((widelimb) in[3]) * tmp2;
499 out[6] = ((widelimb) in[3]) * in[3];
502 /* Multiply two field elements: out = in1 * in2 */
503 static void felem_mul(widefelem out, const felem in1, const felem in2)
505 out[0] = ((widelimb) in1[0]) * in2[0];
506 out[1] = ((widelimb) in1[0]) * in2[1] + ((widelimb) in1[1]) * in2[0];
507 out[2] = ((widelimb) in1[0]) * in2[2] + ((widelimb) in1[1]) * in2[1] +
508 ((widelimb) in1[2]) * in2[0];
509 out[3] = ((widelimb) in1[0]) * in2[3] + ((widelimb) in1[1]) * in2[2] +
510 ((widelimb) in1[2]) * in2[1] + ((widelimb) in1[3]) * in2[0];
511 out[4] = ((widelimb) in1[1]) * in2[3] + ((widelimb) in1[2]) * in2[2] +
512 ((widelimb) in1[3]) * in2[1];
513 out[5] = ((widelimb) in1[2]) * in2[3] + ((widelimb) in1[3]) * in2[2];
514 out[6] = ((widelimb) in1[3]) * in2[3];
517 /* Reduce seven 128-bit coefficients to four 64-bit coefficients.
518 * Requires in[i] < 2^126,
519 * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] <= 2^56 + 2^16 */
520 static void felem_reduce(felem out, const widefelem in)
522 static const widelimb two127p15 = (((widelimb) 1) << 127) +
523 (((widelimb) 1) << 15);
524 static const widelimb two127m71 = (((widelimb) 1) << 127) -
525 (((widelimb) 1) << 71);
526 static const widelimb two127m71m55 = (((widelimb) 1) << 127) -
527 (((widelimb) 1) << 71) - (((widelimb) 1) << 55);
530 /* Add 0 mod 2^224-2^96+1 to ensure all differences are positive */
531 output[0] = in[0] + two127p15;
532 output[1] = in[1] + two127m71m55;
533 output[2] = in[2] + two127m71;
537 /* Eliminate in[4], in[5], in[6] */
538 output[4] += in[6] >> 16;
539 output[3] += (in[6] & 0xffff) << 40;
542 output[3] += in[5] >> 16;
543 output[2] += (in[5] & 0xffff) << 40;
546 output[2] += output[4] >> 16;
547 output[1] += (output[4] & 0xffff) << 40;
548 output[0] -= output[4];
550 /* Carry 2 -> 3 -> 4 */
551 output[3] += output[2] >> 56;
552 output[2] &= 0x00ffffffffffffff;
554 output[4] = output[3] >> 56;
555 output[3] &= 0x00ffffffffffffff;
557 /* Now output[2] < 2^56, output[3] < 2^56, output[4] < 2^72 */
559 /* Eliminate output[4] */
560 output[2] += output[4] >> 16;
561 /* output[2] < 2^56 + 2^56 = 2^57 */
562 output[1] += (output[4] & 0xffff) << 40;
563 output[0] -= output[4];
565 /* Carry 0 -> 1 -> 2 -> 3 */
566 output[1] += output[0] >> 56;
567 out[0] = output[0] & 0x00ffffffffffffff;
569 output[2] += output[1] >> 56;
570 /* output[2] < 2^57 + 2^72 */
571 out[1] = output[1] & 0x00ffffffffffffff;
572 output[3] += output[2] >> 56;
573 /* output[3] <= 2^56 + 2^16 */
574 out[2] = output[2] & 0x00ffffffffffffff;
576 /* out[0] < 2^56, out[1] < 2^56, out[2] < 2^56,
577 * out[3] <= 2^56 + 2^16 (due to final carry),
582 static void felem_square_reduce(felem out, const felem in)
585 felem_square(tmp, in);
586 felem_reduce(out, tmp);
589 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
592 felem_mul(tmp, in1, in2);
593 felem_reduce(out, tmp);
596 /* Reduce to unique minimal representation.
597 * Requires 0 <= in < 2*p (always call felem_reduce first) */
598 static void felem_contract(felem out, const felem in)
600 static const int64_t two56 = ((limb) 1) << 56;
601 /* 0 <= in < 2*p, p = 2^224 - 2^96 + 1 */
602 /* if in > p , reduce in = in - 2^224 + 2^96 - 1 */
608 /* Case 1: a = 1 iff in >= 2^224 */
612 tmp[3] &= 0x00ffffffffffffff;
613 /* Case 2: a = 0 iff p <= in < 2^224, i.e.,
614 * the high 128 bits are all 1 and the lower part is non-zero */
615 a = ((in[3] & in[2] & (in[1] | 0x000000ffffffffff)) + 1) |
616 (((int64_t)(in[0] + (in[1] & 0x000000ffffffffff)) - 1) >> 63);
617 a &= 0x00ffffffffffffff;
618 /* turn a into an all-one mask (if a = 0) or an all-zero mask */
620 /* subtract 2^224 - 2^96 + 1 if a is all-one*/
621 tmp[3] &= a ^ 0xffffffffffffffff;
622 tmp[2] &= a ^ 0xffffffffffffffff;
623 tmp[1] &= (a ^ 0xffffffffffffffff) | 0x000000ffffffffff;
626 /* eliminate negative coefficients: if tmp[0] is negative, tmp[1] must
627 * be non-zero, so we only need one step */
632 /* carry 1 -> 2 -> 3 */
633 tmp[2] += tmp[1] >> 56;
634 tmp[1] &= 0x00ffffffffffffff;
636 tmp[3] += tmp[2] >> 56;
637 tmp[2] &= 0x00ffffffffffffff;
639 /* Now 0 <= out < p */
646 /* Zero-check: returns 1 if input is 0, and 0 otherwise.
647 * We know that field elements are reduced to in < 2^225,
648 * so we only need to check three cases: 0, 2^224 - 2^96 + 1,
649 * and 2^225 - 2^97 + 2 */
650 static limb felem_is_zero(const felem in)
652 limb zero, two224m96p1, two225m97p2;
654 zero = in[0] | in[1] | in[2] | in[3];
655 zero = (((int64_t)(zero) - 1) >> 63) & 1;
656 two224m96p1 = (in[0] ^ 1) | (in[1] ^ 0x00ffff0000000000)
657 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x00ffffffffffffff);
658 two224m96p1 = (((int64_t)(two224m96p1) - 1) >> 63) & 1;
659 two225m97p2 = (in[0] ^ 2) | (in[1] ^ 0x00fffe0000000000)
660 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x01ffffffffffffff);
661 two225m97p2 = (((int64_t)(two225m97p2) - 1) >> 63) & 1;
662 return (zero | two224m96p1 | two225m97p2);
665 static limb felem_is_zero_int(const felem in)
667 return (int) (felem_is_zero(in) & ((limb)1));
670 /* Invert a field element */
671 /* Computation chain copied from djb's code */
672 static void felem_inv(felem out, const felem in)
674 felem ftmp, ftmp2, ftmp3, ftmp4;
678 felem_square(tmp, in); felem_reduce(ftmp, tmp); /* 2 */
679 felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp); /* 2^2 - 1 */
680 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp); /* 2^3 - 2 */
681 felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp); /* 2^3 - 1 */
682 felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp); /* 2^4 - 2 */
683 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp); /* 2^5 - 4 */
684 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp); /* 2^6 - 8 */
685 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp); /* 2^6 - 1 */
686 felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp); /* 2^7 - 2 */
687 for (i = 0; i < 5; ++i) /* 2^12 - 2^6 */
689 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
691 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp2, tmp); /* 2^12 - 1 */
692 felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp); /* 2^13 - 2 */
693 for (i = 0; i < 11; ++i) /* 2^24 - 2^12 */
695 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
697 felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp2, tmp); /* 2^24 - 1 */
698 felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp); /* 2^25 - 2 */
699 for (i = 0; i < 23; ++i) /* 2^48 - 2^24 */
701 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
703 felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp3, tmp); /* 2^48 - 1 */
704 felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp); /* 2^49 - 2 */
705 for (i = 0; i < 47; ++i) /* 2^96 - 2^48 */
707 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
709 felem_mul(tmp, ftmp3, ftmp4); felem_reduce(ftmp3, tmp); /* 2^96 - 1 */
710 felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp); /* 2^97 - 2 */
711 for (i = 0; i < 23; ++i) /* 2^120 - 2^24 */
713 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
715 felem_mul(tmp, ftmp2, ftmp4); felem_reduce(ftmp2, tmp); /* 2^120 - 1 */
716 for (i = 0; i < 6; ++i) /* 2^126 - 2^6 */
718 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
720 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp); /* 2^126 - 1 */
721 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp); /* 2^127 - 2 */
722 felem_mul(tmp, ftmp, in); felem_reduce(ftmp, tmp); /* 2^127 - 1 */
723 for (i = 0; i < 97; ++i) /* 2^224 - 2^97 */
725 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp);
727 felem_mul(tmp, ftmp, ftmp3); felem_reduce(out, tmp); /* 2^224 - 2^96 - 1 */
730 /* Copy in constant time:
731 * if icopy == 1, copy in to out,
732 * if icopy == 0, copy out to itself. */
734 copy_conditional(felem out, const felem in, limb icopy)
737 /* icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one */
738 const limb copy = -icopy;
739 for (i = 0; i < 4; ++i)
741 const limb tmp = copy & (in[i] ^ out[i]);
746 /******************************************************************************/
747 /* ELLIPTIC CURVE POINT OPERATIONS
749 * Points are represented in Jacobian projective coordinates:
750 * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
751 * or to the point at infinity if Z == 0.
755 /* Double an elliptic curve point:
756 * (X', Y', Z') = 2 * (X, Y, Z), where
757 * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
758 * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^2
759 * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
760 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
761 * while x_out == y_in is not (maybe this works, but it's not tested). */
763 point_double(felem x_out, felem y_out, felem z_out,
764 const felem x_in, const felem y_in, const felem z_in)
767 felem delta, gamma, beta, alpha, ftmp, ftmp2;
769 felem_assign(ftmp, x_in);
770 felem_assign(ftmp2, x_in);
773 felem_square(tmp, z_in);
774 felem_reduce(delta, tmp);
777 felem_square(tmp, y_in);
778 felem_reduce(gamma, tmp);
781 felem_mul(tmp, x_in, gamma);
782 felem_reduce(beta, tmp);
784 /* alpha = 3*(x-delta)*(x+delta) */
785 felem_diff(ftmp, delta);
786 /* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
787 felem_sum(ftmp2, delta);
788 /* ftmp2[i] < 2^57 + 2^57 = 2^58 */
789 felem_scalar(ftmp2, 3);
790 /* ftmp2[i] < 3 * 2^58 < 2^60 */
791 felem_mul(tmp, ftmp, ftmp2);
792 /* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
793 felem_reduce(alpha, tmp);
795 /* x' = alpha^2 - 8*beta */
796 felem_square(tmp, alpha);
797 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
798 felem_assign(ftmp, beta);
799 felem_scalar(ftmp, 8);
800 /* ftmp[i] < 8 * 2^57 = 2^60 */
801 felem_diff_128_64(tmp, ftmp);
802 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
803 felem_reduce(x_out, tmp);
805 /* z' = (y + z)^2 - gamma - delta */
806 felem_sum(delta, gamma);
807 /* delta[i] < 2^57 + 2^57 = 2^58 */
808 felem_assign(ftmp, y_in);
809 felem_sum(ftmp, z_in);
810 /* ftmp[i] < 2^57 + 2^57 = 2^58 */
811 felem_square(tmp, ftmp);
812 /* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
813 felem_diff_128_64(tmp, delta);
814 /* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
815 felem_reduce(z_out, tmp);
817 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
818 felem_scalar(beta, 4);
819 /* beta[i] < 4 * 2^57 = 2^59 */
820 felem_diff(beta, x_out);
821 /* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
822 felem_mul(tmp, alpha, beta);
823 /* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
824 felem_square(tmp2, gamma);
825 /* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
826 widefelem_scalar(tmp2, 8);
827 /* tmp2[i] < 8 * 2^116 = 2^119 */
828 widefelem_diff(tmp, tmp2);
829 /* tmp[i] < 2^119 + 2^120 < 2^121 */
830 felem_reduce(y_out, tmp);
833 /* Add two elliptic curve points:
834 * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
835 * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
836 * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
837 * 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) -
838 * Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
839 * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2)
841 * This runs faster if 'mixed' is set, which requires Z_2 = 1 or Z_2 = 0.
844 /* This function is not entirely constant-time:
845 * it includes a branch for checking whether the two input points are equal,
846 * (while not equal to the point at infinity).
847 * This case never happens during single point multiplication,
848 * so there is no timing leak for ECDH or ECDSA signing. */
849 static void point_add(felem x3, felem y3, felem z3,
850 const felem x1, const felem y1, const felem z1,
851 const int mixed, const felem x2, const felem y2, const felem z2)
853 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, x_out, y_out, z_out;
855 limb z1_is_zero, z2_is_zero, x_equal, y_equal;
860 felem_square(tmp, z2);
861 felem_reduce(ftmp2, tmp);
864 felem_mul(tmp, ftmp2, z2);
865 felem_reduce(ftmp4, tmp);
867 /* ftmp4 = z2^3*y1 */
868 felem_mul(tmp2, ftmp4, y1);
869 felem_reduce(ftmp4, tmp2);
871 /* ftmp2 = z2^2*x1 */
872 felem_mul(tmp2, ftmp2, x1);
873 felem_reduce(ftmp2, tmp2);
877 /* We'll assume z2 = 1 (special case z2 = 0 is handled later) */
879 /* ftmp4 = z2^3*y1 */
880 felem_assign(ftmp4, y1);
882 /* ftmp2 = z2^2*x1 */
883 felem_assign(ftmp2, x1);
887 felem_square(tmp, z1);
888 felem_reduce(ftmp, tmp);
891 felem_mul(tmp, ftmp, z1);
892 felem_reduce(ftmp3, tmp);
895 felem_mul(tmp, ftmp3, y2);
896 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
898 /* ftmp3 = z1^3*y2 - z2^3*y1 */
899 felem_diff_128_64(tmp, ftmp4);
900 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
901 felem_reduce(ftmp3, tmp);
904 felem_mul(tmp, ftmp, x2);
905 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
907 /* ftmp = z1^2*x2 - z2^2*x1 */
908 felem_diff_128_64(tmp, ftmp2);
909 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
910 felem_reduce(ftmp, tmp);
912 /* the formulae are incorrect if the points are equal
913 * so we check for this and do doubling if this happens */
914 x_equal = felem_is_zero(ftmp);
915 y_equal = felem_is_zero(ftmp3);
916 z1_is_zero = felem_is_zero(z1);
917 z2_is_zero = felem_is_zero(z2);
918 /* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
919 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero)
921 point_double(x3, y3, z3, x1, y1, z1);
928 felem_mul(tmp, z1, z2);
929 felem_reduce(ftmp5, tmp);
933 /* special case z2 = 0 is handled later */
934 felem_assign(ftmp5, z1);
937 /* z_out = (z1^2*x2 - z2^2*x1)*(z1*z2) */
938 felem_mul(tmp, ftmp, ftmp5);
939 felem_reduce(z_out, tmp);
941 /* ftmp = (z1^2*x2 - z2^2*x1)^2 */
942 felem_assign(ftmp5, ftmp);
943 felem_square(tmp, ftmp);
944 felem_reduce(ftmp, tmp);
946 /* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
947 felem_mul(tmp, ftmp, ftmp5);
948 felem_reduce(ftmp5, tmp);
950 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
951 felem_mul(tmp, ftmp2, ftmp);
952 felem_reduce(ftmp2, tmp);
954 /* tmp = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
955 felem_mul(tmp, ftmp4, ftmp5);
956 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
958 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
959 felem_square(tmp2, ftmp3);
960 /* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
962 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
963 felem_diff_128_64(tmp2, ftmp5);
964 /* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
966 /* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
967 felem_assign(ftmp5, ftmp2);
968 felem_scalar(ftmp5, 2);
969 /* ftmp5[i] < 2 * 2^57 = 2^58 */
971 /* x_out = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
972 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
973 felem_diff_128_64(tmp2, ftmp5);
974 /* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
975 felem_reduce(x_out, tmp2);
977 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out */
978 felem_diff(ftmp2, x_out);
979 /* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
981 /* tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) */
982 felem_mul(tmp2, ftmp3, ftmp2);
983 /* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
985 /* y_out = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) -
986 z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
987 widefelem_diff(tmp2, tmp);
988 /* tmp2[i] < 2^118 + 2^120 < 2^121 */
989 felem_reduce(y_out, tmp2);
991 /* the result (x_out, y_out, z_out) is incorrect if one of the inputs is
992 * the point at infinity, so we need to check for this separately */
994 /* if point 1 is at infinity, copy point 2 to output, and vice versa */
995 copy_conditional(x_out, x2, z1_is_zero);
996 copy_conditional(x_out, x1, z2_is_zero);
997 copy_conditional(y_out, y2, z1_is_zero);
998 copy_conditional(y_out, y1, z2_is_zero);
999 copy_conditional(z_out, z2, z1_is_zero);
1000 copy_conditional(z_out, z1, z2_is_zero);
1001 felem_assign(x3, x_out);
1002 felem_assign(y3, y_out);
1003 felem_assign(z3, z_out);
1006 /* select_point selects the |idx|th point from a precomputation table and
1007 * copies it to out. */
1008 static void select_point(const u64 idx, unsigned int size, const felem pre_comp[/*size*/][3], felem out[3])
1011 limb *outlimbs = &out[0][0];
1012 memset(outlimbs, 0, 3 * sizeof(felem));
1014 for (i = 0; i < size; i++)
1016 const limb *inlimbs = &pre_comp[i][0][0];
1023 for (j = 0; j < 4 * 3; j++)
1024 outlimbs[j] |= inlimbs[j] & mask;
1028 /* get_bit returns the |i|th bit in |in| */
1029 static char get_bit(const felem_bytearray in, unsigned i)
1033 return (in[i >> 3] >> (i & 7)) & 1;
1036 /* Interleaved point multiplication using precomputed point multiples:
1037 * The small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[],
1038 * the scalars in scalars[]. If g_scalar is non-NULL, we also add this multiple
1039 * of the generator, using certain (large) precomputed multiples in g_pre_comp.
1040 * Output point (X, Y, Z) is stored in x_out, y_out, z_out */
1041 static void batch_mul(felem x_out, felem y_out, felem z_out,
1042 const felem_bytearray scalars[], const unsigned num_points, const u8 *g_scalar,
1043 const int mixed, const felem pre_comp[][17][3], const felem g_pre_comp[2][16][3])
1047 unsigned gen_mul = (g_scalar != NULL);
1048 felem nq[3], tmp[4];
1052 /* set nq to the point at infinity */
1053 memset(nq, 0, 3 * sizeof(felem));
1055 /* Loop over all scalars msb-to-lsb, interleaving additions
1056 * of multiples of the generator (two in each of the last 28 rounds)
1057 * and additions of other points multiples (every 5th round).
1059 skip = 1; /* save two point operations in the first round */
1060 for (i = (num_points ? 220 : 27); i >= 0; --i)
1064 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1066 /* add multiples of the generator */
1067 if (gen_mul && (i <= 27))
1069 /* first, look 28 bits upwards */
1070 bits = get_bit(g_scalar, i + 196) << 3;
1071 bits |= get_bit(g_scalar, i + 140) << 2;
1072 bits |= get_bit(g_scalar, i + 84) << 1;
1073 bits |= get_bit(g_scalar, i + 28);
1074 /* select the point to add, in constant time */
1075 select_point(bits, 16, g_pre_comp[1], tmp);
1079 point_add(nq[0], nq[1], nq[2],
1080 nq[0], nq[1], nq[2],
1081 1 /* mixed */, tmp[0], tmp[1], tmp[2]);
1085 memcpy(nq, tmp, 3 * sizeof(felem));
1089 /* second, look at the current position */
1090 bits = get_bit(g_scalar, i + 168) << 3;
1091 bits |= get_bit(g_scalar, i + 112) << 2;
1092 bits |= get_bit(g_scalar, i + 56) << 1;
1093 bits |= get_bit(g_scalar, i);
1094 /* select the point to add, in constant time */
1095 select_point(bits, 16, g_pre_comp[0], tmp);
1096 point_add(nq[0], nq[1], nq[2],
1097 nq[0], nq[1], nq[2],
1098 1 /* mixed */, tmp[0], tmp[1], tmp[2]);
1101 /* do other additions every 5 doublings */
1102 if (num_points && (i % 5 == 0))
1104 /* loop over all scalars */
1105 for (num = 0; num < num_points; ++num)
1107 bits = get_bit(scalars[num], i + 4) << 5;
1108 bits |= get_bit(scalars[num], i + 3) << 4;
1109 bits |= get_bit(scalars[num], i + 2) << 3;
1110 bits |= get_bit(scalars[num], i + 1) << 2;
1111 bits |= get_bit(scalars[num], i) << 1;
1112 bits |= get_bit(scalars[num], i - 1);
1113 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1115 /* select the point to add or subtract */
1116 select_point(digit, 17, pre_comp[num], tmp);
1117 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative point */
1118 copy_conditional(tmp[1], tmp[3], sign);
1122 point_add(nq[0], nq[1], nq[2],
1123 nq[0], nq[1], nq[2],
1124 mixed, tmp[0], tmp[1], tmp[2]);
1128 memcpy(nq, tmp, 3 * sizeof(felem));
1134 felem_assign(x_out, nq[0]);
1135 felem_assign(y_out, nq[1]);
1136 felem_assign(z_out, nq[2]);
1139 /******************************************************************************/
1140 /* FUNCTIONS TO MANAGE PRECOMPUTATION
1143 static NISTP224_PRE_COMP *nistp224_pre_comp_new()
1145 NISTP224_PRE_COMP *ret = NULL;
1146 ret = (NISTP224_PRE_COMP *) OPENSSL_malloc(sizeof *ret);
1149 ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1152 memset(ret->g_pre_comp, 0, sizeof(ret->g_pre_comp));
1153 ret->references = 1;
1157 static void *nistp224_pre_comp_dup(void *src_)
1159 NISTP224_PRE_COMP *src = src_;
1161 /* no need to actually copy, these objects never change! */
1162 CRYPTO_add(&src->references, 1, CRYPTO_LOCK_EC_PRE_COMP);
1167 static void nistp224_pre_comp_free(void *pre_)
1170 NISTP224_PRE_COMP *pre = pre_;
1175 i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1182 static void nistp224_pre_comp_clear_free(void *pre_)
1185 NISTP224_PRE_COMP *pre = pre_;
1190 i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1194 OPENSSL_cleanse(pre, sizeof *pre);
1198 /******************************************************************************/
1199 /* OPENSSL EC_METHOD FUNCTIONS
1202 int ec_GFp_nistp224_group_init(EC_GROUP *group)
1205 ret = ec_GFp_simple_group_init(group);
1206 group->a_is_minus3 = 1;
1210 int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1211 const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
1214 BN_CTX *new_ctx = NULL;
1215 BIGNUM *curve_p, *curve_a, *curve_b;
1218 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1220 if (((curve_p = BN_CTX_get(ctx)) == NULL) ||
1221 ((curve_a = BN_CTX_get(ctx)) == NULL) ||
1222 ((curve_b = BN_CTX_get(ctx)) == NULL)) goto err;
1223 BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
1224 BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
1225 BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
1226 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) ||
1227 (BN_cmp(curve_b, b)))
1229 ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
1230 EC_R_WRONG_CURVE_PARAMETERS);
1233 group->field_mod_func = BN_nist_mod_224;
1234 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1237 if (new_ctx != NULL)
1238 BN_CTX_free(new_ctx);
1242 /* Takes the Jacobian coordinates (X, Y, Z) of a point and returns
1243 * (X', Y') = (X/Z^2, Y/Z^3) */
1244 int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
1245 const EC_POINT *point, BIGNUM *x, BIGNUM *y, BN_CTX *ctx)
1247 felem z1, z2, x_in, y_in, x_out, y_out;
1250 if (EC_POINT_is_at_infinity(group, point))
1252 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1253 EC_R_POINT_AT_INFINITY);
1256 if ((!BN_to_felem(x_in, &point->X)) || (!BN_to_felem(y_in, &point->Y)) ||
1257 (!BN_to_felem(z1, &point->Z))) return 0;
1259 felem_square(tmp, z2); felem_reduce(z1, tmp);
1260 felem_mul(tmp, x_in, z1); felem_reduce(x_in, tmp);
1261 felem_contract(x_out, x_in);
1264 if (!felem_to_BN(x, x_out)) {
1265 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1270 felem_mul(tmp, z1, z2); felem_reduce(z1, tmp);
1271 felem_mul(tmp, y_in, z1); felem_reduce(y_in, tmp);
1272 felem_contract(y_out, y_in);
1275 if (!felem_to_BN(y, y_out)) {
1276 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1284 static void make_points_affine(size_t num, felem points[/*num*/][3], felem tmp_felems[/*num+1*/])
1286 /* Runs in constant time, unless an input is the point at infinity
1287 * (which normally shouldn't happen). */
1288 ec_GFp_nistp_points_make_affine_internal(
1293 (void (*)(void *)) felem_one,
1294 (int (*)(const void *)) felem_is_zero_int,
1295 (void (*)(void *, const void *)) felem_assign,
1296 (void (*)(void *, const void *)) felem_square_reduce,
1297 (void (*)(void *, const void *, const void *)) felem_mul_reduce,
1298 (void (*)(void *, const void *)) felem_inv,
1299 (void (*)(void *, const void *)) felem_contract);
1302 /* Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL values
1303 * Result is stored in r (r can equal one of the inputs). */
1304 int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
1305 const BIGNUM *scalar, size_t num, const EC_POINT *points[],
1306 const BIGNUM *scalars[], BN_CTX *ctx)
1312 BN_CTX *new_ctx = NULL;
1313 BIGNUM *x, *y, *z, *tmp_scalar;
1314 felem_bytearray g_secret;
1315 felem_bytearray *secrets = NULL;
1316 felem (*pre_comp)[17][3] = NULL;
1317 felem *tmp_felems = NULL;
1318 felem_bytearray tmp;
1320 int have_pre_comp = 0;
1321 size_t num_points = num;
1322 felem x_in, y_in, z_in, x_out, y_out, z_out;
1323 NISTP224_PRE_COMP *pre = NULL;
1324 const felem (*g_pre_comp)[16][3] = NULL;
1325 EC_POINT *generator = NULL;
1326 const EC_POINT *p = NULL;
1327 const BIGNUM *p_scalar = NULL;
1330 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1332 if (((x = BN_CTX_get(ctx)) == NULL) ||
1333 ((y = BN_CTX_get(ctx)) == NULL) ||
1334 ((z = BN_CTX_get(ctx)) == NULL) ||
1335 ((tmp_scalar = BN_CTX_get(ctx)) == NULL))
1340 pre = EC_EX_DATA_get_data(group->extra_data,
1341 nistp224_pre_comp_dup, nistp224_pre_comp_free,
1342 nistp224_pre_comp_clear_free);
1344 /* we have precomputation, try to use it */
1345 g_pre_comp = (const felem (*)[16][3]) pre->g_pre_comp;
1347 /* try to use the standard precomputation */
1348 g_pre_comp = &gmul[0];
1349 generator = EC_POINT_new(group);
1350 if (generator == NULL)
1352 /* get the generator from precomputation */
1353 if (!felem_to_BN(x, g_pre_comp[0][1][0]) ||
1354 !felem_to_BN(y, g_pre_comp[0][1][1]) ||
1355 !felem_to_BN(z, g_pre_comp[0][1][2]))
1357 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1360 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1361 generator, x, y, z, ctx))
1363 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1364 /* precomputation matches generator */
1367 /* we don't have valid precomputation:
1368 * treat the generator as a random point */
1369 num_points = num_points + 1;
1374 if (num_points >= 3)
1376 /* unless we precompute multiples for just one or two points,
1377 * converting those into affine form is time well spent */
1380 secrets = OPENSSL_malloc(num_points * sizeof(felem_bytearray));
1381 pre_comp = OPENSSL_malloc(num_points * 17 * 3 * sizeof(felem));
1383 tmp_felems = OPENSSL_malloc((num_points * 17 + 1) * sizeof(felem));
1384 if ((secrets == NULL) || (pre_comp == NULL) || (mixed && (tmp_felems == NULL)))
1386 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1390 /* we treat NULL scalars as 0, and NULL points as points at infinity,
1391 * i.e., they contribute nothing to the linear combination */
1392 memset(secrets, 0, num_points * sizeof(felem_bytearray));
1393 memset(pre_comp, 0, num_points * 17 * 3 * sizeof(felem));
1394 for (i = 0; i < num_points; ++i)
1399 p = EC_GROUP_get0_generator(group);
1403 /* the i^th point */
1406 p_scalar = scalars[i];
1408 if ((p_scalar != NULL) && (p != NULL))
1410 /* reduce scalar to 0 <= scalar < 2^224 */
1411 if ((BN_num_bits(p_scalar) > 224) || (BN_is_negative(p_scalar)))
1413 /* this is an unusual input, and we don't guarantee
1414 * constant-timeness */
1415 if (!BN_nnmod(tmp_scalar, p_scalar, &group->order, ctx))
1417 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1420 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1423 num_bytes = BN_bn2bin(p_scalar, tmp);
1424 flip_endian(secrets[i], tmp, num_bytes);
1425 /* precompute multiples */
1426 if ((!BN_to_felem(x_out, &p->X)) ||
1427 (!BN_to_felem(y_out, &p->Y)) ||
1428 (!BN_to_felem(z_out, &p->Z))) goto err;
1429 felem_assign(pre_comp[i][1][0], x_out);
1430 felem_assign(pre_comp[i][1][1], y_out);
1431 felem_assign(pre_comp[i][1][2], z_out);
1432 for (j = 2; j <= 16; ++j)
1437 pre_comp[i][j][0], pre_comp[i][j][1], pre_comp[i][j][2],
1438 pre_comp[i][1][0], pre_comp[i][1][1], pre_comp[i][1][2],
1439 0, pre_comp[i][j-1][0], pre_comp[i][j-1][1], pre_comp[i][j-1][2]);
1444 pre_comp[i][j][0], pre_comp[i][j][1], pre_comp[i][j][2],
1445 pre_comp[i][j/2][0], pre_comp[i][j/2][1], pre_comp[i][j/2][2]);
1451 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
1454 /* the scalar for the generator */
1455 if ((scalar != NULL) && (have_pre_comp))
1457 memset(g_secret, 0, sizeof g_secret);
1458 /* reduce scalar to 0 <= scalar < 2^224 */
1459 if ((BN_num_bits(scalar) > 224) || (BN_is_negative(scalar)))
1461 /* this is an unusual input, and we don't guarantee
1462 * constant-timeness */
1463 if (!BN_nnmod(tmp_scalar, scalar, &group->order, ctx))
1465 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1468 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1471 num_bytes = BN_bn2bin(scalar, tmp);
1472 flip_endian(g_secret, tmp, num_bytes);
1473 /* do the multiplication with generator precomputation*/
1474 batch_mul(x_out, y_out, z_out,
1475 (const felem_bytearray (*)) secrets, num_points,
1477 mixed, (const felem (*)[17][3]) pre_comp,
1481 /* do the multiplication without generator precomputation */
1482 batch_mul(x_out, y_out, z_out,
1483 (const felem_bytearray (*)) secrets, num_points,
1484 NULL, mixed, (const felem (*)[17][3]) pre_comp, NULL);
1485 /* reduce the output to its unique minimal representation */
1486 felem_contract(x_in, x_out);
1487 felem_contract(y_in, y_out);
1488 felem_contract(z_in, z_out);
1489 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
1490 (!felem_to_BN(z, z_in)))
1492 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1495 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
1499 if (generator != NULL)
1500 EC_POINT_free(generator);
1501 if (new_ctx != NULL)
1502 BN_CTX_free(new_ctx);
1503 if (secrets != NULL)
1504 OPENSSL_free(secrets);
1505 if (pre_comp != NULL)
1506 OPENSSL_free(pre_comp);
1507 if (tmp_felems != NULL)
1508 OPENSSL_free(tmp_felems);
1512 int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
1515 NISTP224_PRE_COMP *pre = NULL;
1517 BN_CTX *new_ctx = NULL;
1519 EC_POINT *generator = NULL;
1520 felem tmp_felems[32];
1522 /* throw away old precomputation */
1523 EC_EX_DATA_free_data(&group->extra_data, nistp224_pre_comp_dup,
1524 nistp224_pre_comp_free, nistp224_pre_comp_clear_free);
1526 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1528 if (((x = BN_CTX_get(ctx)) == NULL) ||
1529 ((y = BN_CTX_get(ctx)) == NULL))
1531 /* get the generator */
1532 if (group->generator == NULL) goto err;
1533 generator = EC_POINT_new(group);
1534 if (generator == NULL)
1536 BN_bin2bn(nistp224_curve_params[3], sizeof (felem_bytearray), x);
1537 BN_bin2bn(nistp224_curve_params[4], sizeof (felem_bytearray), y);
1538 if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
1540 if ((pre = nistp224_pre_comp_new()) == NULL)
1542 /* if the generator is the standard one, use built-in precomputation */
1543 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1545 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
1549 if ((!BN_to_felem(pre->g_pre_comp[0][1][0], &group->generator->X)) ||
1550 (!BN_to_felem(pre->g_pre_comp[0][1][1], &group->generator->Y)) ||
1551 (!BN_to_felem(pre->g_pre_comp[0][1][2], &group->generator->Z)))
1553 /* compute 2^56*G, 2^112*G, 2^168*G for the first table,
1554 * 2^28*G, 2^84*G, 2^140*G, 2^196*G for the second one
1556 for (i = 1; i <= 8; i <<= 1)
1559 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2],
1560 pre->g_pre_comp[0][i][0], pre->g_pre_comp[0][i][1], pre->g_pre_comp[0][i][2]);
1561 for (j = 0; j < 27; ++j)
1564 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2],
1565 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1570 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2],
1571 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1572 for (j = 0; j < 27; ++j)
1575 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2],
1576 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2]);
1579 for (i = 0; i < 2; i++)
1581 /* g_pre_comp[i][0] is the point at infinity */
1582 memset(pre->g_pre_comp[i][0], 0, sizeof(pre->g_pre_comp[i][0]));
1583 /* the remaining multiples */
1584 /* 2^56*G + 2^112*G resp. 2^84*G + 2^140*G */
1586 pre->g_pre_comp[i][6][0], pre->g_pre_comp[i][6][1],
1587 pre->g_pre_comp[i][6][2], pre->g_pre_comp[i][4][0],
1588 pre->g_pre_comp[i][4][1], pre->g_pre_comp[i][4][2],
1589 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1590 pre->g_pre_comp[i][2][2]);
1591 /* 2^56*G + 2^168*G resp. 2^84*G + 2^196*G */
1593 pre->g_pre_comp[i][10][0], pre->g_pre_comp[i][10][1],
1594 pre->g_pre_comp[i][10][2], pre->g_pre_comp[i][8][0],
1595 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1596 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1597 pre->g_pre_comp[i][2][2]);
1598 /* 2^112*G + 2^168*G resp. 2^140*G + 2^196*G */
1600 pre->g_pre_comp[i][12][0], pre->g_pre_comp[i][12][1],
1601 pre->g_pre_comp[i][12][2], pre->g_pre_comp[i][8][0],
1602 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][2],
1603 0, pre->g_pre_comp[i][4][0], pre->g_pre_comp[i][4][1],
1604 pre->g_pre_comp[i][4][2]);
1605 /* 2^56*G + 2^112*G + 2^168*G resp. 2^84*G + 2^140*G + 2^196*G */
1607 pre->g_pre_comp[i][14][0], pre->g_pre_comp[i][14][1],
1608 pre->g_pre_comp[i][14][2], pre->g_pre_comp[i][12][0],
1609 pre->g_pre_comp[i][12][1], pre->g_pre_comp[i][12][2],
1610 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1611 pre->g_pre_comp[i][2][2]);
1612 for (j = 1; j < 8; ++j)
1614 /* odd multiples: add G resp. 2^28*G */
1616 pre->g_pre_comp[i][2*j+1][0], pre->g_pre_comp[i][2*j+1][1],
1617 pre->g_pre_comp[i][2*j+1][2], pre->g_pre_comp[i][2*j][0],
1618 pre->g_pre_comp[i][2*j][1], pre->g_pre_comp[i][2*j][2],
1619 0, pre->g_pre_comp[i][1][0], pre->g_pre_comp[i][1][1],
1620 pre->g_pre_comp[i][1][2]);
1623 make_points_affine(31, &(pre->g_pre_comp[0][1]), tmp_felems);
1625 if (!EC_EX_DATA_set_data(&group->extra_data, pre, nistp224_pre_comp_dup,
1626 nistp224_pre_comp_free, nistp224_pre_comp_clear_free))
1632 if (generator != NULL)
1633 EC_POINT_free(generator);
1634 if (new_ctx != NULL)
1635 BN_CTX_free(new_ctx);
1637 nistp224_pre_comp_free(pre);
1641 int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
1643 if (EC_EX_DATA_get_data(group->extra_data, nistp224_pre_comp_dup,
1644 nistp224_pre_comp_free, nistp224_pre_comp_clear_free)
1652 static void *dummy=&dummy;