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 /******************************************************************************/
346 * Field operations, using the internal representation of field elements.
347 * NB! These operations are specific to our point multiplication and cannot be
348 * expected to be correct in general - e.g., multiplication with a large scalar
349 * will cause an overflow.
353 static void felem_one(felem out)
361 static void felem_assign(felem out, const felem in)
369 /* Sum two field elements: out += in */
370 static void felem_sum(felem out, const felem in)
378 /* Get negative value: out = -in */
379 /* Assumes in[i] < 2^57 */
380 static void felem_neg(felem out, const felem in)
382 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
383 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
384 static const limb two58m42m2 = (((limb) 1) << 58) -
385 (((limb) 1) << 42) - (((limb) 1) << 2);
387 /* Set to 0 mod 2^224-2^96+1 to ensure out > in */
388 out[0] = two58p2 - in[0];
389 out[1] = two58m42m2 - in[1];
390 out[2] = two58m2 - in[2];
391 out[3] = two58m2 - in[3];
394 /* Subtract field elements: out -= in */
395 /* Assumes in[i] < 2^57 */
396 static void felem_diff(felem out, const felem in)
398 static const limb two58p2 = (((limb) 1) << 58) + (((limb) 1) << 2);
399 static const limb two58m2 = (((limb) 1) << 58) - (((limb) 1) << 2);
400 static const limb two58m42m2 = (((limb) 1) << 58) -
401 (((limb) 1) << 42) - (((limb) 1) << 2);
403 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
405 out[1] += two58m42m2;
415 /* Subtract in unreduced 128-bit mode: out -= in */
416 /* Assumes in[i] < 2^119 */
417 static void widefelem_diff(widefelem out, const widefelem in)
419 static const widelimb two120 = ((widelimb) 1) << 120;
420 static const widelimb two120m64 = (((widelimb) 1) << 120) -
421 (((widelimb) 1) << 64);
422 static const widelimb two120m104m64 = (((widelimb) 1) << 120) -
423 (((widelimb) 1) << 104) - (((widelimb) 1) << 64);
425 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
430 out[4] += two120m104m64;
443 /* Subtract in mixed mode: out128 -= in64 */
445 static void felem_diff_128_64(widefelem out, const felem in)
447 static const widelimb two64p8 = (((widelimb) 1) << 64) +
448 (((widelimb) 1) << 8);
449 static const widelimb two64m8 = (((widelimb) 1) << 64) -
450 (((widelimb) 1) << 8);
451 static const widelimb two64m48m8 = (((widelimb) 1) << 64) -
452 (((widelimb) 1) << 48) - (((widelimb) 1) << 8);
454 /* Add 0 mod 2^224-2^96+1 to ensure out > in */
456 out[1] += two64m48m8;
466 /* Multiply a field element by a scalar: out = out * scalar
467 * The scalars we actually use are small, so results fit without overflow */
468 static void felem_scalar(felem out, const limb scalar)
476 /* Multiply an unreduced field element by a scalar: out = out * scalar
477 * The scalars we actually use are small, so results fit without overflow */
478 static void widefelem_scalar(widefelem out, const widelimb scalar)
489 /* Square a field element: out = in^2 */
490 static void felem_square(widefelem out, const felem in)
492 limb tmp0, tmp1, tmp2;
493 tmp0 = 2 * in[0]; tmp1 = 2 * in[1]; tmp2 = 2 * in[2];
494 out[0] = ((widelimb) in[0]) * in[0];
495 out[1] = ((widelimb) in[0]) * tmp1;
496 out[2] = ((widelimb) in[0]) * tmp2 + ((widelimb) in[1]) * in[1];
497 out[3] = ((widelimb) in[3]) * tmp0 +
498 ((widelimb) in[1]) * tmp2;
499 out[4] = ((widelimb) in[3]) * tmp1 + ((widelimb) in[2]) * in[2];
500 out[5] = ((widelimb) in[3]) * tmp2;
501 out[6] = ((widelimb) in[3]) * in[3];
504 /* Multiply two field elements: out = in1 * in2 */
505 static void felem_mul(widefelem out, const felem in1, const felem in2)
507 out[0] = ((widelimb) in1[0]) * in2[0];
508 out[1] = ((widelimb) in1[0]) * in2[1] + ((widelimb) in1[1]) * in2[0];
509 out[2] = ((widelimb) in1[0]) * in2[2] + ((widelimb) in1[1]) * in2[1] +
510 ((widelimb) in1[2]) * in2[0];
511 out[3] = ((widelimb) in1[0]) * in2[3] + ((widelimb) in1[1]) * in2[2] +
512 ((widelimb) in1[2]) * in2[1] + ((widelimb) in1[3]) * in2[0];
513 out[4] = ((widelimb) in1[1]) * in2[3] + ((widelimb) in1[2]) * in2[2] +
514 ((widelimb) in1[3]) * in2[1];
515 out[5] = ((widelimb) in1[2]) * in2[3] + ((widelimb) in1[3]) * in2[2];
516 out[6] = ((widelimb) in1[3]) * in2[3];
519 /* Reduce seven 128-bit coefficients to four 64-bit coefficients.
520 * Requires in[i] < 2^126,
521 * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^56, out[3] <= 2^56 + 2^16 */
522 static void felem_reduce(felem out, const widefelem in)
524 static const widelimb two127p15 = (((widelimb) 1) << 127) +
525 (((widelimb) 1) << 15);
526 static const widelimb two127m71 = (((widelimb) 1) << 127) -
527 (((widelimb) 1) << 71);
528 static const widelimb two127m71m55 = (((widelimb) 1) << 127) -
529 (((widelimb) 1) << 71) - (((widelimb) 1) << 55);
532 /* Add 0 mod 2^224-2^96+1 to ensure all differences are positive */
533 output[0] = in[0] + two127p15;
534 output[1] = in[1] + two127m71m55;
535 output[2] = in[2] + two127m71;
539 /* Eliminate in[4], in[5], in[6] */
540 output[4] += in[6] >> 16;
541 output[3] += (in[6] & 0xffff) << 40;
544 output[3] += in[5] >> 16;
545 output[2] += (in[5] & 0xffff) << 40;
548 output[2] += output[4] >> 16;
549 output[1] += (output[4] & 0xffff) << 40;
550 output[0] -= output[4];
552 /* Carry 2 -> 3 -> 4 */
553 output[3] += output[2] >> 56;
554 output[2] &= 0x00ffffffffffffff;
556 output[4] = output[3] >> 56;
557 output[3] &= 0x00ffffffffffffff;
559 /* Now output[2] < 2^56, output[3] < 2^56, output[4] < 2^72 */
561 /* Eliminate output[4] */
562 output[2] += output[4] >> 16;
563 /* output[2] < 2^56 + 2^56 = 2^57 */
564 output[1] += (output[4] & 0xffff) << 40;
565 output[0] -= output[4];
567 /* Carry 0 -> 1 -> 2 -> 3 */
568 output[1] += output[0] >> 56;
569 out[0] = output[0] & 0x00ffffffffffffff;
571 output[2] += output[1] >> 56;
572 /* output[2] < 2^57 + 2^72 */
573 out[1] = output[1] & 0x00ffffffffffffff;
574 output[3] += output[2] >> 56;
575 /* output[3] <= 2^56 + 2^16 */
576 out[2] = output[2] & 0x00ffffffffffffff;
579 * out[0] < 2^56, out[1] < 2^56, out[2] < 2^56,
580 * out[3] <= 2^56 + 2^16 (due to final carry),
586 static void felem_square_reduce(felem out, const felem in)
589 felem_square(tmp, in);
590 felem_reduce(out, tmp);
593 static void felem_mul_reduce(felem out, const felem in1, const felem in2)
596 felem_mul(tmp, in1, in2);
597 felem_reduce(out, tmp);
600 /* Reduce to unique minimal representation.
601 * Requires 0 <= in < 2*p (always call felem_reduce first) */
602 static void felem_contract(felem out, const felem in)
604 static const int64_t two56 = ((limb) 1) << 56;
605 /* 0 <= in < 2*p, p = 2^224 - 2^96 + 1 */
606 /* if in > p , reduce in = in - 2^224 + 2^96 - 1 */
612 /* Case 1: a = 1 iff in >= 2^224 */
616 tmp[3] &= 0x00ffffffffffffff;
617 /* Case 2: a = 0 iff p <= in < 2^224, i.e.,
618 * the high 128 bits are all 1 and the lower part is non-zero */
619 a = ((in[3] & in[2] & (in[1] | 0x000000ffffffffff)) + 1) |
620 (((int64_t)(in[0] + (in[1] & 0x000000ffffffffff)) - 1) >> 63);
621 a &= 0x00ffffffffffffff;
622 /* turn a into an all-one mask (if a = 0) or an all-zero mask */
624 /* subtract 2^224 - 2^96 + 1 if a is all-one*/
625 tmp[3] &= a ^ 0xffffffffffffffff;
626 tmp[2] &= a ^ 0xffffffffffffffff;
627 tmp[1] &= (a ^ 0xffffffffffffffff) | 0x000000ffffffffff;
630 /* eliminate negative coefficients: if tmp[0] is negative, tmp[1] must
631 * be non-zero, so we only need one step */
636 /* carry 1 -> 2 -> 3 */
637 tmp[2] += tmp[1] >> 56;
638 tmp[1] &= 0x00ffffffffffffff;
640 tmp[3] += tmp[2] >> 56;
641 tmp[2] &= 0x00ffffffffffffff;
643 /* Now 0 <= out < p */
650 /* Zero-check: returns 1 if input is 0, and 0 otherwise.
651 * We know that field elements are reduced to in < 2^225,
652 * so we only need to check three cases: 0, 2^224 - 2^96 + 1,
653 * and 2^225 - 2^97 + 2 */
654 static limb felem_is_zero(const felem in)
656 limb zero, two224m96p1, two225m97p2;
658 zero = in[0] | in[1] | in[2] | in[3];
659 zero = (((int64_t)(zero) - 1) >> 63) & 1;
660 two224m96p1 = (in[0] ^ 1) | (in[1] ^ 0x00ffff0000000000)
661 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x00ffffffffffffff);
662 two224m96p1 = (((int64_t)(two224m96p1) - 1) >> 63) & 1;
663 two225m97p2 = (in[0] ^ 2) | (in[1] ^ 0x00fffe0000000000)
664 | (in[2] ^ 0x00ffffffffffffff) | (in[3] ^ 0x01ffffffffffffff);
665 two225m97p2 = (((int64_t)(two225m97p2) - 1) >> 63) & 1;
666 return (zero | two224m96p1 | two225m97p2);
669 static limb felem_is_zero_int(const felem in)
671 return (int) (felem_is_zero(in) & ((limb)1));
674 /* Invert a field element */
675 /* Computation chain copied from djb's code */
676 static void felem_inv(felem out, const felem in)
678 felem ftmp, ftmp2, ftmp3, ftmp4;
682 felem_square(tmp, in); felem_reduce(ftmp, tmp); /* 2 */
683 felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp); /* 2^2 - 1 */
684 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp); /* 2^3 - 2 */
685 felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp); /* 2^3 - 1 */
686 felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp); /* 2^4 - 2 */
687 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp); /* 2^5 - 4 */
688 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp); /* 2^6 - 8 */
689 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp); /* 2^6 - 1 */
690 felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp); /* 2^7 - 2 */
691 for (i = 0; i < 5; ++i) /* 2^12 - 2^6 */
693 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
695 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp2, tmp); /* 2^12 - 1 */
696 felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp); /* 2^13 - 2 */
697 for (i = 0; i < 11; ++i) /* 2^24 - 2^12 */
699 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
701 felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp2, tmp); /* 2^24 - 1 */
702 felem_square(tmp, ftmp2); felem_reduce(ftmp3, tmp); /* 2^25 - 2 */
703 for (i = 0; i < 23; ++i) /* 2^48 - 2^24 */
705 felem_square(tmp, ftmp3); felem_reduce(ftmp3, tmp);
707 felem_mul(tmp, ftmp3, ftmp2); felem_reduce(ftmp3, tmp); /* 2^48 - 1 */
708 felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp); /* 2^49 - 2 */
709 for (i = 0; i < 47; ++i) /* 2^96 - 2^48 */
711 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
713 felem_mul(tmp, ftmp3, ftmp4); felem_reduce(ftmp3, tmp); /* 2^96 - 1 */
714 felem_square(tmp, ftmp3); felem_reduce(ftmp4, tmp); /* 2^97 - 2 */
715 for (i = 0; i < 23; ++i) /* 2^120 - 2^24 */
717 felem_square(tmp, ftmp4); felem_reduce(ftmp4, tmp);
719 felem_mul(tmp, ftmp2, ftmp4); felem_reduce(ftmp2, tmp); /* 2^120 - 1 */
720 for (i = 0; i < 6; ++i) /* 2^126 - 2^6 */
722 felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
724 felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp); /* 2^126 - 1 */
725 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp); /* 2^127 - 2 */
726 felem_mul(tmp, ftmp, in); felem_reduce(ftmp, tmp); /* 2^127 - 1 */
727 for (i = 0; i < 97; ++i) /* 2^224 - 2^97 */
729 felem_square(tmp, ftmp); felem_reduce(ftmp, tmp);
731 felem_mul(tmp, ftmp, ftmp3); felem_reduce(out, tmp); /* 2^224 - 2^96 - 1 */
734 /* Copy in constant time:
735 * if icopy == 1, copy in to out,
736 * if icopy == 0, copy out to itself. */
738 copy_conditional(felem out, const felem in, limb icopy)
741 /* icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one */
742 const limb copy = -icopy;
743 for (i = 0; i < 4; ++i)
745 const limb tmp = copy & (in[i] ^ out[i]);
750 /******************************************************************************/
751 /* ELLIPTIC CURVE POINT OPERATIONS
753 * Points are represented in Jacobian projective coordinates:
754 * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
755 * or to the point at infinity if Z == 0.
760 * Double an elliptic curve point:
761 * (X', Y', Z') = 2 * (X, Y, Z), where
762 * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
763 * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^2
764 * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
765 * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
766 * while x_out == y_in is not (maybe this works, but it's not tested).
769 point_double(felem x_out, felem y_out, felem z_out,
770 const felem x_in, const felem y_in, const felem z_in)
773 felem delta, gamma, beta, alpha, ftmp, ftmp2;
775 felem_assign(ftmp, x_in);
776 felem_assign(ftmp2, x_in);
779 felem_square(tmp, z_in);
780 felem_reduce(delta, tmp);
783 felem_square(tmp, y_in);
784 felem_reduce(gamma, tmp);
787 felem_mul(tmp, x_in, gamma);
788 felem_reduce(beta, tmp);
790 /* alpha = 3*(x-delta)*(x+delta) */
791 felem_diff(ftmp, delta);
792 /* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
793 felem_sum(ftmp2, delta);
794 /* ftmp2[i] < 2^57 + 2^57 = 2^58 */
795 felem_scalar(ftmp2, 3);
796 /* ftmp2[i] < 3 * 2^58 < 2^60 */
797 felem_mul(tmp, ftmp, ftmp2);
798 /* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
799 felem_reduce(alpha, tmp);
801 /* x' = alpha^2 - 8*beta */
802 felem_square(tmp, alpha);
803 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
804 felem_assign(ftmp, beta);
805 felem_scalar(ftmp, 8);
806 /* ftmp[i] < 8 * 2^57 = 2^60 */
807 felem_diff_128_64(tmp, ftmp);
808 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
809 felem_reduce(x_out, tmp);
811 /* z' = (y + z)^2 - gamma - delta */
812 felem_sum(delta, gamma);
813 /* delta[i] < 2^57 + 2^57 = 2^58 */
814 felem_assign(ftmp, y_in);
815 felem_sum(ftmp, z_in);
816 /* ftmp[i] < 2^57 + 2^57 = 2^58 */
817 felem_square(tmp, ftmp);
818 /* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
819 felem_diff_128_64(tmp, delta);
820 /* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
821 felem_reduce(z_out, tmp);
823 /* y' = alpha*(4*beta - x') - 8*gamma^2 */
824 felem_scalar(beta, 4);
825 /* beta[i] < 4 * 2^57 = 2^59 */
826 felem_diff(beta, x_out);
827 /* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
828 felem_mul(tmp, alpha, beta);
829 /* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
830 felem_square(tmp2, gamma);
831 /* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
832 widefelem_scalar(tmp2, 8);
833 /* tmp2[i] < 8 * 2^116 = 2^119 */
834 widefelem_diff(tmp, tmp2);
835 /* tmp[i] < 2^119 + 2^120 < 2^121 */
836 felem_reduce(y_out, tmp);
840 * Add two elliptic curve points:
841 * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
842 * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
843 * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
844 * 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) -
845 * Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
846 * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2)
848 * This runs faster if 'mixed' is set, which requires Z_2 = 1 or Z_2 = 0.
851 /* This function is not entirely constant-time:
852 * it includes a branch for checking whether the two input points are equal,
853 * (while not equal to the point at infinity).
854 * This case never happens during single point multiplication,
855 * so there is no timing leak for ECDH or ECDSA signing. */
856 static void point_add(felem x3, felem y3, felem z3,
857 const felem x1, const felem y1, const felem z1,
858 const int mixed, const felem x2, const felem y2, const felem z2)
860 felem ftmp, ftmp2, ftmp3, ftmp4, ftmp5, x_out, y_out, z_out;
862 limb z1_is_zero, z2_is_zero, x_equal, y_equal;
867 felem_square(tmp, z2);
868 felem_reduce(ftmp2, tmp);
871 felem_mul(tmp, ftmp2, z2);
872 felem_reduce(ftmp4, tmp);
874 /* ftmp4 = z2^3*y1 */
875 felem_mul(tmp2, ftmp4, y1);
876 felem_reduce(ftmp4, tmp2);
878 /* ftmp2 = z2^2*x1 */
879 felem_mul(tmp2, ftmp2, x1);
880 felem_reduce(ftmp2, tmp2);
884 /* We'll assume z2 = 1 (special case z2 = 0 is handled later) */
886 /* ftmp4 = z2^3*y1 */
887 felem_assign(ftmp4, y1);
889 /* ftmp2 = z2^2*x1 */
890 felem_assign(ftmp2, x1);
894 felem_square(tmp, z1);
895 felem_reduce(ftmp, tmp);
898 felem_mul(tmp, ftmp, z1);
899 felem_reduce(ftmp3, tmp);
902 felem_mul(tmp, ftmp3, y2);
903 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
905 /* ftmp3 = z1^3*y2 - z2^3*y1 */
906 felem_diff_128_64(tmp, ftmp4);
907 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
908 felem_reduce(ftmp3, tmp);
911 felem_mul(tmp, ftmp, x2);
912 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
914 /* ftmp = z1^2*x2 - z2^2*x1 */
915 felem_diff_128_64(tmp, ftmp2);
916 /* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
917 felem_reduce(ftmp, tmp);
919 /* the formulae are incorrect if the points are equal
920 * so we check for this and do doubling if this happens */
921 x_equal = felem_is_zero(ftmp);
922 y_equal = felem_is_zero(ftmp3);
923 z1_is_zero = felem_is_zero(z1);
924 z2_is_zero = felem_is_zero(z2);
925 /* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
926 if (x_equal && y_equal && !z1_is_zero && !z2_is_zero)
928 point_double(x3, y3, z3, x1, y1, z1);
935 felem_mul(tmp, z1, z2);
936 felem_reduce(ftmp5, tmp);
940 /* special case z2 = 0 is handled later */
941 felem_assign(ftmp5, z1);
944 /* z_out = (z1^2*x2 - z2^2*x1)*(z1*z2) */
945 felem_mul(tmp, ftmp, ftmp5);
946 felem_reduce(z_out, tmp);
948 /* ftmp = (z1^2*x2 - z2^2*x1)^2 */
949 felem_assign(ftmp5, ftmp);
950 felem_square(tmp, ftmp);
951 felem_reduce(ftmp, tmp);
953 /* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
954 felem_mul(tmp, ftmp, ftmp5);
955 felem_reduce(ftmp5, tmp);
957 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
958 felem_mul(tmp, ftmp2, ftmp);
959 felem_reduce(ftmp2, tmp);
961 /* tmp = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
962 felem_mul(tmp, ftmp4, ftmp5);
963 /* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
965 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
966 felem_square(tmp2, ftmp3);
967 /* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
969 /* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
970 felem_diff_128_64(tmp2, ftmp5);
971 /* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
973 /* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
974 felem_assign(ftmp5, ftmp2);
975 felem_scalar(ftmp5, 2);
976 /* ftmp5[i] < 2 * 2^57 = 2^58 */
979 * x_out = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
980 * 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2
982 felem_diff_128_64(tmp2, ftmp5);
983 /* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
984 felem_reduce(x_out, tmp2);
986 /* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out */
987 felem_diff(ftmp2, x_out);
988 /* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
990 /* tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) */
991 felem_mul(tmp2, ftmp3, ftmp2);
992 /* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
995 * y_out = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - x_out) -
996 * z2^3*y1*(z1^2*x2 - z2^2*x1)^3
998 widefelem_diff(tmp2, tmp);
999 /* tmp2[i] < 2^118 + 2^120 < 2^121 */
1000 felem_reduce(y_out, tmp2);
1002 /* the result (x_out, y_out, z_out) is incorrect if one of the inputs is
1003 * the point at infinity, so we need to check for this separately */
1005 /* if point 1 is at infinity, copy point 2 to output, and vice versa */
1006 copy_conditional(x_out, x2, z1_is_zero);
1007 copy_conditional(x_out, x1, z2_is_zero);
1008 copy_conditional(y_out, y2, z1_is_zero);
1009 copy_conditional(y_out, y1, z2_is_zero);
1010 copy_conditional(z_out, z2, z1_is_zero);
1011 copy_conditional(z_out, z1, z2_is_zero);
1012 felem_assign(x3, x_out);
1013 felem_assign(y3, y_out);
1014 felem_assign(z3, z_out);
1017 /* select_point selects the |idx|th point from a precomputation table and
1018 * copies it to out. */
1019 static void select_point(const u64 idx, unsigned int size, const felem pre_comp[/*size*/][3], felem out[3])
1022 limb *outlimbs = &out[0][0];
1023 memset(outlimbs, 0, 3 * sizeof(felem));
1025 for (i = 0; i < size; i++)
1027 const limb *inlimbs = &pre_comp[i][0][0];
1034 for (j = 0; j < 4 * 3; j++)
1035 outlimbs[j] |= inlimbs[j] & mask;
1039 /* get_bit returns the |i|th bit in |in| */
1040 static char get_bit(const felem_bytearray in, unsigned i)
1044 return (in[i >> 3] >> (i & 7)) & 1;
1047 /* Interleaved point multiplication using precomputed point multiples:
1048 * The small point multiples 0*P, 1*P, ..., 16*P are in pre_comp[],
1049 * the scalars in scalars[]. If g_scalar is non-NULL, we also add this multiple
1050 * of the generator, using certain (large) precomputed multiples in g_pre_comp.
1051 * Output point (X, Y, Z) is stored in x_out, y_out, z_out */
1052 static void batch_mul(felem x_out, felem y_out, felem z_out,
1053 const felem_bytearray scalars[], const unsigned num_points, const u8 *g_scalar,
1054 const int mixed, const felem pre_comp[][17][3], const felem g_pre_comp[2][16][3])
1058 unsigned gen_mul = (g_scalar != NULL);
1059 felem nq[3], tmp[4];
1063 /* set nq to the point at infinity */
1064 memset(nq, 0, 3 * sizeof(felem));
1066 /* Loop over all scalars msb-to-lsb, interleaving additions
1067 * of multiples of the generator (two in each of the last 28 rounds)
1068 * and additions of other points multiples (every 5th round).
1070 skip = 1; /* save two point operations in the first round */
1071 for (i = (num_points ? 220 : 27); i >= 0; --i)
1075 point_double(nq[0], nq[1], nq[2], nq[0], nq[1], nq[2]);
1077 /* add multiples of the generator */
1078 if (gen_mul && (i <= 27))
1080 /* first, look 28 bits upwards */
1081 bits = get_bit(g_scalar, i + 196) << 3;
1082 bits |= get_bit(g_scalar, i + 140) << 2;
1083 bits |= get_bit(g_scalar, i + 84) << 1;
1084 bits |= get_bit(g_scalar, i + 28);
1085 /* select the point to add, in constant time */
1086 select_point(bits, 16, g_pre_comp[1], tmp);
1090 point_add(nq[0], nq[1], nq[2],
1091 nq[0], nq[1], nq[2],
1092 1 /* mixed */, tmp[0], tmp[1], tmp[2]);
1096 memcpy(nq, tmp, 3 * sizeof(felem));
1100 /* second, look at the current position */
1101 bits = get_bit(g_scalar, i + 168) << 3;
1102 bits |= get_bit(g_scalar, i + 112) << 2;
1103 bits |= get_bit(g_scalar, i + 56) << 1;
1104 bits |= get_bit(g_scalar, i);
1105 /* select the point to add, in constant time */
1106 select_point(bits, 16, g_pre_comp[0], tmp);
1107 point_add(nq[0], nq[1], nq[2],
1108 nq[0], nq[1], nq[2],
1109 1 /* mixed */, tmp[0], tmp[1], tmp[2]);
1112 /* do other additions every 5 doublings */
1113 if (num_points && (i % 5 == 0))
1115 /* loop over all scalars */
1116 for (num = 0; num < num_points; ++num)
1118 bits = get_bit(scalars[num], i + 4) << 5;
1119 bits |= get_bit(scalars[num], i + 3) << 4;
1120 bits |= get_bit(scalars[num], i + 2) << 3;
1121 bits |= get_bit(scalars[num], i + 1) << 2;
1122 bits |= get_bit(scalars[num], i) << 1;
1123 bits |= get_bit(scalars[num], i - 1);
1124 ec_GFp_nistp_recode_scalar_bits(&sign, &digit, bits);
1126 /* select the point to add or subtract */
1127 select_point(digit, 17, pre_comp[num], tmp);
1128 felem_neg(tmp[3], tmp[1]); /* (X, -Y, Z) is the negative point */
1129 copy_conditional(tmp[1], tmp[3], sign);
1133 point_add(nq[0], nq[1], nq[2],
1134 nq[0], nq[1], nq[2],
1135 mixed, tmp[0], tmp[1], tmp[2]);
1139 memcpy(nq, tmp, 3 * sizeof(felem));
1145 felem_assign(x_out, nq[0]);
1146 felem_assign(y_out, nq[1]);
1147 felem_assign(z_out, nq[2]);
1150 /******************************************************************************/
1151 /* FUNCTIONS TO MANAGE PRECOMPUTATION
1154 static NISTP224_PRE_COMP *nistp224_pre_comp_new()
1156 NISTP224_PRE_COMP *ret = NULL;
1157 ret = (NISTP224_PRE_COMP *) OPENSSL_malloc(sizeof *ret);
1160 ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
1163 memset(ret->g_pre_comp, 0, sizeof(ret->g_pre_comp));
1164 ret->references = 1;
1168 static void *nistp224_pre_comp_dup(void *src_)
1170 NISTP224_PRE_COMP *src = src_;
1172 /* no need to actually copy, these objects never change! */
1173 CRYPTO_add(&src->references, 1, CRYPTO_LOCK_EC_PRE_COMP);
1178 static void nistp224_pre_comp_free(void *pre_)
1181 NISTP224_PRE_COMP *pre = pre_;
1186 i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1193 static void nistp224_pre_comp_clear_free(void *pre_)
1196 NISTP224_PRE_COMP *pre = pre_;
1201 i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
1205 OPENSSL_cleanse(pre, sizeof *pre);
1209 /******************************************************************************/
1210 /* OPENSSL EC_METHOD FUNCTIONS
1213 int ec_GFp_nistp224_group_init(EC_GROUP *group)
1216 ret = ec_GFp_simple_group_init(group);
1217 group->a_is_minus3 = 1;
1221 int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
1222 const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
1225 BN_CTX *new_ctx = NULL;
1226 BIGNUM *curve_p, *curve_a, *curve_b;
1229 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1231 if (((curve_p = BN_CTX_get(ctx)) == NULL) ||
1232 ((curve_a = BN_CTX_get(ctx)) == NULL) ||
1233 ((curve_b = BN_CTX_get(ctx)) == NULL)) goto err;
1234 BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
1235 BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
1236 BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
1237 if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) ||
1238 (BN_cmp(curve_b, b)))
1240 ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
1241 EC_R_WRONG_CURVE_PARAMETERS);
1244 group->field_mod_func = BN_nist_mod_224;
1245 ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
1248 if (new_ctx != NULL)
1249 BN_CTX_free(new_ctx);
1253 /* Takes the Jacobian coordinates (X, Y, Z) of a point and returns
1254 * (X', Y') = (X/Z^2, Y/Z^3) */
1255 int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
1256 const EC_POINT *point, BIGNUM *x, BIGNUM *y, BN_CTX *ctx)
1258 felem z1, z2, x_in, y_in, x_out, y_out;
1261 if (EC_POINT_is_at_infinity(group, point))
1263 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1264 EC_R_POINT_AT_INFINITY);
1267 if ((!BN_to_felem(x_in, &point->X)) || (!BN_to_felem(y_in, &point->Y)) ||
1268 (!BN_to_felem(z1, &point->Z))) return 0;
1270 felem_square(tmp, z2); felem_reduce(z1, tmp);
1271 felem_mul(tmp, x_in, z1); felem_reduce(x_in, tmp);
1272 felem_contract(x_out, x_in);
1275 if (!felem_to_BN(x, x_out)) {
1276 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1281 felem_mul(tmp, z1, z2); felem_reduce(z1, tmp);
1282 felem_mul(tmp, y_in, z1); felem_reduce(y_in, tmp);
1283 felem_contract(y_out, y_in);
1286 if (!felem_to_BN(y, y_out)) {
1287 ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
1295 static void make_points_affine(size_t num, felem points[/*num*/][3], felem tmp_felems[/*num+1*/])
1297 /* Runs in constant time, unless an input is the point at infinity
1298 * (which normally shouldn't happen). */
1299 ec_GFp_nistp_points_make_affine_internal(
1304 (void (*)(void *)) felem_one,
1305 (int (*)(const void *)) felem_is_zero_int,
1306 (void (*)(void *, const void *)) felem_assign,
1307 (void (*)(void *, const void *)) felem_square_reduce,
1308 (void (*)(void *, const void *, const void *)) felem_mul_reduce,
1309 (void (*)(void *, const void *)) felem_inv,
1310 (void (*)(void *, const void *)) felem_contract);
1313 /* Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL values
1314 * Result is stored in r (r can equal one of the inputs). */
1315 int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
1316 const BIGNUM *scalar, size_t num, const EC_POINT *points[],
1317 const BIGNUM *scalars[], BN_CTX *ctx)
1323 BN_CTX *new_ctx = NULL;
1324 BIGNUM *x, *y, *z, *tmp_scalar;
1325 felem_bytearray g_secret;
1326 felem_bytearray *secrets = NULL;
1327 felem (*pre_comp)[17][3] = NULL;
1328 felem *tmp_felems = NULL;
1329 felem_bytearray tmp;
1331 int have_pre_comp = 0;
1332 size_t num_points = num;
1333 felem x_in, y_in, z_in, x_out, y_out, z_out;
1334 NISTP224_PRE_COMP *pre = NULL;
1335 const felem (*g_pre_comp)[16][3] = NULL;
1336 EC_POINT *generator = NULL;
1337 const EC_POINT *p = NULL;
1338 const BIGNUM *p_scalar = NULL;
1341 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1343 if (((x = BN_CTX_get(ctx)) == NULL) ||
1344 ((y = BN_CTX_get(ctx)) == NULL) ||
1345 ((z = BN_CTX_get(ctx)) == NULL) ||
1346 ((tmp_scalar = BN_CTX_get(ctx)) == NULL))
1351 pre = EC_EX_DATA_get_data(group->extra_data,
1352 nistp224_pre_comp_dup, nistp224_pre_comp_free,
1353 nistp224_pre_comp_clear_free);
1355 /* we have precomputation, try to use it */
1356 g_pre_comp = (const felem (*)[16][3]) pre->g_pre_comp;
1358 /* try to use the standard precomputation */
1359 g_pre_comp = &gmul[0];
1360 generator = EC_POINT_new(group);
1361 if (generator == NULL)
1363 /* get the generator from precomputation */
1364 if (!felem_to_BN(x, g_pre_comp[0][1][0]) ||
1365 !felem_to_BN(y, g_pre_comp[0][1][1]) ||
1366 !felem_to_BN(z, g_pre_comp[0][1][2]))
1368 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1371 if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
1372 generator, x, y, z, ctx))
1374 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1375 /* precomputation matches generator */
1378 /* we don't have valid precomputation:
1379 * treat the generator as a random point */
1380 num_points = num_points + 1;
1385 if (num_points >= 3)
1387 /* unless we precompute multiples for just one or two points,
1388 * converting those into affine form is time well spent */
1391 secrets = OPENSSL_malloc(num_points * sizeof(felem_bytearray));
1392 pre_comp = OPENSSL_malloc(num_points * 17 * 3 * sizeof(felem));
1394 tmp_felems = OPENSSL_malloc((num_points * 17 + 1) * sizeof(felem));
1395 if ((secrets == NULL) || (pre_comp == NULL) || (mixed && (tmp_felems == NULL)))
1397 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
1401 /* we treat NULL scalars as 0, and NULL points as points at infinity,
1402 * i.e., they contribute nothing to the linear combination */
1403 memset(secrets, 0, num_points * sizeof(felem_bytearray));
1404 memset(pre_comp, 0, num_points * 17 * 3 * sizeof(felem));
1405 for (i = 0; i < num_points; ++i)
1410 p = EC_GROUP_get0_generator(group);
1414 /* the i^th point */
1417 p_scalar = scalars[i];
1419 if ((p_scalar != NULL) && (p != NULL))
1421 /* reduce scalar to 0 <= scalar < 2^224 */
1422 if ((BN_num_bits(p_scalar) > 224) || (BN_is_negative(p_scalar)))
1424 /* this is an unusual input, and we don't guarantee
1425 * constant-timeness */
1426 if (!BN_nnmod(tmp_scalar, p_scalar, &group->order, ctx))
1428 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1431 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1434 num_bytes = BN_bn2bin(p_scalar, tmp);
1435 flip_endian(secrets[i], tmp, num_bytes);
1436 /* precompute multiples */
1437 if ((!BN_to_felem(x_out, &p->X)) ||
1438 (!BN_to_felem(y_out, &p->Y)) ||
1439 (!BN_to_felem(z_out, &p->Z))) goto err;
1440 felem_assign(pre_comp[i][1][0], x_out);
1441 felem_assign(pre_comp[i][1][1], y_out);
1442 felem_assign(pre_comp[i][1][2], z_out);
1443 for (j = 2; j <= 16; ++j)
1448 pre_comp[i][j][0], pre_comp[i][j][1], pre_comp[i][j][2],
1449 pre_comp[i][1][0], pre_comp[i][1][1], pre_comp[i][1][2],
1450 0, pre_comp[i][j-1][0], pre_comp[i][j-1][1], pre_comp[i][j-1][2]);
1455 pre_comp[i][j][0], pre_comp[i][j][1], pre_comp[i][j][2],
1456 pre_comp[i][j/2][0], pre_comp[i][j/2][1], pre_comp[i][j/2][2]);
1462 make_points_affine(num_points * 17, pre_comp[0], tmp_felems);
1465 /* the scalar for the generator */
1466 if ((scalar != NULL) && (have_pre_comp))
1468 memset(g_secret, 0, sizeof g_secret);
1469 /* reduce scalar to 0 <= scalar < 2^224 */
1470 if ((BN_num_bits(scalar) > 224) || (BN_is_negative(scalar)))
1472 /* this is an unusual input, and we don't guarantee
1473 * constant-timeness */
1474 if (!BN_nnmod(tmp_scalar, scalar, &group->order, ctx))
1476 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1479 num_bytes = BN_bn2bin(tmp_scalar, tmp);
1482 num_bytes = BN_bn2bin(scalar, tmp);
1483 flip_endian(g_secret, tmp, num_bytes);
1484 /* do the multiplication with generator precomputation*/
1485 batch_mul(x_out, y_out, z_out,
1486 (const felem_bytearray (*)) secrets, num_points,
1488 mixed, (const felem (*)[17][3]) pre_comp,
1492 /* do the multiplication without generator precomputation */
1493 batch_mul(x_out, y_out, z_out,
1494 (const felem_bytearray (*)) secrets, num_points,
1495 NULL, mixed, (const felem (*)[17][3]) pre_comp, NULL);
1496 /* reduce the output to its unique minimal representation */
1497 felem_contract(x_in, x_out);
1498 felem_contract(y_in, y_out);
1499 felem_contract(z_in, z_out);
1500 if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
1501 (!felem_to_BN(z, z_in)))
1503 ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
1506 ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
1510 if (generator != NULL)
1511 EC_POINT_free(generator);
1512 if (new_ctx != NULL)
1513 BN_CTX_free(new_ctx);
1514 if (secrets != NULL)
1515 OPENSSL_free(secrets);
1516 if (pre_comp != NULL)
1517 OPENSSL_free(pre_comp);
1518 if (tmp_felems != NULL)
1519 OPENSSL_free(tmp_felems);
1523 int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
1526 NISTP224_PRE_COMP *pre = NULL;
1528 BN_CTX *new_ctx = NULL;
1530 EC_POINT *generator = NULL;
1531 felem tmp_felems[32];
1533 /* throw away old precomputation */
1534 EC_EX_DATA_free_data(&group->extra_data, nistp224_pre_comp_dup,
1535 nistp224_pre_comp_free, nistp224_pre_comp_clear_free);
1537 if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
1539 if (((x = BN_CTX_get(ctx)) == NULL) ||
1540 ((y = BN_CTX_get(ctx)) == NULL))
1542 /* get the generator */
1543 if (group->generator == NULL) goto err;
1544 generator = EC_POINT_new(group);
1545 if (generator == NULL)
1547 BN_bin2bn(nistp224_curve_params[3], sizeof (felem_bytearray), x);
1548 BN_bin2bn(nistp224_curve_params[4], sizeof (felem_bytearray), y);
1549 if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
1551 if ((pre = nistp224_pre_comp_new()) == NULL)
1553 /* if the generator is the standard one, use built-in precomputation */
1554 if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
1556 memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
1560 if ((!BN_to_felem(pre->g_pre_comp[0][1][0], &group->generator->X)) ||
1561 (!BN_to_felem(pre->g_pre_comp[0][1][1], &group->generator->Y)) ||
1562 (!BN_to_felem(pre->g_pre_comp[0][1][2], &group->generator->Z)))
1564 /* compute 2^56*G, 2^112*G, 2^168*G for the first table,
1565 * 2^28*G, 2^84*G, 2^140*G, 2^196*G for the second one
1567 for (i = 1; i <= 8; i <<= 1)
1570 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2],
1571 pre->g_pre_comp[0][i][0], pre->g_pre_comp[0][i][1], pre->g_pre_comp[0][i][2]);
1572 for (j = 0; j < 27; ++j)
1575 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2],
1576 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1581 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2],
1582 pre->g_pre_comp[1][i][0], pre->g_pre_comp[1][i][1], pre->g_pre_comp[1][i][2]);
1583 for (j = 0; j < 27; ++j)
1586 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2],
1587 pre->g_pre_comp[0][2*i][0], pre->g_pre_comp[0][2*i][1], pre->g_pre_comp[0][2*i][2]);
1590 for (i = 0; i < 2; i++)
1592 /* g_pre_comp[i][0] is the point at infinity */
1593 memset(pre->g_pre_comp[i][0], 0, sizeof(pre->g_pre_comp[i][0]));
1594 /* the remaining multiples */
1595 /* 2^56*G + 2^112*G resp. 2^84*G + 2^140*G */
1597 pre->g_pre_comp[i][6][0], pre->g_pre_comp[i][6][1],
1598 pre->g_pre_comp[i][6][2], pre->g_pre_comp[i][4][0],
1599 pre->g_pre_comp[i][4][1], pre->g_pre_comp[i][4][2],
1600 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1601 pre->g_pre_comp[i][2][2]);
1602 /* 2^56*G + 2^168*G resp. 2^84*G + 2^196*G */
1604 pre->g_pre_comp[i][10][0], pre->g_pre_comp[i][10][1],
1605 pre->g_pre_comp[i][10][2], pre->g_pre_comp[i][8][0],
1606 pre->g_pre_comp[i][8][1], pre->g_pre_comp[i][8][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^112*G + 2^168*G resp. 2^140*G + 2^196*G */
1611 pre->g_pre_comp[i][12][0], pre->g_pre_comp[i][12][1],
1612 pre->g_pre_comp[i][12][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][4][0], pre->g_pre_comp[i][4][1],
1615 pre->g_pre_comp[i][4][2]);
1616 /* 2^56*G + 2^112*G + 2^168*G resp. 2^84*G + 2^140*G + 2^196*G */
1618 pre->g_pre_comp[i][14][0], pre->g_pre_comp[i][14][1],
1619 pre->g_pre_comp[i][14][2], pre->g_pre_comp[i][12][0],
1620 pre->g_pre_comp[i][12][1], pre->g_pre_comp[i][12][2],
1621 0, pre->g_pre_comp[i][2][0], pre->g_pre_comp[i][2][1],
1622 pre->g_pre_comp[i][2][2]);
1623 for (j = 1; j < 8; ++j)
1625 /* odd multiples: add G resp. 2^28*G */
1627 pre->g_pre_comp[i][2*j+1][0], pre->g_pre_comp[i][2*j+1][1],
1628 pre->g_pre_comp[i][2*j+1][2], pre->g_pre_comp[i][2*j][0],
1629 pre->g_pre_comp[i][2*j][1], pre->g_pre_comp[i][2*j][2],
1630 0, pre->g_pre_comp[i][1][0], pre->g_pre_comp[i][1][1],
1631 pre->g_pre_comp[i][1][2]);
1634 make_points_affine(31, &(pre->g_pre_comp[0][1]), tmp_felems);
1636 if (!EC_EX_DATA_set_data(&group->extra_data, pre, nistp224_pre_comp_dup,
1637 nistp224_pre_comp_free, nistp224_pre_comp_clear_free))
1643 if (generator != NULL)
1644 EC_POINT_free(generator);
1645 if (new_ctx != NULL)
1646 BN_CTX_free(new_ctx);
1648 nistp224_pre_comp_free(pre);
1652 int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
1654 if (EC_EX_DATA_get_data(group->extra_data, nistp224_pre_comp_dup,
1655 nistp224_pre_comp_free, nistp224_pre_comp_clear_free)
1663 static void *dummy=&dummy;