PageRenderTime 48ms CodeModel.GetById 13ms app.highlight 29ms RepoModel.GetById 2ms app.codeStats 0ms

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

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