The Tor Browser / file comparison

comparison: security/nss/lib/freebl/ecl/ecp_fp.c

security/nss/lib/freebl/ecl/ecp_fp.c

branch
TOR_BUG_9701
changeset 15
b8a032363ba2
equal deleted inserted replaced
-1:000000000000 0:e23a0103e589
1 /* This Source Code Form is subject to the terms of the Mozilla Public
2 * License, v. 2.0. If a copy of the MPL was not distributed with this
3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
4
5 #include "ecp_fp.h"
6 #include "ecl-priv.h"
7 #include <stdlib.h>
8
9 /* Performs tidying on a short multi-precision floating point integer (the
10 * lower group->numDoubles floats). */
11 void
12 ecfp_tidyShort(double *t, const EC_group_fp * group)
13 {
14 group->ecfp_tidy(t, group->alpha, group);
15 }
16
17 /* Performs tidying on only the upper float digits of a multi-precision
18 * floating point integer, i.e. the digits beyond the regular length which
19 * are removed in the reduction step. */
20 void
21 ecfp_tidyUpper(double *t, const EC_group_fp * group)
22 {
23 group->ecfp_tidy(t + group->numDoubles,
24 group->alpha + group->numDoubles, group);
25 }
26
27 /* Performs a "tidy" operation, which performs carrying, moving excess
28 * bits from one double to the next double, so that the precision of the
29 * doubles is reduced to the regular precision group->doubleBitSize. This
30 * might result in some float digits being negative. Alternative C version
31 * for portability. */
32 void
33 ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group)
34 {
35 double q;
36 int i;
37
38 /* Do carrying */
39 for (i = 0; i < group->numDoubles - 1; i++) {
40 q = t[i] + alpha[i + 1];
41 q -= alpha[i + 1];
42 t[i] -= q;
43 t[i + 1] += q;
44
45 /* If we don't assume that truncation rounding is used, then q
46 * might be 2^n bigger than expected (if it rounds up), then t[0]
47 * could be negative and t[1] 2^n larger than expected. */
48 }
49 }
50
51 /* Performs a more mathematically precise "tidying" so that each term is
52 * positive. This is slower than the regular tidying, and is used for
53 * conversion from floating point to integer. */
54 void
55 ecfp_positiveTidy(double *t, const EC_group_fp * group)
56 {
57 double q;
58 int i;
59
60 /* Do carrying */
61 for (i = 0; i < group->numDoubles - 1; i++) {
62 /* Subtract beta to force rounding down */
63 q = t[i] - ecfp_beta[i + 1];
64 q += group->alpha[i + 1];
65 q -= group->alpha[i + 1];
66 t[i] -= q;
67 t[i + 1] += q;
68
69 /* Due to subtracting ecfp_beta, we should have each term a
70 * non-negative int */
71 ECFP_ASSERT(t[i] / ecfp_exp[i] ==
72 (unsigned long long) (t[i] / ecfp_exp[i]));
73 ECFP_ASSERT(t[i] >= 0);
74 }
75 }
76
77 /* Converts from a floating point representation into an mp_int. Expects
78 * that d is already reduced. */
79 void
80 ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup)
81 {
82 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
83 unsigned short i16[(group->primeBitSize + 15) / 16];
84 double q = 1;
85
86 #ifdef ECL_THIRTY_TWO_BIT
87 /* TEST uint32_t z = 0; */
88 unsigned int z = 0;
89 #else
90 uint64_t z = 0;
91 #endif
92 int zBits = 0;
93 int copiedBits = 0;
94 int i = 0;
95 int j = 0;
96
97 mp_digit *out;
98
99 /* Result should always be >= 0, so set sign accordingly */
100 MP_SIGN(mpout) = MP_ZPOS;
101
102 /* Tidy up so we're just dealing with positive numbers */
103 ecfp_positiveTidy(d, group);
104
105 /* We might need to do this reduction step more than once if the
106 * reduction adds smaller terms which carry-over to cause another
107 * reduction. However, this should happen very rarely, if ever,
108 * depending on the elliptic curve. */
109 do {
110 /* Init loop data */
111 z = 0;
112 zBits = 0;
113 q = 1;
114 i = 0;
115 j = 0;
116 copiedBits = 0;
117
118 /* Might have to do a bit more reduction */
119 group->ecfp_singleReduce(d, group);
120
121 /* Grow the size of the mpint if it's too small */
122 s_mp_grow(mpout, group->numInts);
123 MP_USED(mpout) = group->numInts;
124 out = MP_DIGITS(mpout);
125
126 /* Convert double to 16 bit integers */
127 while (copiedBits < group->primeBitSize) {
128 if (zBits < 16) {
129 z += d[i] * q;
130 i++;
131 ECFP_ASSERT(i < (group->primeBitSize + 15) / 16);
132 zBits += group->doubleBitSize;
133 }
134 i16[j] = z;
135 j++;
136 z >>= 16;
137 zBits -= 16;
138 q *= ecfp_twom16;
139 copiedBits += 16;
140 }
141 } while (z != 0);
142
143 /* Convert 16 bit integers to mp_digit */
144 #ifdef ECL_THIRTY_TWO_BIT
145 for (i = 0; i < (group->primeBitSize + 15) / 16; i += 2) {
146 *out = 0;
147 if (i + 1 < (group->primeBitSize + 15) / 16) {
148 *out = i16[i + 1];
149 *out <<= 16;
150 }
151 *out++ += i16[i];
152 }
153 #else /* 64 bit */
154 for (i = 0; i < (group->primeBitSize + 15) / 16; i += 4) {
155 *out = 0;
156 if (i + 3 < (group->primeBitSize + 15) / 16) {
157 *out = i16[i + 3];
158 *out <<= 16;
159 }
160 if (i + 2 < (group->primeBitSize + 15) / 16) {
161 *out += i16[i + 2];
162 *out <<= 16;
163 }
164 if (i + 1 < (group->primeBitSize + 15) / 16) {
165 *out += i16[i + 1];
166 *out <<= 16;
167 }
168 *out++ += i16[i];
169 }
170 #endif
171
172 /* Perform final reduction. mpout should already be the same number
173 * of bits as p, but might not be less than p. Make it so. Since
174 * mpout has the same number of bits as p, and 2p has a larger bit
175 * size, then mpout < 2p, so a single subtraction of p will suffice. */
176 if (mp_cmp(mpout, &ecgroup->meth->irr) >= 0) {
177 mp_sub(mpout, &ecgroup->meth->irr, mpout);
178 }
179
180 /* Shrink the size of the mp_int to the actual used size (required for
181 * mp_cmp_z == 0) */
182 out = MP_DIGITS(mpout);
183 for (i = group->numInts - 1; i > 0; i--) {
184 if (out[i] != 0)
185 break;
186 }
187 MP_USED(mpout) = i + 1;
188
189 /* Should be between 0 and p-1 */
190 ECFP_ASSERT(mp_cmp(mpout, &ecgroup->meth->irr) < 0);
191 ECFP_ASSERT(mp_cmp_z(mpout) >= 0);
192 }
193
194 /* Converts from an mpint into a floating point representation. */
195 void
196 ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup)
197 {
198 int i;
199 int j = 0;
200 int size;
201 double shift = 1;
202 mp_digit *in;
203 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
204
205 #ifdef ECL_DEBUG
206 /* if debug mode, convert result back using ecfp_fp2i into cmp, then
207 * compare to x. */
208 mp_int cmp;
209
210 MP_DIGITS(&cmp) = NULL;
211 mp_init(&cmp);
212 #endif
213
214 ECFP_ASSERT(group != NULL);
215
216 /* init output to 0 (since we skip over some terms) */
217 for (i = 0; i < group->numDoubles; i++)
218 out[i] = 0;
219 i = 0;
220
221 size = MP_USED(x);
222 in = MP_DIGITS(x);
223
224 /* Copy from int into doubles */
225 #ifdef ECL_THIRTY_TWO_BIT
226 while (j < size) {
227 while (group->doubleBitSize * (i + 1) <= 32 * j) {
228 i++;
229 }
230 ECFP_ASSERT(group->doubleBitSize * i <= 32 * j);
231 out[i] = in[j];
232 out[i] *= shift;
233 shift *= ecfp_two32;
234 j++;
235 }
236 #else
237 while (j < size) {
238 while (group->doubleBitSize * (i + 1) <= 64 * j) {
239 i++;
240 }
241 ECFP_ASSERT(group->doubleBitSize * i <= 64 * j);
242 out[i] = (in[j] & 0x00000000FFFFFFFF) * shift;
243
244 while (group->doubleBitSize * (i + 1) <= 64 * j + 32) {
245 i++;
246 }
247 ECFP_ASSERT(24 * i <= 64 * j + 32);
248 out[i] = (in[j] & 0xFFFFFFFF00000000) * shift;
249
250 shift *= ecfp_two64;
251 j++;
252 }
253 #endif
254 /* Realign bits to match double boundaries */
255 ecfp_tidyShort(out, group);
256
257 #ifdef ECL_DEBUG
258 /* Convert result back to mp_int, compare to original */
259 ecfp_fp2i(&cmp, out, ecgroup);
260 ECFP_ASSERT(mp_cmp(&cmp, x) == 0);
261 mp_clear(&cmp);
262 #endif
263 }
264
265 /* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters
266 * a, b and p are the elliptic curve coefficients and the prime that
267 * determines the field GFp. Elliptic curve points P and R can be
268 * identical. Uses Jacobian coordinates. Uses 4-bit window method. */
269 mp_err
270 ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px,
271 const mp_int *py, mp_int *rx, mp_int *ry,
272 const ECGroup *ecgroup)
273 {
274 mp_err res = MP_OKAY;
275 ecfp_jac_pt precomp[16], r;
276 ecfp_aff_pt p;
277 EC_group_fp *group;
278
279 mp_int rz;
280 int i, ni, d;
281
282 ARGCHK(ecgroup != NULL, MP_BADARG);
283 ARGCHK((n != NULL) && (px != NULL) && (py != NULL), MP_BADARG);
284
285 group = (EC_group_fp *) ecgroup->extra1;
286 MP_DIGITS(&rz) = 0;
287 MP_CHECKOK(mp_init(&rz));
288
289 /* init p, da */
290 ecfp_i2fp(p.x, px, ecgroup);
291 ecfp_i2fp(p.y, py, ecgroup);
292 ecfp_i2fp(group->curvea, &ecgroup->curvea, ecgroup);
293
294 /* Do precomputation */
295 group->precompute_jac(precomp, &p, group);
296
297 /* Do main body of calculations */
298 d = (mpl_significant_bits(n) + 3) / 4;
299
300 /* R = inf */
301 for (i = 0; i < group->numDoubles; i++) {
302 r.z[i] = 0;
303 }
304
305 for (i = d - 1; i >= 0; i--) {
306 /* compute window ni */
307 ni = MP_GET_BIT(n, 4 * i + 3);
308 ni <<= 1;
309 ni |= MP_GET_BIT(n, 4 * i + 2);
310 ni <<= 1;
311 ni |= MP_GET_BIT(n, 4 * i + 1);
312 ni <<= 1;
313 ni |= MP_GET_BIT(n, 4 * i);
314
315 /* R = 2^4 * R */
316 group->pt_dbl_jac(&r, &r, group);
317 group->pt_dbl_jac(&r, &r, group);
318 group->pt_dbl_jac(&r, &r, group);
319 group->pt_dbl_jac(&r, &r, group);
320
321 /* R = R + (ni * P) */
322 group->pt_add_jac(&r, &precomp[ni], &r, group);
323 }
324
325 /* Convert back to integer */
326 ecfp_fp2i(rx, r.x, ecgroup);
327 ecfp_fp2i(ry, r.y, ecgroup);
328 ecfp_fp2i(&rz, r.z, ecgroup);
329
330 /* convert result S to affine coordinates */
331 MP_CHECKOK(ec_GFp_pt_jac2aff(rx, ry, &rz, rx, ry, ecgroup));
332
333 CLEANUP:
334 mp_clear(&rz);
335 return res;
336 }
337
338 /* Uses mixed Jacobian-affine coordinates to perform a point
339 * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine
340 * coordinates (Jacobian coordinates for doubles and affine coordinates
341 * for additions; based on recommendation from Brown et al.). Not very
342 * time efficient but quite space efficient, no precomputation needed.
343 * group contains the elliptic curve coefficients and the prime that
344 * determines the field GFp. Elliptic curve points P and R can be
345 * identical. Performs calculations in floating point number format, since
346 * this is faster than the integer operations on the ULTRASPARC III.
347 * Uses left-to-right binary method (double & add) (algorithm 9) for
348 * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes.
349 * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */
350 mp_err
351 ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py,
352 mp_int *rx, mp_int *ry, const ECGroup *ecgroup)
353 {
354 mp_err res;
355 mp_int sx, sy, sz;
356
357 ecfp_aff_pt p;
358 ecfp_jac_pt r;
359 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
360
361 int i, l;
362
363 MP_DIGITS(&sx) = 0;
364 MP_DIGITS(&sy) = 0;
365 MP_DIGITS(&sz) = 0;
366 MP_CHECKOK(mp_init(&sx));
367 MP_CHECKOK(mp_init(&sy));
368 MP_CHECKOK(mp_init(&sz));
369
370 /* if n = 0 then r = inf */
371 if (mp_cmp_z(n) == 0) {
372 mp_zero(rx);
373 mp_zero(ry);
374 res = MP_OKAY;
375 goto CLEANUP;
376 /* if n < 0 then out of range error */
377 } else if (mp_cmp_z(n) < 0) {
378 res = MP_RANGE;
379 goto CLEANUP;
380 }
381
382 /* Convert from integer to floating point */
383 ecfp_i2fp(p.x, px, ecgroup);
384 ecfp_i2fp(p.y, py, ecgroup);
385 ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
386
387 /* Init r to point at infinity */
388 for (i = 0; i < group->numDoubles; i++) {
389 r.z[i] = 0;
390 }
391
392 /* double and add method */
393 l = mpl_significant_bits(n) - 1;
394
395 for (i = l; i >= 0; i--) {
396 /* R = 2R */
397 group->pt_dbl_jac(&r, &r, group);
398
399 /* if n_i = 1, then R = R + Q */
400 if (MP_GET_BIT(n, i) != 0) {
401 group->pt_add_jac_aff(&r, &p, &r, group);
402 }
403 }
404
405 /* Convert from floating point to integer */
406 ecfp_fp2i(&sx, r.x, ecgroup);
407 ecfp_fp2i(&sy, r.y, ecgroup);
408 ecfp_fp2i(&sz, r.z, ecgroup);
409
410 /* convert result R to affine coordinates */
411 MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
412
413 CLEANUP:
414 mp_clear(&sx);
415 mp_clear(&sy);
416 mp_clear(&sz);
417 return res;
418 }
419
420 /* Computes R = nP where R is (rx, ry) and P is the base point. Elliptic
421 * curve points P and R can be identical. Uses mixed Modified-Jacobian
422 * co-ordinates for doubling and Chudnovsky Jacobian coordinates for
423 * additions. Uses 5-bit window NAF method (algorithm 11) for scalar-point
424 * multiplication from Brown, Hankerson, Lopez, Menezes. Software
425 * Implementation of the NIST Elliptic Curves Over Prime Fields. */
426 mp_err
427 ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px,
428 const mp_int *py, mp_int *rx, mp_int *ry,
429 const ECGroup *ecgroup)
430 {
431 mp_err res = MP_OKAY;
432 mp_int sx, sy, sz;
433 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
434 ecfp_chud_pt precomp[16];
435
436 ecfp_aff_pt p;
437 ecfp_jm_pt r;
438
439 signed char naf[group->orderBitSize + 1];
440 int i;
441
442 MP_DIGITS(&sx) = 0;
443 MP_DIGITS(&sy) = 0;
444 MP_DIGITS(&sz) = 0;
445 MP_CHECKOK(mp_init(&sx));
446 MP_CHECKOK(mp_init(&sy));
447 MP_CHECKOK(mp_init(&sz));
448
449 /* if n = 0 then r = inf */
450 if (mp_cmp_z(n) == 0) {
451 mp_zero(rx);
452 mp_zero(ry);
453 res = MP_OKAY;
454 goto CLEANUP;
455 /* if n < 0 then out of range error */
456 } else if (mp_cmp_z(n) < 0) {
457 res = MP_RANGE;
458 goto CLEANUP;
459 }
460
461 /* Convert from integer to floating point */
462 ecfp_i2fp(p.x, px, ecgroup);
463 ecfp_i2fp(p.y, py, ecgroup);
464 ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
465
466 /* Perform precomputation */
467 group->precompute_chud(precomp, &p, group);
468
469 /* Compute 5NAF */
470 ec_compute_wNAF(naf, group->orderBitSize, n, 5);
471
472 /* Init R = pt at infinity */
473 for (i = 0; i < group->numDoubles; i++) {
474 r.z[i] = 0;
475 }
476
477 /* wNAF method */
478 for (i = group->orderBitSize; i >= 0; i--) {
479 /* R = 2R */
480 group->pt_dbl_jm(&r, &r, group);
481
482 if (naf[i] != 0) {
483 group->pt_add_jm_chud(&r, &precomp[(naf[i] + 15) / 2], &r,
484 group);
485 }
486 }
487
488 /* Convert from floating point to integer */
489 ecfp_fp2i(&sx, r.x, ecgroup);
490 ecfp_fp2i(&sy, r.y, ecgroup);
491 ecfp_fp2i(&sz, r.z, ecgroup);
492
493 /* convert result R to affine coordinates */
494 MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
495
496 CLEANUP:
497 mp_clear(&sx);
498 mp_clear(&sy);
499 mp_clear(&sz);
500 return res;
501 }
502
503 /* Cleans up extra memory allocated in ECGroup for this implementation. */
504 void
505 ec_GFp_extra_free_fp(ECGroup *group)
506 {
507 if (group->extra1 != NULL) {
508 free(group->extra1);
509 group->extra1 = NULL;
510 }
511 }
512
513 /* Tests what precision floating point arithmetic is set to. This should
514 * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa
515 * (extended precision on x86) and sets it into the EC_group_fp. Returns
516 * either 53 or 64 accordingly. */
517 int
518 ec_set_fp_precision(EC_group_fp * group)
519 {
520 double a = 9007199254740992.0; /* 2^53 */
521 double b = a + 1;
522
523 if (a == b) {
524 group->fpPrecision = 53;
525 group->alpha = ecfp_alpha_53;
526 return 53;
527 }
528 group->fpPrecision = 64;
529 group->alpha = ecfp_alpha_64;
530 return 64;
531 }

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial