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: /* This source file is meant to be included by other source files
michael@0:  * (ecp_fp###.c, where ### is one of 160, 192, 224) and should not
michael@0:  * constitute an independent compilation unit. It requires the following
michael@0:  * preprocessor definitions be made: ECFP_BSIZE - the number of bits in
michael@0:  * the field's prime 
michael@0:  * ECFP_NUMDOUBLES - the number of doubles to store one
michael@0:  * multi-precision integer in floating point 
michael@0: 
michael@0: /* Adds a prefix to a given token to give a unique token name. Prefixes
michael@0:  * with "ecfp" + ECFP_BSIZE + "_". e.g. if ECFP_BSIZE = 160, then
michael@0:  * PREFIX(hello) = ecfp160_hello This optimization allows static function
michael@0:  * linking and compiler loop unrolling without code duplication. */
michael@0: #ifndef PREFIX
michael@0: #define PREFIX(b) PREFIX1(ECFP_BSIZE, b)
michael@0: #define PREFIX1(bsize, b) PREFIX2(bsize, b)
michael@0: #define PREFIX2(bsize, b) ecfp ## bsize ## _ ## b
michael@0: #endif
michael@0: 
michael@0: /* Returns true iff every double in d is 0. (If d == 0 and it is tidied,
michael@0:  * this will be true.) */
michael@0: mp_err PREFIX(isZero) (const double *d) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		if (d[i] != 0)
michael@0: 			return MP_NO;
michael@0: 	}
michael@0: 	return MP_YES;
michael@0: }
michael@0: 
michael@0: /* Sets the multi-precision floating point number at t = 0 */
michael@0: void PREFIX(zero) (double *t) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		t[i] = 0;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Sets the multi-precision floating point number at t = 1 */
michael@0: void PREFIX(one) (double *t) {
michael@0: 	int i;
michael@0: 
michael@0: 	t[0] = 1;
michael@0: 	for (i = 1; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		t[i] = 0;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Checks if point P(x, y, z) is at infinity. Uses Jacobian coordinates. */
michael@0: mp_err PREFIX(pt_is_inf_jac) (const ecfp_jac_pt * p) {
michael@0: 	return PREFIX(isZero) (p->z);
michael@0: }
michael@0: 
michael@0: /* Sets the Jacobian point P to be at infinity. */
michael@0: void PREFIX(set_pt_inf_jac) (ecfp_jac_pt * p) {
michael@0: 	PREFIX(zero) (p->z);
michael@0: }
michael@0: 
michael@0: /* Checks if point P(x, y) is at infinity. Uses Affine coordinates. */
michael@0: mp_err PREFIX(pt_is_inf_aff) (const ecfp_aff_pt * p) {
michael@0: 	if (PREFIX(isZero) (p->x) == MP_YES && PREFIX(isZero) (p->y) == MP_YES)
michael@0: 		return MP_YES;
michael@0: 	return MP_NO;
michael@0: }
michael@0: 
michael@0: /* Sets the affine point P to be at infinity. */
michael@0: void PREFIX(set_pt_inf_aff) (ecfp_aff_pt * p) {
michael@0: 	PREFIX(zero) (p->x);
michael@0: 	PREFIX(zero) (p->y);
michael@0: }
michael@0: 
michael@0: /* Checks if point P(x, y, z, a*z^4) is at infinity. Uses Modified
michael@0:  * Jacobian coordinates. */
michael@0: mp_err PREFIX(pt_is_inf_jm) (const ecfp_jm_pt * p) {
michael@0: 	return PREFIX(isZero) (p->z);
michael@0: }
michael@0: 
michael@0: /* Sets the Modified Jacobian point P to be at infinity. */
michael@0: void PREFIX(set_pt_inf_jm) (ecfp_jm_pt * p) {
michael@0: 	PREFIX(zero) (p->z);
michael@0: }
michael@0: 
michael@0: /* Checks if point P(x, y, z, z^2, z^3) is at infinity. Uses Chudnovsky
michael@0:  * Jacobian coordinates */
michael@0: mp_err PREFIX(pt_is_inf_chud) (const ecfp_chud_pt * p) {
michael@0: 	return PREFIX(isZero) (p->z);
michael@0: }
michael@0: 
michael@0: /* Sets the Chudnovsky Jacobian point P to be at infinity. */
michael@0: void PREFIX(set_pt_inf_chud) (ecfp_chud_pt * p) {
michael@0: 	PREFIX(zero) (p->z);
michael@0: }
michael@0: 
michael@0: /* Copies a multi-precision floating point number, Setting dest = src */
michael@0: void PREFIX(copy) (double *dest, const double *src) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		dest[i] = src[i];
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Sets dest = -src */
michael@0: void PREFIX(negLong) (double *dest, const double *src) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
michael@0: 		dest[i] = -src[i];
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Sets r = -p p = (x, y, z, z2, z3) r = (x, -y, z, z2, z3) Uses
michael@0:  * Chudnovsky Jacobian coordinates. */
michael@0: /* TODO reverse order */
michael@0: void PREFIX(pt_neg_chud) (const ecfp_chud_pt * p, ecfp_chud_pt * r) {
michael@0: 	int i;
michael@0: 
michael@0: 	PREFIX(copy) (r->x, p->x);
michael@0: 	PREFIX(copy) (r->z, p->z);
michael@0: 	PREFIX(copy) (r->z2, p->z2);
michael@0: 	PREFIX(copy) (r->z3, p->z3);
michael@0: 	for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		r->y[i] = -p->y[i];
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Computes r = x + y. Does not tidy or reduce. Any combinations of r, x,
michael@0:  * y can point to the same data. Componentwise adds first ECFP_NUMDOUBLES
michael@0:  * doubles of x and y and stores the result in r. */
michael@0: void PREFIX(addShort) (double *r, const double *x, const double *y) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		*r++ = *x++ + *y++;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Computes r = x + y. Does not tidy or reduce. Any combinations of r, x,
michael@0:  * y can point to the same data. Componentwise adds first
michael@0:  * 2*ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
michael@0: void PREFIX(addLong) (double *r, const double *x, const double *y) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
michael@0: 		*r++ = *x++ + *y++;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Computes r = x - y. Does not tidy or reduce. Any combinations of r, x,
michael@0:  * y can point to the same data. Componentwise subtracts first
michael@0:  * ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
michael@0: void PREFIX(subtractShort) (double *r, const double *x, const double *y) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0: 		*r++ = *x++ - *y++;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Computes r = x - y. Does not tidy or reduce. Any combinations of r, x,
michael@0:  * y can point to the same data. Componentwise subtracts first
michael@0:  * 2*ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
michael@0: void PREFIX(subtractLong) (double *r, const double *x, const double *y) {
michael@0: 	int i;
michael@0: 
michael@0: 	for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
michael@0: 		*r++ = *x++ - *y++;
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Computes r = x*y.  Both x and y should be tidied and reduced,
michael@0:  * r must be different (point to different memory) than x and y.
michael@0:  * Does not tidy or reduce. */
michael@0: void PREFIX(multiply)(double *r, const double *x, const double *y) {
michael@0: 	int i, j;
michael@0: 
michael@0: 	for(j=0;j<ECFP_NUMDOUBLES-1;j++) {
michael@0: 		r[j] = x[0] * y[j];
michael@0: 		r[j+(ECFP_NUMDOUBLES-1)] = x[ECFP_NUMDOUBLES-1] * y[j];
michael@0: 	}
michael@0: 	r[ECFP_NUMDOUBLES-1] = x[0] * y[ECFP_NUMDOUBLES-1];
michael@0: 	r[ECFP_NUMDOUBLES-1] += x[ECFP_NUMDOUBLES-1] * y[0];
michael@0: 	r[2*ECFP_NUMDOUBLES-2] = x[ECFP_NUMDOUBLES-1] * y[ECFP_NUMDOUBLES-1];
michael@0: 	r[2*ECFP_NUMDOUBLES-1] = 0;
michael@0: 	
michael@0: 	for(i=1;i<ECFP_NUMDOUBLES-1;i++) {
michael@0: 		for(j=0;j<ECFP_NUMDOUBLES;j++) {
michael@0: 			r[i+j] += (x[i] * y[j]);
michael@0: 		}
michael@0: 	}
michael@0: }
michael@0: 
michael@0: /* Computes the square of x and stores the result in r.  x should be
michael@0:  * tidied & reduced, r will be neither tidied nor reduced. 
michael@0:  * r should point to different memory than x */
michael@0: void PREFIX(square) (double *r, const double *x) {
michael@0: 	PREFIX(multiply) (r, x, x);
michael@0: }
michael@0: 
michael@0: /* Perform a point doubling in Jacobian coordinates. Input and output
michael@0:  * should be multi-precision floating point integers. */
michael@0: void PREFIX(pt_dbl_jac) (const ecfp_jac_pt * dp, ecfp_jac_pt * dr,
michael@0: 						 const EC_group_fp * group) {
michael@0: 	double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0: 		M[2 * ECFP_NUMDOUBLES], S[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity */
michael@0: 	if (PREFIX(pt_is_inf_jac) (dp) == MP_YES) {
michael@0: 		/* Set r = pt at infinity */
michael@0: 		PREFIX(set_pt_inf_jac) (dr);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* Perform typical point doubling operations */
michael@0: 
michael@0: 	/* TODO? is it worthwhile to do optimizations for when pz = 1? */
michael@0: 
michael@0: 	if (group->aIsM3) {
michael@0: 		/* When a = -3, M = 3(px - pz^2)(px + pz^2) */
michael@0: 		PREFIX(square) (t1, dp->z);
michael@0: 		group->ecfp_reduce(t1, t1, group);	/* 2^23 since the negative
michael@0: 											 * rounding buys another bit */
michael@0: 		PREFIX(addShort) (t0, dp->x, t1);	/* 2*2^23 */
michael@0: 		PREFIX(subtractShort) (t1, dp->x, t1);	/* 2 * 2^23 */
michael@0: 		PREFIX(multiply) (M, t0, t1);	/* 40 * 2^46 */
michael@0: 		PREFIX(addLong) (t0, M, M);	/* 80 * 2^46 */
michael@0: 		PREFIX(addLong) (M, t0, M);	/* 120 * 2^46 < 2^53 */
michael@0: 		group->ecfp_reduce(M, M, group);
michael@0: 	} else {
michael@0: 		/* Generic case */
michael@0: 		/* M = 3 (px^2) + a*(pz^4) */
michael@0: 		PREFIX(square) (t0, dp->x);
michael@0: 		PREFIX(addLong) (M, t0, t0);
michael@0: 		PREFIX(addLong) (t0, t0, M);	/* t0 = 3(px^2) */
michael@0: 		PREFIX(square) (M, dp->z);
michael@0: 		group->ecfp_reduce(M, M, group);
michael@0: 		PREFIX(square) (t1, M);
michael@0: 		group->ecfp_reduce(t1, t1, group);
michael@0: 		PREFIX(multiply) (M, t1, group->curvea);	/* M = a(pz^4) */
michael@0: 		PREFIX(addLong) (M, M, t0);
michael@0: 		group->ecfp_reduce(M, M, group);
michael@0: 	}
michael@0: 
michael@0: 	/* rz = 2 * py * pz */
michael@0: 	PREFIX(multiply) (t1, dp->y, dp->z);
michael@0: 	PREFIX(addLong) (t1, t1, t1);
michael@0: 	group->ecfp_reduce(dr->z, t1, group);
michael@0: 
michael@0: 	/* t0 = 2y^2 */
michael@0: 	PREFIX(square) (t0, dp->y);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(addShort) (t0, t0, t0);
michael@0: 
michael@0: 	/* S = 4 * px * py^2 = 2 * px * t0 */
michael@0: 	PREFIX(multiply) (S, dp->x, t0);
michael@0: 	PREFIX(addLong) (S, S, S);
michael@0: 	group->ecfp_reduce(S, S, group);
michael@0: 
michael@0: 	/* rx = M^2 - 2 * S */
michael@0: 	PREFIX(square) (t1, M);
michael@0: 	PREFIX(subtractShort) (t1, t1, S);
michael@0: 	PREFIX(subtractShort) (t1, t1, S);
michael@0: 	group->ecfp_reduce(dr->x, t1, group);
michael@0: 
michael@0: 	/* ry = M * (S - rx) - 8 * py^4 */
michael@0: 	PREFIX(square) (t1, t0);	/* t1 = 4y^4 */
michael@0: 	PREFIX(subtractShort) (S, S, dr->x);
michael@0: 	PREFIX(multiply) (t0, M, S);
michael@0: 	PREFIX(subtractLong) (t0, t0, t1);
michael@0: 	PREFIX(subtractLong) (t0, t0, t1);
michael@0: 	group->ecfp_reduce(dr->y, t0, group);
michael@0: 
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Perform a point addition using coordinate system Jacobian + Affine ->
michael@0:  * Jacobian. Input and output should be multi-precision floating point
michael@0:  * integers. */
michael@0: void PREFIX(pt_add_jac_aff) (const ecfp_jac_pt * p, const ecfp_aff_pt * q,
michael@0: 							 ecfp_jac_pt * r, const EC_group_fp * group) {
michael@0: 	/* Temporary storage */
michael@0: 	double A[2 * ECFP_NUMDOUBLES], B[2 * ECFP_NUMDOUBLES],
michael@0: 		C[2 * ECFP_NUMDOUBLES], C2[2 * ECFP_NUMDOUBLES],
michael@0: 		D[2 * ECFP_NUMDOUBLES], C3[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity for p or q */
michael@0: 	if (PREFIX(pt_is_inf_aff) (q) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, p->x);
michael@0: 		PREFIX(copy) (r->y, p->y);
michael@0: 		PREFIX(copy) (r->z, p->z);
michael@0: 		goto CLEANUP;
michael@0: 	} else if (PREFIX(pt_is_inf_jac) (p) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, q->x);
michael@0: 		PREFIX(copy) (r->y, q->y);
michael@0: 		/* Since the affine point is not infinity, we can set r->z = 1 */
michael@0: 		PREFIX(one) (r->z);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* Calculates c = qx * pz^2 - px d = (qy * b - py) rx = d^2 - c^3 + 2
michael@0: 	 * (px * c^2) ry = d * (c-rx) - py*c^3 rz = c * pz */
michael@0: 
michael@0: 	/* A = pz^2, B = pz^3 */
michael@0: 	PREFIX(square) (A, p->z);
michael@0: 	group->ecfp_reduce(A, A, group);
michael@0: 	PREFIX(multiply) (B, A, p->z);
michael@0: 	group->ecfp_reduce(B, B, group);
michael@0: 
michael@0: 	/* C = qx * A - px */
michael@0: 	PREFIX(multiply) (C, q->x, A);
michael@0: 	PREFIX(subtractShort) (C, C, p->x);
michael@0: 	group->ecfp_reduce(C, C, group);
michael@0: 
michael@0: 	/* D = qy * B - py */
michael@0: 	PREFIX(multiply) (D, q->y, B);
michael@0: 	PREFIX(subtractShort) (D, D, p->y);
michael@0: 	group->ecfp_reduce(D, D, group);
michael@0: 
michael@0: 	/* C2 = C^2, C3 = C^3 */
michael@0: 	PREFIX(square) (C2, C);
michael@0: 	group->ecfp_reduce(C2, C2, group);
michael@0: 	PREFIX(multiply) (C3, C2, C);
michael@0: 	group->ecfp_reduce(C3, C3, group);
michael@0: 
michael@0: 	/* rz = A = pz * C */
michael@0: 	PREFIX(multiply) (A, p->z, C);
michael@0: 	group->ecfp_reduce(r->z, A, group);
michael@0: 
michael@0: 	/* C = px * C^2, untidied, unreduced */
michael@0: 	PREFIX(multiply) (C, p->x, C2);
michael@0: 
michael@0: 	/* A = D^2, untidied, unreduced */
michael@0: 	PREFIX(square) (A, D);
michael@0: 
michael@0: 	/* rx = B = A - C3 - C - C = D^2 - (C^3 + 2 * (px * C^2) */
michael@0: 	PREFIX(subtractShort) (A, A, C3);
michael@0: 	PREFIX(subtractLong) (A, A, C);
michael@0: 	PREFIX(subtractLong) (A, A, C);
michael@0: 	group->ecfp_reduce(r->x, A, group);
michael@0: 
michael@0: 	/* B = py * C3, untidied, unreduced */
michael@0: 	PREFIX(multiply) (B, p->y, C3);
michael@0: 
michael@0: 	/* C = px * C^2 - rx */
michael@0: 	PREFIX(subtractShort) (C, C, r->x);
michael@0: 	group->ecfp_reduce(C, C, group);
michael@0: 
michael@0: 	/* ry = A = D * C - py * C^3 */
michael@0: 	PREFIX(multiply) (A, D, C);
michael@0: 	PREFIX(subtractLong) (A, A, B);
michael@0: 	group->ecfp_reduce(r->y, A, group);
michael@0: 
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Perform a point addition using Jacobian coordinate system. Input and
michael@0:  * output should be multi-precision floating point integers. */
michael@0: void PREFIX(pt_add_jac) (const ecfp_jac_pt * p, const ecfp_jac_pt * q,
michael@0: 						 ecfp_jac_pt * r, const EC_group_fp * group) {
michael@0: 
michael@0: 	/* Temporary Storage */
michael@0: 	double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0: 		U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
michael@0: 		S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
michael@0: 		H3[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity for p, if so set r = q */
michael@0: 	if (PREFIX(pt_is_inf_jac) (p) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, q->x);
michael@0: 		PREFIX(copy) (r->y, q->y);
michael@0: 		PREFIX(copy) (r->z, q->z);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* Check for point at infinity for p, if so set r = q */
michael@0: 	if (PREFIX(pt_is_inf_jac) (q) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, p->x);
michael@0: 		PREFIX(copy) (r->y, p->y);
michael@0: 		PREFIX(copy) (r->z, p->z);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* U = px * qz^2 , S = py * qz^3 */
michael@0: 	PREFIX(square) (t0, q->z);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (U, p->x, t0);
michael@0: 	group->ecfp_reduce(U, U, group);
michael@0: 	PREFIX(multiply) (t1, t0, q->z);
michael@0: 	group->ecfp_reduce(t1, t1, group);
michael@0: 	PREFIX(multiply) (t0, p->y, t1);
michael@0: 	group->ecfp_reduce(S, t0, group);
michael@0: 
michael@0: 	/* H = qx*(pz)^2 - U , R = (qy * pz^3 - S) */
michael@0: 	PREFIX(square) (t0, p->z);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (H, q->x, t0);
michael@0: 	PREFIX(subtractShort) (H, H, U);
michael@0: 	group->ecfp_reduce(H, H, group);
michael@0: 	PREFIX(multiply) (t1, t0, p->z);	/* t1 = pz^3 */
michael@0: 	group->ecfp_reduce(t1, t1, group);
michael@0: 	PREFIX(multiply) (t0, t1, q->y);	/* t0 = qy * pz^3 */
michael@0: 	PREFIX(subtractShort) (t0, t0, S);
michael@0: 	group->ecfp_reduce(R, t0, group);
michael@0: 
michael@0: 	/* U = U*H^2, H3 = H^3 */
michael@0: 	PREFIX(square) (t0, H);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (t1, U, t0);
michael@0: 	group->ecfp_reduce(U, t1, group);
michael@0: 	PREFIX(multiply) (H3, t0, H);
michael@0: 	group->ecfp_reduce(H3, H3, group);
michael@0: 
michael@0: 	/* rz = pz * qz * H */
michael@0: 	PREFIX(multiply) (t0, q->z, H);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (t1, t0, p->z);
michael@0: 	group->ecfp_reduce(r->z, t1, group);
michael@0: 
michael@0: 	/* rx = R^2 - H^3 - 2 * U */
michael@0: 	PREFIX(square) (t0, R);
michael@0: 	PREFIX(subtractShort) (t0, t0, H3);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	group->ecfp_reduce(r->x, t0, group);
michael@0: 
michael@0: 	/* ry = R(U - rx) - S*H3 */
michael@0: 	PREFIX(subtractShort) (t1, U, r->x);
michael@0: 	PREFIX(multiply) (t0, t1, R);
michael@0: 	PREFIX(multiply) (t1, S, H3);
michael@0: 	PREFIX(subtractLong) (t1, t0, t1);
michael@0: 	group->ecfp_reduce(r->y, t1, group);
michael@0: 
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Perform a point doubling in Modified Jacobian coordinates. Input and
michael@0:  * output should be multi-precision floating point integers. */
michael@0: void PREFIX(pt_dbl_jm) (const ecfp_jm_pt * p, ecfp_jm_pt * r,
michael@0: 						const EC_group_fp * group) {
michael@0: 
michael@0: 	/* Temporary storage */
michael@0: 	double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0: 		M[2 * ECFP_NUMDOUBLES], S[2 * ECFP_NUMDOUBLES],
michael@0: 		U[2 * ECFP_NUMDOUBLES], T[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity */
michael@0: 	if (PREFIX(pt_is_inf_jm) (p) == MP_YES) {
michael@0: 		/* Set r = pt at infinity by setting rz = 0 */
michael@0: 		PREFIX(set_pt_inf_jm) (r);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* M = 3 (px^2) + a*(pz^4) */
michael@0: 	PREFIX(square) (t0, p->x);
michael@0: 	PREFIX(addLong) (M, t0, t0);
michael@0: 	PREFIX(addLong) (t0, t0, M);	/* t0 = 3(px^2) */
michael@0: 	PREFIX(addShort) (t0, t0, p->az4);
michael@0: 	group->ecfp_reduce(M, t0, group);
michael@0: 
michael@0: 	/* rz = 2 * py * pz */
michael@0: 	PREFIX(multiply) (t1, p->y, p->z);
michael@0: 	PREFIX(addLong) (t1, t1, t1);
michael@0: 	group->ecfp_reduce(r->z, t1, group);
michael@0: 
michael@0: 	/* t0 = 2y^2, U = 8y^4 */
michael@0: 	PREFIX(square) (t0, p->y);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(addShort) (t0, t0, t0);
michael@0: 	PREFIX(square) (U, t0);
michael@0: 	group->ecfp_reduce(U, U, group);
michael@0: 	PREFIX(addShort) (U, U, U);
michael@0: 
michael@0: 	/* S = 4 * px * py^2 = 2 * px * t0 */
michael@0: 	PREFIX(multiply) (S, p->x, t0);
michael@0: 	group->ecfp_reduce(S, S, group);
michael@0: 	PREFIX(addShort) (S, S, S);
michael@0: 
michael@0: 	/* rx = M^2 - 2S */
michael@0: 	PREFIX(square) (T, M);
michael@0: 	PREFIX(subtractShort) (T, T, S);
michael@0: 	PREFIX(subtractShort) (T, T, S);
michael@0: 	group->ecfp_reduce(r->x, T, group);
michael@0: 
michael@0: 	/* ry = M * (S - rx) - U */
michael@0: 	PREFIX(subtractShort) (S, S, r->x);
michael@0: 	PREFIX(multiply) (t0, M, S);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	group->ecfp_reduce(r->y, t0, group);
michael@0: 
michael@0: 	/* ra*z^4 = 2*U*(apz4) */
michael@0: 	PREFIX(multiply) (t1, U, p->az4);
michael@0: 	PREFIX(addLong) (t1, t1, t1);
michael@0: 	group->ecfp_reduce(r->az4, t1, group);
michael@0: 
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Perform a point doubling using coordinates Affine -> Chudnovsky
michael@0:  * Jacobian. Input and output should be multi-precision floating point
michael@0:  * integers. */
michael@0: void PREFIX(pt_dbl_aff2chud) (const ecfp_aff_pt * p, ecfp_chud_pt * r,
michael@0: 							  const EC_group_fp * group) {
michael@0: 	double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0: 		M[2 * ECFP_NUMDOUBLES], twoY2[2 * ECFP_NUMDOUBLES],
michael@0: 		S[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity for p, if so set r = O */
michael@0: 	if (PREFIX(pt_is_inf_aff) (p) == MP_YES) {
michael@0: 		PREFIX(set_pt_inf_chud) (r);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* M = 3(px)^2 + a */
michael@0: 	PREFIX(square) (t0, p->x);
michael@0: 	PREFIX(addLong) (t1, t0, t0);
michael@0: 	PREFIX(addLong) (t1, t1, t0);
michael@0: 	PREFIX(addShort) (t1, t1, group->curvea);
michael@0: 	group->ecfp_reduce(M, t1, group);
michael@0: 
michael@0: 	/* twoY2 = 2*(py)^2, S = 4(px)(py)^2 */
michael@0: 	PREFIX(square) (twoY2, p->y);
michael@0: 	PREFIX(addLong) (twoY2, twoY2, twoY2);
michael@0: 	group->ecfp_reduce(twoY2, twoY2, group);
michael@0: 	PREFIX(multiply) (S, p->x, twoY2);
michael@0: 	PREFIX(addLong) (S, S, S);
michael@0: 	group->ecfp_reduce(S, S, group);
michael@0: 
michael@0: 	/* rx = M^2 - 2S */
michael@0: 	PREFIX(square) (t0, M);
michael@0: 	PREFIX(subtractShort) (t0, t0, S);
michael@0: 	PREFIX(subtractShort) (t0, t0, S);
michael@0: 	group->ecfp_reduce(r->x, t0, group);
michael@0: 
michael@0: 	/* ry = M(S-rx) - 8y^4 */
michael@0: 	PREFIX(subtractShort) (t0, S, r->x);
michael@0: 	PREFIX(multiply) (t1, t0, M);
michael@0: 	PREFIX(square) (t0, twoY2);
michael@0: 	PREFIX(subtractLong) (t1, t1, t0);
michael@0: 	PREFIX(subtractLong) (t1, t1, t0);
michael@0: 	group->ecfp_reduce(r->y, t1, group);
michael@0: 
michael@0: 	/* rz = 2py */
michael@0: 	PREFIX(addShort) (r->z, p->y, p->y);
michael@0: 
michael@0: 	/* rz2 = rz^2 */
michael@0: 	PREFIX(square) (t0, r->z);
michael@0: 	group->ecfp_reduce(r->z2, t0, group);
michael@0: 
michael@0: 	/* rz3 = rz^3 */
michael@0: 	PREFIX(multiply) (t0, r->z, r->z2);
michael@0: 	group->ecfp_reduce(r->z3, t0, group);
michael@0: 
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Perform a point addition using coordinates: Modified Jacobian +
michael@0:  * Chudnovsky Jacobian -> Modified Jacobian. Input and output should be
michael@0:  * multi-precision floating point integers. */
michael@0: void PREFIX(pt_add_jm_chud) (ecfp_jm_pt * p, ecfp_chud_pt * q,
michael@0: 							 ecfp_jm_pt * r, const EC_group_fp * group) {
michael@0: 
michael@0: 	double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0: 		U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
michael@0: 		S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
michael@0: 		H3[2 * ECFP_NUMDOUBLES], pz2[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity for p, if so set r = q need to convert
michael@0: 	 * from Chudnovsky form to Modified Jacobian form */
michael@0: 	if (PREFIX(pt_is_inf_jm) (p) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, q->x);
michael@0: 		PREFIX(copy) (r->y, q->y);
michael@0: 		PREFIX(copy) (r->z, q->z);
michael@0: 		PREFIX(square) (t0, q->z2);
michael@0: 		group->ecfp_reduce(t0, t0, group);
michael@0: 		PREFIX(multiply) (t1, t0, group->curvea);
michael@0: 		group->ecfp_reduce(r->az4, t1, group);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 	/* Check for point at infinity for q, if so set r = p */
michael@0: 	if (PREFIX(pt_is_inf_chud) (q) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, p->x);
michael@0: 		PREFIX(copy) (r->y, p->y);
michael@0: 		PREFIX(copy) (r->z, p->z);
michael@0: 		PREFIX(copy) (r->az4, p->az4);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* U = px * qz^2 */
michael@0: 	PREFIX(multiply) (U, p->x, q->z2);
michael@0: 	group->ecfp_reduce(U, U, group);
michael@0: 
michael@0: 	/* H = qx*(pz)^2 - U */
michael@0: 	PREFIX(square) (t0, p->z);
michael@0: 	group->ecfp_reduce(pz2, t0, group);
michael@0: 	PREFIX(multiply) (H, pz2, q->x);
michael@0: 	group->ecfp_reduce(H, H, group);
michael@0: 	PREFIX(subtractShort) (H, H, U);
michael@0: 
michael@0: 	/* U = U*H^2, H3 = H^3 */
michael@0: 	PREFIX(square) (t0, H);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (t1, U, t0);
michael@0: 	group->ecfp_reduce(U, t1, group);
michael@0: 	PREFIX(multiply) (H3, t0, H);
michael@0: 	group->ecfp_reduce(H3, H3, group);
michael@0: 
michael@0: 	/* S = py * qz^3 */
michael@0: 	PREFIX(multiply) (S, p->y, q->z3);
michael@0: 	group->ecfp_reduce(S, S, group);
michael@0: 
michael@0: 	/* R = (qy * z1^3 - s) */
michael@0: 	PREFIX(multiply) (t0, pz2, p->z);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (R, t0, q->y);
michael@0: 	PREFIX(subtractShort) (R, R, S);
michael@0: 	group->ecfp_reduce(R, R, group);
michael@0: 
michael@0: 	/* rz = pz * qz * H */
michael@0: 	PREFIX(multiply) (t1, q->z, H);
michael@0: 	group->ecfp_reduce(t1, t1, group);
michael@0: 	PREFIX(multiply) (t0, p->z, t1);
michael@0: 	group->ecfp_reduce(r->z, t0, group);
michael@0: 
michael@0: 	/* rx = R^2 - H^3 - 2 * U */
michael@0: 	PREFIX(square) (t0, R);
michael@0: 	PREFIX(subtractShort) (t0, t0, H3);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	group->ecfp_reduce(r->x, t0, group);
michael@0: 
michael@0: 	/* ry = R(U - rx) - S*H3 */
michael@0: 	PREFIX(subtractShort) (t1, U, r->x);
michael@0: 	PREFIX(multiply) (t0, t1, R);
michael@0: 	PREFIX(multiply) (t1, S, H3);
michael@0: 	PREFIX(subtractLong) (t1, t0, t1);
michael@0: 	group->ecfp_reduce(r->y, t1, group);
michael@0: 
michael@0: 	if (group->aIsM3) {			/* a == -3 */
michael@0: 		/* a(rz^4) = -3 * ((rz^2)^2) */
michael@0: 		PREFIX(square) (t0, r->z);
michael@0: 		group->ecfp_reduce(t0, t0, group);
michael@0: 		PREFIX(square) (t1, t0);
michael@0: 		PREFIX(addLong) (t0, t1, t1);
michael@0: 		PREFIX(addLong) (t0, t0, t1);
michael@0: 		PREFIX(negLong) (t0, t0);
michael@0: 		group->ecfp_reduce(r->az4, t0, group);
michael@0: 	} else {					/* Generic case */
michael@0: 		/* a(rz^4) = a * ((rz^2)^2) */
michael@0: 		PREFIX(square) (t0, r->z);
michael@0: 		group->ecfp_reduce(t0, t0, group);
michael@0: 		PREFIX(square) (t1, t0);
michael@0: 		group->ecfp_reduce(t1, t1, group);
michael@0: 		PREFIX(multiply) (t0, group->curvea, t1);
michael@0: 		group->ecfp_reduce(r->az4, t0, group);
michael@0: 	}
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Perform a point addition using Chudnovsky Jacobian coordinates. Input
michael@0:  * and output should be multi-precision floating point integers. */
michael@0: void PREFIX(pt_add_chud) (const ecfp_chud_pt * p, const ecfp_chud_pt * q,
michael@0: 						  ecfp_chud_pt * r, const EC_group_fp * group) {
michael@0: 
michael@0: 	/* Temporary Storage */
michael@0: 	double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0: 		U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
michael@0: 		S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
michael@0: 		H3[2 * ECFP_NUMDOUBLES];
michael@0: 
michael@0: 	/* Check for point at infinity for p, if so set r = q */
michael@0: 	if (PREFIX(pt_is_inf_chud) (p) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, q->x);
michael@0: 		PREFIX(copy) (r->y, q->y);
michael@0: 		PREFIX(copy) (r->z, q->z);
michael@0: 		PREFIX(copy) (r->z2, q->z2);
michael@0: 		PREFIX(copy) (r->z3, q->z3);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* Check for point at infinity for p, if so set r = q */
michael@0: 	if (PREFIX(pt_is_inf_chud) (q) == MP_YES) {
michael@0: 		PREFIX(copy) (r->x, p->x);
michael@0: 		PREFIX(copy) (r->y, p->y);
michael@0: 		PREFIX(copy) (r->z, p->z);
michael@0: 		PREFIX(copy) (r->z2, p->z2);
michael@0: 		PREFIX(copy) (r->z3, p->z3);
michael@0: 		goto CLEANUP;
michael@0: 	}
michael@0: 
michael@0: 	/* U = px * qz^2 */
michael@0: 	PREFIX(multiply) (U, p->x, q->z2);
michael@0: 	group->ecfp_reduce(U, U, group);
michael@0: 
michael@0: 	/* H = qx*(pz)^2 - U */
michael@0: 	PREFIX(multiply) (H, q->x, p->z2);
michael@0: 	PREFIX(subtractShort) (H, H, U);
michael@0: 	group->ecfp_reduce(H, H, group);
michael@0: 
michael@0: 	/* U = U*H^2, H3 = H^3 */
michael@0: 	PREFIX(square) (t0, H);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (t1, U, t0);
michael@0: 	group->ecfp_reduce(U, t1, group);
michael@0: 	PREFIX(multiply) (H3, t0, H);
michael@0: 	group->ecfp_reduce(H3, H3, group);
michael@0: 
michael@0: 	/* S = py * qz^3 */
michael@0: 	PREFIX(multiply) (S, p->y, q->z3);
michael@0: 	group->ecfp_reduce(S, S, group);
michael@0: 
michael@0: 	/* rz = pz * qz * H */
michael@0: 	PREFIX(multiply) (t0, q->z, H);
michael@0: 	group->ecfp_reduce(t0, t0, group);
michael@0: 	PREFIX(multiply) (t1, t0, p->z);
michael@0: 	group->ecfp_reduce(r->z, t1, group);
michael@0: 
michael@0: 	/* R = (qy * z1^3 - s) */
michael@0: 	PREFIX(multiply) (t0, q->y, p->z3);
michael@0: 	PREFIX(subtractShort) (t0, t0, S);
michael@0: 	group->ecfp_reduce(R, t0, group);
michael@0: 
michael@0: 	/* rx = R^2 - H^3 - 2 * U */
michael@0: 	PREFIX(square) (t0, R);
michael@0: 	PREFIX(subtractShort) (t0, t0, H3);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	PREFIX(subtractShort) (t0, t0, U);
michael@0: 	group->ecfp_reduce(r->x, t0, group);
michael@0: 
michael@0: 	/* ry = R(U - rx) - S*H3 */
michael@0: 	PREFIX(subtractShort) (t1, U, r->x);
michael@0: 	PREFIX(multiply) (t0, t1, R);
michael@0: 	PREFIX(multiply) (t1, S, H3);
michael@0: 	PREFIX(subtractLong) (t1, t0, t1);
michael@0: 	group->ecfp_reduce(r->y, t1, group);
michael@0: 
michael@0: 	/* rz2 = rz^2 */
michael@0: 	PREFIX(square) (t0, r->z);
michael@0: 	group->ecfp_reduce(r->z2, t0, group);
michael@0: 
michael@0: 	/* rz3 = rz^3 */
michael@0: 	PREFIX(multiply) (t0, r->z, r->z2);
michael@0: 	group->ecfp_reduce(r->z3, t0, group);
michael@0: 
michael@0:   CLEANUP:
michael@0: 	return;
michael@0: }
michael@0: 
michael@0: /* Expects out to be an array of size 16 of Chudnovsky Jacobian points.
michael@0:  * Fills in Chudnovsky Jacobian form (x, y, z, z^2, z^3), for -15P, -13P,
michael@0:  * -11P, -9P, -7P, -5P, -3P, -P, P, 3P, 5P, 7P, 9P, 11P, 13P, 15P */
michael@0: void PREFIX(precompute_chud) (ecfp_chud_pt * out, const ecfp_aff_pt * p,
michael@0: 							  const EC_group_fp * group) {
michael@0: 
michael@0: 	ecfp_chud_pt p2;
michael@0: 
michael@0: 	/* Set out[8] = P */
michael@0: 	PREFIX(copy) (out[8].x, p->x);
michael@0: 	PREFIX(copy) (out[8].y, p->y);
michael@0: 	PREFIX(one) (out[8].z);
michael@0: 	PREFIX(one) (out[8].z2);
michael@0: 	PREFIX(one) (out[8].z3);
michael@0: 
michael@0: 	/* Set p2 = 2P */
michael@0: 	PREFIX(pt_dbl_aff2chud) (p, &p2, group);
michael@0: 
michael@0: 	/* Set 3P, 5P, ..., 15P */
michael@0: 	PREFIX(pt_add_chud) (&out[8], &p2, &out[9], group);
michael@0: 	PREFIX(pt_add_chud) (&out[9], &p2, &out[10], group);
michael@0: 	PREFIX(pt_add_chud) (&out[10], &p2, &out[11], group);
michael@0: 	PREFIX(pt_add_chud) (&out[11], &p2, &out[12], group);
michael@0: 	PREFIX(pt_add_chud) (&out[12], &p2, &out[13], group);
michael@0: 	PREFIX(pt_add_chud) (&out[13], &p2, &out[14], group);
michael@0: 	PREFIX(pt_add_chud) (&out[14], &p2, &out[15], group);
michael@0: 
michael@0: 	/* Set -15P, -13P, ..., -P */
michael@0: 	PREFIX(pt_neg_chud) (&out[8], &out[7]);
michael@0: 	PREFIX(pt_neg_chud) (&out[9], &out[6]);
michael@0: 	PREFIX(pt_neg_chud) (&out[10], &out[5]);
michael@0: 	PREFIX(pt_neg_chud) (&out[11], &out[4]);
michael@0: 	PREFIX(pt_neg_chud) (&out[12], &out[3]);
michael@0: 	PREFIX(pt_neg_chud) (&out[13], &out[2]);
michael@0: 	PREFIX(pt_neg_chud) (&out[14], &out[1]);
michael@0: 	PREFIX(pt_neg_chud) (&out[15], &out[0]);
michael@0: }
michael@0: 
michael@0: /* Expects out to be an array of size 16 of Jacobian points. Fills in
michael@0:  * Jacobian form (x, y, z), for O, P, 2P, ... 15P */
michael@0: void PREFIX(precompute_jac) (ecfp_jac_pt * precomp, const ecfp_aff_pt * p,
michael@0: 							 const EC_group_fp * group) {
michael@0: 	int i;
michael@0: 
michael@0: 	/* fill precomputation table */
michael@0: 	/* set precomp[0] */
michael@0: 	PREFIX(set_pt_inf_jac) (&precomp[0]);
michael@0: 	/* set precomp[1] */
michael@0: 	PREFIX(copy) (precomp[1].x, p->x);
michael@0: 	PREFIX(copy) (precomp[1].y, p->y);
michael@0: 	if (PREFIX(pt_is_inf_aff) (p) == MP_YES) {
michael@0: 		PREFIX(zero) (precomp[1].z);
michael@0: 	} else {
michael@0: 		PREFIX(one) (precomp[1].z);
michael@0: 	}
michael@0: 	/* set precomp[2] */
michael@0: 	group->pt_dbl_jac(&precomp[1], &precomp[2], group);
michael@0: 
michael@0: 	/* set rest of precomp */
michael@0: 	for (i = 3; i < 16; i++) {
michael@0: 		group->pt_add_jac_aff(&precomp[i - 1], p, &precomp[i], group);
michael@0: 	}
michael@0: }