michael@0: /* This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  * License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
michael@0: 
michael@0: #include "ecp_fp.h"
michael@0: #include "ecl-priv.h"
michael@0: #include <stdlib.h>
michael@0: 
michael@0: /* Performs tidying on a short multi-precision floating point integer (the 
michael@0:  * lower group->numDoubles floats). */
michael@0: void
michael@0: ecfp_tidyShort(double *t, const EC_group_fp * group)
michael@0: {
michael@0: 	group->ecfp_tidy(t, group->alpha, group);
michael@0: }
michael@0: 
michael@0: /* Performs tidying on only the upper float digits of a multi-precision
michael@0:  * floating point integer, i.e. the digits beyond the regular length which 
michael@0:  * are removed in the reduction step. */
michael@0: void
michael@0: ecfp_tidyUpper(double *t, const EC_group_fp * group)
michael@0: {
michael@0: 	group->ecfp_tidy(t + group->numDoubles,
michael@0: 					 group->alpha + group->numDoubles, group);
michael@0: }
michael@0: 
michael@0: /* Performs a "tidy" operation, which performs carrying, moving excess
michael@0:  * bits from one double to the next double, so that the precision of the
michael@0:  * doubles is reduced to the regular precision group->doubleBitSize. This
michael@0:  * might result in some float digits being negative. Alternative C version 
michael@0:  * for portability. */
michael@0: void
michael@0: ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group)
michael@0: {
michael@0: 	double q;
michael@0: 	int i;
michael@0: 
michael@0: 	/* Do carrying */
michael@0: 	for (i = 0; i < group->numDoubles - 1; i++) {
michael@0: 		q = t[i] + alpha[i + 1];
michael@0: 		q -= alpha[i + 1];
michael@0: 		t[i] -= q;
michael@0: 		t[i + 1] += q;
michael@0: 
michael@0: 		/* If we don't assume that truncation rounding is used, then q
michael@0: 		 * might be 2^n bigger than expected (if it rounds up), then t[0]
michael@0: 		 * could be negative and t[1] 2^n larger than expected. */
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Performs a more mathematically precise "tidying" so that each term is
michael@0:  * positive.  This is slower than the regular tidying, and is used for
michael@0:  * conversion from floating point to integer. */
michael@0: void
michael@0: ecfp_positiveTidy(double *t, const EC_group_fp * group)
michael@0: {
michael@0: 	double q;
michael@0: 	int i;
michael@0: 
michael@0: 	/* Do carrying */
michael@0: 	for (i = 0; i < group->numDoubles - 1; i++) {
michael@0: 		/* Subtract beta to force rounding down */
michael@0: 		q = t[i] - ecfp_beta[i + 1];
michael@0: 		q += group->alpha[i + 1];
michael@0: 		q -= group->alpha[i + 1];
michael@0: 		t[i] -= q;
michael@0: 		t[i + 1] += q;
michael@0: 
michael@0: 		/* Due to subtracting ecfp_beta, we should have each term a
michael@0: 		 * non-negative int */
michael@0: 		ECFP_ASSERT(t[i] / ecfp_exp[i] ==
michael@0: 					(unsigned long long) (t[i] / ecfp_exp[i]));
michael@0: 		ECFP_ASSERT(t[i] >= 0);
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Converts from a floating point representation into an mp_int. Expects
michael@0:  * that d is already reduced. */
michael@0: void
michael@0: ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup)
michael@0: {
michael@0: 	EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
michael@0: 	unsigned short i16[(group->primeBitSize + 15) / 16];
michael@0: 	double q = 1;
michael@0: 
michael@0: #ifdef ECL_THIRTY_TWO_BIT
michael@0: 	/* TEST uint32_t z = 0; */
michael@0: 	unsigned int z = 0;
michael@0: #else
michael@0: 	uint64_t z = 0;
michael@0: #endif
michael@0: 	int zBits = 0;
michael@0: 	int copiedBits = 0;
michael@0: 	int i = 0;
michael@0: 	int j = 0;
michael@0: 
michael@0: 	mp_digit *out;
michael@0: 
michael@0: 	/* Result should always be >= 0, so set sign accordingly */
michael@0: 	MP_SIGN(mpout) = MP_ZPOS;
michael@0: 
michael@0: 	/* Tidy up so we're just dealing with positive numbers */
michael@0: 	ecfp_positiveTidy(d, group);
michael@0: 
michael@0: 	/* We might need to do this reduction step more than once if the
michael@0: 	 * reduction adds smaller terms which carry-over to cause another
michael@0: 	 * reduction. However, this should happen very rarely, if ever,
michael@0: 	 * depending on the elliptic curve. */
michael@0: 	do {
michael@0: 		/* Init loop data */
michael@0: 		z = 0;
michael@0: 		zBits = 0;
michael@0: 		q = 1;
michael@0: 		i = 0;
michael@0: 		j = 0;
michael@0: 		copiedBits = 0;
michael@0: 
michael@0: 		/* Might have to do a bit more reduction */
michael@0: 		group->ecfp_singleReduce(d, group);
michael@0: 
michael@0: 		/* Grow the size of the mpint if it's too small */
michael@0: 		s_mp_grow(mpout, group->numInts);
michael@0: 		MP_USED(mpout) = group->numInts;
michael@0: 		out = MP_DIGITS(mpout);
michael@0: 
michael@0: 		/* Convert double to 16 bit integers */
michael@0: 		while (copiedBits < group->primeBitSize) {
michael@0: 			if (zBits < 16) {
michael@0: 				z += d[i] * q;
michael@0: 				i++;
michael@0: 				ECFP_ASSERT(i < (group->primeBitSize + 15) / 16);
michael@0: 				zBits += group->doubleBitSize;
michael@0: 			}
michael@0: 			i16[j] = z;
michael@0: 			j++;
michael@0: 			z >>= 16;
michael@0: 			zBits -= 16;
michael@0: 			q *= ecfp_twom16;
michael@0: 			copiedBits += 16;
michael@0: 		}
michael@0: 	} while (z != 0);
michael@0: 
michael@0: 	/* Convert 16 bit integers to mp_digit */
michael@0: #ifdef ECL_THIRTY_TWO_BIT
michael@0: 	for (i = 0; i < (group->primeBitSize + 15) / 16; i += 2) {
michael@0: 		*out = 0;
michael@0: 		if (i + 1 < (group->primeBitSize + 15) / 16) {
michael@0: 			*out = i16[i + 1];
michael@0: 			*out <<= 16;
michael@0: 		}
michael@0: 		*out++ += i16[i];
michael@0: 	}
michael@0: #else							/* 64 bit */
michael@0: 	for (i = 0; i < (group->primeBitSize + 15) / 16; i += 4) {
michael@0: 		*out = 0;
michael@0: 		if (i + 3 < (group->primeBitSize + 15) / 16) {
michael@0: 			*out = i16[i + 3];
michael@0: 			*out <<= 16;
michael@0: 		}
michael@0: 		if (i + 2 < (group->primeBitSize + 15) / 16) {
michael@0: 			*out += i16[i + 2];
michael@0: 			*out <<= 16;
michael@0: 		}
michael@0: 		if (i + 1 < (group->primeBitSize + 15) / 16) {
michael@0: 			*out += i16[i + 1];
michael@0: 			*out <<= 16;
michael@0: 		}
michael@0: 		*out++ += i16[i];
michael@0: 	}
michael@0: #endif
michael@0: 
michael@0: 	/* Perform final reduction.  mpout should already be the same number
michael@0: 	 * of bits as p, but might not be less than p.  Make it so. Since
michael@0: 	 * mpout has the same number of bits as p, and 2p has a larger bit
michael@0: 	 * size, then mpout < 2p, so a single subtraction of p will suffice. */
michael@0: 	if (mp_cmp(mpout, &ecgroup->meth->irr) >= 0) {
michael@0: 		mp_sub(mpout, &ecgroup->meth->irr, mpout);
michael@0: 	}
michael@0: 
michael@0: 	/* Shrink the size of the mp_int to the actual used size (required for 
michael@0: 	 * mp_cmp_z == 0) */
michael@0: 	out = MP_DIGITS(mpout);
michael@0: 	for (i = group->numInts - 1; i > 0; i--) {
michael@0: 		if (out[i] != 0)
michael@0: 			break;
michael@0: 	}
michael@0: 	MP_USED(mpout) = i + 1;
michael@0: 
michael@0: 	/* Should be between 0 and p-1 */
michael@0: 	ECFP_ASSERT(mp_cmp(mpout, &ecgroup->meth->irr) < 0);
michael@0: 	ECFP_ASSERT(mp_cmp_z(mpout) >= 0);
michael@0: }
michael@0: 
michael@0: /* Converts from an mpint into a floating point representation. */
michael@0: void
michael@0: ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup)
michael@0: {
michael@0: 	int i;
michael@0: 	int j = 0;
michael@0: 	int size;
michael@0: 	double shift = 1;
michael@0: 	mp_digit *in;
michael@0: 	EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
michael@0: 
michael@0: #ifdef ECL_DEBUG
michael@0: 	/* if debug mode, convert result back using ecfp_fp2i into cmp, then
michael@0: 	 * compare to x. */
michael@0: 	mp_int cmp;
michael@0: 
michael@0: 	MP_DIGITS(&cmp) = NULL;
michael@0: 	mp_init(&cmp);
michael@0: #endif
michael@0: 
michael@0: 	ECFP_ASSERT(group != NULL);
michael@0: 
michael@0: 	/* init output to 0 (since we skip over some terms) */
michael@0: 	for (i = 0; i < group->numDoubles; i++)
michael@0: 		out[i] = 0;
michael@0: 	i = 0;
michael@0: 
michael@0: 	size = MP_USED(x);
michael@0: 	in = MP_DIGITS(x);
michael@0: 
michael@0: 	/* Copy from int into doubles */
michael@0: #ifdef ECL_THIRTY_TWO_BIT
michael@0: 	while (j < size) {
michael@0: 		while (group->doubleBitSize * (i + 1) <= 32 * j) {
michael@0: 			i++;
michael@0: 		}
michael@0: 		ECFP_ASSERT(group->doubleBitSize * i <= 32 * j);
michael@0: 		out[i] = in[j];
michael@0: 		out[i] *= shift;
michael@0: 		shift *= ecfp_two32;
michael@0: 		j++;
michael@0: 	}
michael@0: #else
michael@0: 	while (j < size) {
michael@0: 		while (group->doubleBitSize * (i + 1) <= 64 * j) {
michael@0: 			i++;
michael@0: 		}
michael@0: 		ECFP_ASSERT(group->doubleBitSize * i <= 64 * j);
michael@0: 		out[i] = (in[j] & 0x00000000FFFFFFFF) * shift;
michael@0: 
michael@0: 		while (group->doubleBitSize * (i + 1) <= 64 * j + 32) {
michael@0: 			i++;
michael@0: 		}
michael@0: 		ECFP_ASSERT(24 * i <= 64 * j + 32);
michael@0: 		out[i] = (in[j] & 0xFFFFFFFF00000000) * shift;
michael@0: 
michael@0: 		shift *= ecfp_two64;
michael@0: 		j++;
michael@0: 	}
michael@0: #endif
michael@0: 	/* Realign bits to match double boundaries */
michael@0: 	ecfp_tidyShort(out, group);
michael@0: 
michael@0: #ifdef ECL_DEBUG
michael@0: 	/* Convert result back to mp_int, compare to original */
michael@0: 	ecfp_fp2i(&cmp, out, ecgroup);
michael@0: 	ECFP_ASSERT(mp_cmp(&cmp, x) == 0);
michael@0: 	mp_clear(&cmp);
michael@0: #endif
michael@0: }
michael@0: 
michael@0: /* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters
michael@0:  * a, b and p are the elliptic curve coefficients and the prime that
michael@0:  * determines the field GFp.  Elliptic curve points P and R can be
michael@0:  * identical.  Uses Jacobian coordinates. Uses 4-bit window method. */
michael@0: mp_err
michael@0: ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px,
michael@0: 						   const mp_int *py, mp_int *rx, mp_int *ry,
michael@0: 						   const ECGroup *ecgroup)
michael@0: {
michael@0: 	mp_err res = MP_OKAY;
michael@0: 	ecfp_jac_pt precomp[16], r;
michael@0: 	ecfp_aff_pt p;
michael@0: 	EC_group_fp *group;
michael@0: 
michael@0: 	mp_int rz;
michael@0: 	int i, ni, d;
michael@0: 
michael@0: 	ARGCHK(ecgroup != NULL, MP_BADARG);
michael@0: 	ARGCHK((n != NULL) && (px != NULL) && (py != NULL), MP_BADARG);
michael@0: 
michael@0: 	group = (EC_group_fp *) ecgroup->extra1;
michael@0: 	MP_DIGITS(&rz) = 0;
michael@0: 	MP_CHECKOK(mp_init(&rz));
michael@0: 
michael@0: 	/* init p, da */
michael@0: 	ecfp_i2fp(p.x, px, ecgroup);
michael@0: 	ecfp_i2fp(p.y, py, ecgroup);
michael@0: 	ecfp_i2fp(group->curvea, &ecgroup->curvea, ecgroup);
michael@0: 
michael@0: 	/* Do precomputation */
michael@0: 	group->precompute_jac(precomp, &p, group);
michael@0: 
michael@0: 	/* Do main body of calculations */
michael@0: 	d = (mpl_significant_bits(n) + 3) / 4;
michael@0: 
michael@0: 	/* R = inf */
michael@0: 	for (i = 0; i < group->numDoubles; i++) {
michael@0: 		r.z[i] = 0;
michael@0: 	}
michael@0: 
michael@0: 	for (i = d - 1; i >= 0; i--) {
michael@0: 		/* compute window ni */
michael@0: 		ni = MP_GET_BIT(n, 4 * i + 3);
michael@0: 		ni <<= 1;
michael@0: 		ni |= MP_GET_BIT(n, 4 * i + 2);
michael@0: 		ni <<= 1;
michael@0: 		ni |= MP_GET_BIT(n, 4 * i + 1);
michael@0: 		ni <<= 1;
michael@0: 		ni |= MP_GET_BIT(n, 4 * i);
michael@0: 
michael@0: 		/* R = 2^4 * R */
michael@0: 		group->pt_dbl_jac(&r, &r, group);
michael@0: 		group->pt_dbl_jac(&r, &r, group);
michael@0: 		group->pt_dbl_jac(&r, &r, group);
michael@0: 		group->pt_dbl_jac(&r, &r, group);
michael@0: 
michael@0: 		/* R = R + (ni * P) */
michael@0: 		group->pt_add_jac(&r, &precomp[ni], &r, group);
michael@0: 	}
michael@0: 
michael@0: 	/* Convert back to integer */
michael@0: 	ecfp_fp2i(rx, r.x, ecgroup);
michael@0: 	ecfp_fp2i(ry, r.y, ecgroup);
michael@0: 	ecfp_fp2i(&rz, r.z, ecgroup);
michael@0: 
michael@0: 	/* convert result S to affine coordinates */
michael@0: 	MP_CHECKOK(ec_GFp_pt_jac2aff(rx, ry, &rz, rx, ry, ecgroup));
michael@0: 
michael@0:   CLEANUP:
michael@0: 	mp_clear(&rz);
michael@0: 	return res;
michael@0: }
michael@0: 
michael@0: /* Uses mixed Jacobian-affine coordinates to perform a point
michael@0:  * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine
michael@0:  * coordinates (Jacobian coordinates for doubles and affine coordinates
michael@0:  * for additions; based on recommendation from Brown et al.). Not very
michael@0:  * time efficient but quite space efficient, no precomputation needed.
michael@0:  * group contains the elliptic curve coefficients and the prime that
michael@0:  * determines the field GFp.  Elliptic curve points P and R can be
michael@0:  * identical. Performs calculations in floating point number format, since 
michael@0:  * this is faster than the integer operations on the ULTRASPARC III.
michael@0:  * Uses left-to-right binary method (double & add) (algorithm 9) for
michael@0:  * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes.
michael@0:  * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */
michael@0: mp_err
michael@0: ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py,
michael@0: 					 mp_int *rx, mp_int *ry, const ECGroup *ecgroup)
michael@0: {
michael@0: 	mp_err res;
michael@0: 	mp_int sx, sy, sz;
michael@0: 
michael@0: 	ecfp_aff_pt p;
michael@0: 	ecfp_jac_pt r;
michael@0: 	EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
michael@0: 
michael@0: 	int i, l;
michael@0: 
michael@0: 	MP_DIGITS(&sx) = 0;
michael@0: 	MP_DIGITS(&sy) = 0;
michael@0: 	MP_DIGITS(&sz) = 0;
michael@0: 	MP_CHECKOK(mp_init(&sx));
michael@0: 	MP_CHECKOK(mp_init(&sy));
michael@0: 	MP_CHECKOK(mp_init(&sz));
michael@0: 
michael@0: 	/* if n = 0 then r = inf */
michael@0: 	if (mp_cmp_z(n) == 0) {
michael@0: 		mp_zero(rx);
michael@0: 		mp_zero(ry);
michael@0: 		res = MP_OKAY;
michael@0: 		goto CLEANUP;
michael@0: 		/* if n < 0 then out of range error */
michael@0: 	} else if (mp_cmp_z(n) < 0) {
michael@0: 		res = MP_RANGE;
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* Convert from integer to floating point */
michael@0: 	ecfp_i2fp(p.x, px, ecgroup);
michael@0: 	ecfp_i2fp(p.y, py, ecgroup);
michael@0: 	ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
michael@0: 
michael@0: 	/* Init r to point at infinity */
michael@0: 	for (i = 0; i < group->numDoubles; i++) {
michael@0: 		r.z[i] = 0;
michael@0: 	}
michael@0: 
michael@0: 	/* double and add method */
michael@0: 	l = mpl_significant_bits(n) - 1;
michael@0: 
michael@0: 	for (i = l; i >= 0; i--) {
michael@0: 		/* R = 2R */
michael@0: 		group->pt_dbl_jac(&r, &r, group);
michael@0: 
michael@0: 		/* if n_i = 1, then R = R + Q */
michael@0: 		if (MP_GET_BIT(n, i) != 0) {
michael@0: 			group->pt_add_jac_aff(&r, &p, &r, group);
michael@0: 		}
michael@0: 	}
michael@0: 
michael@0: 	/* Convert from floating point to integer */
michael@0: 	ecfp_fp2i(&sx, r.x, ecgroup);
michael@0: 	ecfp_fp2i(&sy, r.y, ecgroup);
michael@0: 	ecfp_fp2i(&sz, r.z, ecgroup);
michael@0: 
michael@0: 	/* convert result R to affine coordinates */
michael@0: 	MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
michael@0: 
michael@0:   CLEANUP:
michael@0: 	mp_clear(&sx);
michael@0: 	mp_clear(&sy);
michael@0: 	mp_clear(&sz);
michael@0: 	return res;
michael@0: }
michael@0: 
michael@0: /* Computes R = nP where R is (rx, ry) and P is the base point. Elliptic
michael@0:  * curve points P and R can be identical. Uses mixed Modified-Jacobian
michael@0:  * co-ordinates for doubling and Chudnovsky Jacobian coordinates for
michael@0:  * additions. Uses 5-bit window NAF method (algorithm 11) for scalar-point 
michael@0:  * multiplication from Brown, Hankerson, Lopez, Menezes. Software
michael@0:  * Implementation of the NIST Elliptic Curves Over Prime Fields. */
michael@0: mp_err
michael@0: ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px,
michael@0: 						 const mp_int *py, mp_int *rx, mp_int *ry,
michael@0: 						 const ECGroup *ecgroup)
michael@0: {
michael@0: 	mp_err res = MP_OKAY;
michael@0: 	mp_int sx, sy, sz;
michael@0: 	EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
michael@0: 	ecfp_chud_pt precomp[16];
michael@0: 
michael@0: 	ecfp_aff_pt p;
michael@0: 	ecfp_jm_pt r;
michael@0: 
michael@0: 	signed char naf[group->orderBitSize + 1];
michael@0: 	int i;
michael@0: 
michael@0: 	MP_DIGITS(&sx) = 0;
michael@0: 	MP_DIGITS(&sy) = 0;
michael@0: 	MP_DIGITS(&sz) = 0;
michael@0: 	MP_CHECKOK(mp_init(&sx));
michael@0: 	MP_CHECKOK(mp_init(&sy));
michael@0: 	MP_CHECKOK(mp_init(&sz));
michael@0: 
michael@0: 	/* if n = 0 then r = inf */
michael@0: 	if (mp_cmp_z(n) == 0) {
michael@0: 		mp_zero(rx);
michael@0: 		mp_zero(ry);
michael@0: 		res = MP_OKAY;
michael@0: 		goto CLEANUP;
michael@0: 		/* if n < 0 then out of range error */
michael@0: 	} else if (mp_cmp_z(n) < 0) {
michael@0: 		res = MP_RANGE;
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* Convert from integer to floating point */
michael@0: 	ecfp_i2fp(p.x, px, ecgroup);
michael@0: 	ecfp_i2fp(p.y, py, ecgroup);
michael@0: 	ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
michael@0: 
michael@0: 	/* Perform precomputation */
michael@0: 	group->precompute_chud(precomp, &p, group);
michael@0: 
michael@0: 	/* Compute 5NAF */
michael@0: 	ec_compute_wNAF(naf, group->orderBitSize, n, 5);
michael@0: 
michael@0: 	/* Init R = pt at infinity */
michael@0: 	for (i = 0; i < group->numDoubles; i++) {
michael@0: 		r.z[i] = 0;
michael@0: 	}
michael@0: 
michael@0: 	/* wNAF method */
michael@0: 	for (i = group->orderBitSize; i >= 0; i--) {
michael@0: 		/* R = 2R */
michael@0: 		group->pt_dbl_jm(&r, &r, group);
michael@0: 
michael@0: 		if (naf[i] != 0) {
michael@0: 			group->pt_add_jm_chud(&r, &precomp[(naf[i] + 15) / 2], &r,
michael@0: 								  group);
michael@0: 		}
michael@0: 	}
michael@0: 
michael@0: 	/* Convert from floating point to integer */
michael@0: 	ecfp_fp2i(&sx, r.x, ecgroup);
michael@0: 	ecfp_fp2i(&sy, r.y, ecgroup);
michael@0: 	ecfp_fp2i(&sz, r.z, ecgroup);
michael@0: 
michael@0: 	/* convert result R to affine coordinates */
michael@0: 	MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
michael@0: 
michael@0:   CLEANUP:
michael@0: 	mp_clear(&sx);
michael@0: 	mp_clear(&sy);
michael@0: 	mp_clear(&sz);
michael@0: 	return res;
michael@0: }
michael@0: 
michael@0: /* Cleans up extra memory allocated in ECGroup for this implementation. */
michael@0: void
michael@0: ec_GFp_extra_free_fp(ECGroup *group)
michael@0: {
michael@0: 	if (group->extra1 != NULL) {
michael@0: 		free(group->extra1);
michael@0: 		group->extra1 = NULL;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Tests what precision floating point arithmetic is set to. This should
michael@0:  * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa
michael@0:  * (extended precision on x86) and sets it into the EC_group_fp. Returns
michael@0:  * either 53 or 64 accordingly. */
michael@0: int
michael@0: ec_set_fp_precision(EC_group_fp * group)
michael@0: {
michael@0: 	double a = 9007199254740992.0;	/* 2^53 */
michael@0: 	double b = a + 1;
michael@0: 
michael@0: 	if (a == b) {
michael@0: 		group->fpPrecision = 53;
michael@0: 		group->alpha = ecfp_alpha_53;
michael@0: 		return 53;
michael@0: 	}
michael@0: 	group->fpPrecision = 64;
michael@0: 	group->alpha = ecfp_alpha_64;
michael@0: 	return 64;
michael@0: }