Back to index

lightning-sunbird  0.9+nobinonly
ecp_fpinc.c
Go to the documentation of this file.
00001 /* 
00002  * ***** BEGIN LICENSE BLOCK *****
00003  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00004  *
00005  * The contents of this file are subject to the Mozilla Public License Version
00006  * 1.1 (the "License"); you may not use this file except in compliance with
00007  * the License. You may obtain a copy of the License at
00008  * http://www.mozilla.org/MPL/
00009  *
00010  * Software distributed under the License is distributed on an "AS IS" basis,
00011  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00012  * for the specific language governing rights and limitations under the
00013  * License.
00014  *
00015  * The Original Code is the elliptic curve math library for prime field curves using floating point operations.
00016  *
00017  * The Initial Developer of the Original Code is
00018  * Sun Microsystems, Inc.
00019  * Portions created by the Initial Developer are Copyright (C) 2003
00020  * the Initial Developer. All Rights Reserved.
00021  *
00022  * Contributor(s):
00023  *   Stephen Fung <fungstep@hotmail.com>, Sun Microsystems Laboratories
00024  *
00025  * Alternatively, the contents of this file may be used under the terms of
00026  * either the GNU General Public License Version 2 or later (the "GPL"), or
00027  * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00028  * in which case the provisions of the GPL or the LGPL are applicable instead
00029  * of those above. If you wish to allow use of your version of this file only
00030  * under the terms of either the GPL or the LGPL, and not to allow others to
00031  * use your version of this file under the terms of the MPL, indicate your
00032  * decision by deleting the provisions above and replace them with the notice
00033  * and other provisions required by the GPL or the LGPL. If you do not delete
00034  * the provisions above, a recipient may use your version of this file under
00035  * the terms of any one of the MPL, the GPL or the LGPL.
00036  *
00037  * ***** END LICENSE BLOCK ***** */
00038 
00039 /* This source file is meant to be included by other source files
00040  * (ecp_fp###.c, where ### is one of 160, 192, 224) and should not
00041  * constitute an independent compilation unit. It requires the following
00042  * preprocessor definitions be made: ECFP_BSIZE - the number of bits in
00043  * the field's prime 
00044  * ECFP_NUMDOUBLES - the number of doubles to store one
00045  * multi-precision integer in floating point 
00046 
00047 /* Adds a prefix to a given token to give a unique token name. Prefixes
00048  * with "ecfp" + ECFP_BSIZE + "_". e.g. if ECFP_BSIZE = 160, then
00049  * PREFIX(hello) = ecfp160_hello This optimization allows static function
00050  * linking and compiler loop unrolling without code duplication. */
00051 #ifndef PREFIX
00052 #define PREFIX(b) PREFIX1(ECFP_BSIZE, b)
00053 #define PREFIX1(bsize, b) PREFIX2(bsize, b)
00054 #define PREFIX2(bsize, b) ecfp ## bsize ## _ ## b
00055 #endif
00056 
00057 /* Returns true iff every double in d is 0. (If d == 0 and it is tidied,
00058  * this will be true.) */
00059 mp_err PREFIX(isZero) (const double *d) {
00060        int i;
00061 
00062        for (i = 0; i < ECFP_NUMDOUBLES; i++) {
00063               if (d[i] != 0)
00064                      return MP_NO;
00065        }
00066        return MP_YES;
00067 }
00068 
00069 /* Sets the multi-precision floating point number at t = 0 */
00070 void PREFIX(zero) (double *t) {
00071        int i;
00072 
00073        for (i = 0; i < ECFP_NUMDOUBLES; i++) {
00074               t[i] = 0;
00075        }
00076 }
00077 
00078 /* Sets the multi-precision floating point number at t = 1 */
00079 void PREFIX(one) (double *t) {
00080        int i;
00081 
00082        t[0] = 1;
00083        for (i = 1; i < ECFP_NUMDOUBLES; i++) {
00084               t[i] = 0;
00085        }
00086 }
00087 
00088 /* Checks if point P(x, y, z) is at infinity. Uses Jacobian coordinates. */
00089 mp_err PREFIX(pt_is_inf_jac) (const ecfp_jac_pt * p) {
00090        return PREFIX(isZero) (p->z);
00091 }
00092 
00093 /* Sets the Jacobian point P to be at infinity. */
00094 void PREFIX(set_pt_inf_jac) (ecfp_jac_pt * p) {
00095        PREFIX(zero) (p->z);
00096 }
00097 
00098 /* Checks if point P(x, y) is at infinity. Uses Affine coordinates. */
00099 mp_err PREFIX(pt_is_inf_aff) (const ecfp_aff_pt * p) {
00100        if (PREFIX(isZero) (p->x) == MP_YES && PREFIX(isZero) (p->y) == MP_YES)
00101               return MP_YES;
00102        return MP_NO;
00103 }
00104 
00105 /* Sets the affine point P to be at infinity. */
00106 void PREFIX(set_pt_inf_aff) (ecfp_aff_pt * p) {
00107        PREFIX(zero) (p->x);
00108        PREFIX(zero) (p->y);
00109 }
00110 
00111 /* Checks if point P(x, y, z, a*z^4) is at infinity. Uses Modified
00112  * Jacobian coordinates. */
00113 mp_err PREFIX(pt_is_inf_jm) (const ecfp_jm_pt * p) {
00114        return PREFIX(isZero) (p->z);
00115 }
00116 
00117 /* Sets the Modified Jacobian point P to be at infinity. */
00118 void PREFIX(set_pt_inf_jm) (ecfp_jm_pt * p) {
00119        PREFIX(zero) (p->z);
00120 }
00121 
00122 /* Checks if point P(x, y, z, z^2, z^3) is at infinity. Uses Chudnovsky
00123  * Jacobian coordinates */
00124 mp_err PREFIX(pt_is_inf_chud) (const ecfp_chud_pt * p) {
00125        return PREFIX(isZero) (p->z);
00126 }
00127 
00128 /* Sets the Chudnovsky Jacobian point P to be at infinity. */
00129 void PREFIX(set_pt_inf_chud) (ecfp_chud_pt * p) {
00130        PREFIX(zero) (p->z);
00131 }
00132 
00133 /* Copies a multi-precision floating point number, Setting dest = src */
00134 void PREFIX(copy) (double *dest, const double *src) {
00135        int i;
00136 
00137        for (i = 0; i < ECFP_NUMDOUBLES; i++) {
00138               dest[i] = src[i];
00139        }
00140 }
00141 
00142 /* Sets dest = -src */
00143 void PREFIX(negLong) (double *dest, const double *src) {
00144        int i;
00145 
00146        for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
00147               dest[i] = -src[i];
00148        }
00149 }
00150 
00151 /* Sets r = -p p = (x, y, z, z2, z3) r = (x, -y, z, z2, z3) Uses
00152  * Chudnovsky Jacobian coordinates. */
00153 /* TODO reverse order */
00154 void PREFIX(pt_neg_chud) (const ecfp_chud_pt * p, ecfp_chud_pt * r) {
00155        int i;
00156 
00157        PREFIX(copy) (r->x, p->x);
00158        PREFIX(copy) (r->z, p->z);
00159        PREFIX(copy) (r->z2, p->z2);
00160        PREFIX(copy) (r->z3, p->z3);
00161        for (i = 0; i < ECFP_NUMDOUBLES; i++) {
00162               r->y[i] = -p->y[i];
00163        }
00164 }
00165 
00166 /* Computes r = x + y. Does not tidy or reduce. Any combinations of r, x,
00167  * y can point to the same data. Componentwise adds first ECFP_NUMDOUBLES
00168  * doubles of x and y and stores the result in r. */
00169 void PREFIX(addShort) (double *r, const double *x, const double *y) {
00170        int i;
00171 
00172        for (i = 0; i < ECFP_NUMDOUBLES; i++) {
00173               *r++ = *x++ + *y++;
00174        }
00175 }
00176 
00177 /* Computes r = x + y. Does not tidy or reduce. Any combinations of r, x,
00178  * y can point to the same data. Componentwise adds first
00179  * 2*ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
00180 void PREFIX(addLong) (double *r, const double *x, const double *y) {
00181        int i;
00182 
00183        for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
00184               *r++ = *x++ + *y++;
00185        }
00186 }
00187 
00188 /* Computes r = x - y. Does not tidy or reduce. Any combinations of r, x,
00189  * y can point to the same data. Componentwise subtracts first
00190  * ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
00191 void PREFIX(subtractShort) (double *r, const double *x, const double *y) {
00192        int i;
00193 
00194        for (i = 0; i < ECFP_NUMDOUBLES; i++) {
00195               *r++ = *x++ - *y++;
00196        }
00197 }
00198 
00199 /* Computes r = x - y. Does not tidy or reduce. Any combinations of r, x,
00200  * y can point to the same data. Componentwise subtracts first
00201  * 2*ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
00202 void PREFIX(subtractLong) (double *r, const double *x, const double *y) {
00203        int i;
00204 
00205        for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
00206               *r++ = *x++ - *y++;
00207        }
00208 }
00209 
00210 /* Computes r = x*y.  Both x and y should be tidied and reduced,
00211  * r must be different (point to different memory) than x and y.
00212  * Does not tidy or reduce. */
00213 void PREFIX(multiply)(double *r, const double *x, const double *y) {
00214        int i, j;
00215 
00216        for(j=0;j<ECFP_NUMDOUBLES-1;j++) {
00217               r[j] = x[0] * y[j];
00218               r[j+(ECFP_NUMDOUBLES-1)] = x[ECFP_NUMDOUBLES-1] * y[j];
00219        }
00220        r[ECFP_NUMDOUBLES-1] = x[0] * y[ECFP_NUMDOUBLES-1];
00221        r[ECFP_NUMDOUBLES-1] += x[ECFP_NUMDOUBLES-1] * y[0];
00222        r[2*ECFP_NUMDOUBLES-2] = x[ECFP_NUMDOUBLES-1] * y[ECFP_NUMDOUBLES-1];
00223        r[2*ECFP_NUMDOUBLES-1] = 0;
00224        
00225        for(i=1;i<ECFP_NUMDOUBLES-1;i++) {
00226               for(j=0;j<ECFP_NUMDOUBLES;j++) {
00227                      r[i+j] += (x[i] * y[j]);
00228               }
00229        }
00230 }
00231 
00232 /* Computes the square of x and stores the result in r.  x should be
00233  * tidied & reduced, r will be neither tidied nor reduced. 
00234  * r should point to different memory than x */
00235 void PREFIX(square) (double *r, const double *x) {
00236        PREFIX(multiply) (r, x, x);
00237 }
00238 
00239 /* Perform a point doubling in Jacobian coordinates. Input and output
00240  * should be multi-precision floating point integers. */
00241 void PREFIX(pt_dbl_jac) (const ecfp_jac_pt * dp, ecfp_jac_pt * dr,
00242                                            const EC_group_fp * group) {
00243        double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
00244               M[2 * ECFP_NUMDOUBLES], S[2 * ECFP_NUMDOUBLES];
00245 
00246        /* Check for point at infinity */
00247        if (PREFIX(pt_is_inf_jac) (dp) == MP_YES) {
00248               /* Set r = pt at infinity */
00249               PREFIX(set_pt_inf_jac) (dr);
00250               goto CLEANUP;
00251        }
00252 
00253        /* Perform typical point doubling operations */
00254 
00255        /* TODO? is it worthwhile to do optimizations for when pz = 1? */
00256 
00257        if (group->aIsM3) {
00258               /* When a = -3, M = 3(px - pz^2)(px + pz^2) */
00259               PREFIX(square) (t1, dp->z);
00260               group->ecfp_reduce(t1, t1, group); /* 2^23 since the negative
00261                                                                               * rounding buys another bit */
00262               PREFIX(addShort) (t0, dp->x, t1);  /* 2*2^23 */
00263               PREFIX(subtractShort) (t1, dp->x, t1);    /* 2 * 2^23 */
00264               PREFIX(multiply) (M, t0, t1);      /* 40 * 2^46 */
00265               PREFIX(addLong) (t0, M, M); /* 80 * 2^46 */
00266               PREFIX(addLong) (M, t0, M); /* 120 * 2^46 < 2^53 */
00267               group->ecfp_reduce(M, M, group);
00268        } else {
00269               /* Generic case */
00270               /* M = 3 (px^2) + a*(pz^4) */
00271               PREFIX(square) (t0, dp->x);
00272               PREFIX(addLong) (M, t0, t0);
00273               PREFIX(addLong) (t0, t0, M);       /* t0 = 3(px^2) */
00274               PREFIX(square) (M, dp->z);
00275               group->ecfp_reduce(M, M, group);
00276               PREFIX(square) (t1, M);
00277               group->ecfp_reduce(t1, t1, group);
00278               PREFIX(multiply) (M, t1, group->curvea);  /* M = a(pz^4) */
00279               PREFIX(addLong) (M, M, t0);
00280               group->ecfp_reduce(M, M, group);
00281        }
00282 
00283        /* rz = 2 * py * pz */
00284        PREFIX(multiply) (t1, dp->y, dp->z);
00285        PREFIX(addLong) (t1, t1, t1);
00286        group->ecfp_reduce(dr->z, t1, group);
00287 
00288        /* t0 = 2y^2 */
00289        PREFIX(square) (t0, dp->y);
00290        group->ecfp_reduce(t0, t0, group);
00291        PREFIX(addShort) (t0, t0, t0);
00292 
00293        /* S = 4 * px * py^2 = 2 * px * t0 */
00294        PREFIX(multiply) (S, dp->x, t0);
00295        PREFIX(addLong) (S, S, S);
00296        group->ecfp_reduce(S, S, group);
00297 
00298        /* rx = M^2 - 2 * S */
00299        PREFIX(square) (t1, M);
00300        PREFIX(subtractShort) (t1, t1, S);
00301        PREFIX(subtractShort) (t1, t1, S);
00302        group->ecfp_reduce(dr->x, t1, group);
00303 
00304        /* ry = M * (S - rx) - 8 * py^4 */
00305        PREFIX(square) (t1, t0);    /* t1 = 4y^4 */
00306        PREFIX(subtractShort) (S, S, dr->x);
00307        PREFIX(multiply) (t0, M, S);
00308        PREFIX(subtractLong) (t0, t0, t1);
00309        PREFIX(subtractLong) (t0, t0, t1);
00310        group->ecfp_reduce(dr->y, t0, group);
00311 
00312   CLEANUP:
00313        return;
00314 }
00315 
00316 /* Perform a point addition using coordinate system Jacobian + Affine ->
00317  * Jacobian. Input and output should be multi-precision floating point
00318  * integers. */
00319 void PREFIX(pt_add_jac_aff) (const ecfp_jac_pt * p, const ecfp_aff_pt * q,
00320                                                   ecfp_jac_pt * r, const EC_group_fp * group) {
00321        /* Temporary storage */
00322        double A[2 * ECFP_NUMDOUBLES], B[2 * ECFP_NUMDOUBLES],
00323               C[2 * ECFP_NUMDOUBLES], C2[2 * ECFP_NUMDOUBLES],
00324               D[2 * ECFP_NUMDOUBLES], C3[2 * ECFP_NUMDOUBLES];
00325 
00326        /* Check for point at infinity for p or q */
00327        if (PREFIX(pt_is_inf_aff) (q) == MP_YES) {
00328               PREFIX(copy) (r->x, p->x);
00329               PREFIX(copy) (r->y, p->y);
00330               PREFIX(copy) (r->z, p->z);
00331               goto CLEANUP;
00332        } else if (PREFIX(pt_is_inf_jac) (p) == MP_YES) {
00333               PREFIX(copy) (r->x, q->x);
00334               PREFIX(copy) (r->y, q->y);
00335               /* Since the affine point is not infinity, we can set r->z = 1 */
00336               PREFIX(one) (r->z);
00337               goto CLEANUP;
00338        }
00339 
00340        /* Calculates c = qx * pz^2 - px d = (qy * b - py) rx = d^2 - c^3 + 2
00341         * (px * c^2) ry = d * (c-rx) - py*c^3 rz = c * pz */
00342 
00343        /* A = pz^2, B = pz^3 */
00344        PREFIX(square) (A, p->z);
00345        group->ecfp_reduce(A, A, group);
00346        PREFIX(multiply) (B, A, p->z);
00347        group->ecfp_reduce(B, B, group);
00348 
00349        /* C = qx * A - px */
00350        PREFIX(multiply) (C, q->x, A);
00351        PREFIX(subtractShort) (C, C, p->x);
00352        group->ecfp_reduce(C, C, group);
00353 
00354        /* D = qy * B - py */
00355        PREFIX(multiply) (D, q->y, B);
00356        PREFIX(subtractShort) (D, D, p->y);
00357        group->ecfp_reduce(D, D, group);
00358 
00359        /* C2 = C^2, C3 = C^3 */
00360        PREFIX(square) (C2, C);
00361        group->ecfp_reduce(C2, C2, group);
00362        PREFIX(multiply) (C3, C2, C);
00363        group->ecfp_reduce(C3, C3, group);
00364 
00365        /* rz = A = pz * C */
00366        PREFIX(multiply) (A, p->z, C);
00367        group->ecfp_reduce(r->z, A, group);
00368 
00369        /* C = px * C^2, untidied, unreduced */
00370        PREFIX(multiply) (C, p->x, C2);
00371 
00372        /* A = D^2, untidied, unreduced */
00373        PREFIX(square) (A, D);
00374 
00375        /* rx = B = A - C3 - C - C = D^2 - (C^3 + 2 * (px * C^2) */
00376        PREFIX(subtractShort) (A, A, C3);
00377        PREFIX(subtractLong) (A, A, C);
00378        PREFIX(subtractLong) (A, A, C);
00379        group->ecfp_reduce(r->x, A, group);
00380 
00381        /* B = py * C3, untidied, unreduced */
00382        PREFIX(multiply) (B, p->y, C3);
00383 
00384        /* C = px * C^2 - rx */
00385        PREFIX(subtractShort) (C, C, r->x);
00386        group->ecfp_reduce(C, C, group);
00387 
00388        /* ry = A = D * C - py * C^3 */
00389        PREFIX(multiply) (A, D, C);
00390        PREFIX(subtractLong) (A, A, B);
00391        group->ecfp_reduce(r->y, A, group);
00392 
00393   CLEANUP:
00394        return;
00395 }
00396 
00397 /* Perform a point addition using Jacobian coordinate system. Input and
00398  * output should be multi-precision floating point integers. */
00399 void PREFIX(pt_add_jac) (const ecfp_jac_pt * p, const ecfp_jac_pt * q,
00400                                            ecfp_jac_pt * r, const EC_group_fp * group) {
00401 
00402        /* Temporary Storage */
00403        double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
00404               U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
00405               S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
00406               H3[2 * ECFP_NUMDOUBLES];
00407 
00408        /* Check for point at infinity for p, if so set r = q */
00409        if (PREFIX(pt_is_inf_jac) (p) == MP_YES) {
00410               PREFIX(copy) (r->x, q->x);
00411               PREFIX(copy) (r->y, q->y);
00412               PREFIX(copy) (r->z, q->z);
00413               goto CLEANUP;
00414        }
00415 
00416        /* Check for point at infinity for p, if so set r = q */
00417        if (PREFIX(pt_is_inf_jac) (q) == MP_YES) {
00418               PREFIX(copy) (r->x, p->x);
00419               PREFIX(copy) (r->y, p->y);
00420               PREFIX(copy) (r->z, p->z);
00421               goto CLEANUP;
00422        }
00423 
00424        /* U = px * qz^2 , S = py * qz^3 */
00425        PREFIX(square) (t0, q->z);
00426        group->ecfp_reduce(t0, t0, group);
00427        PREFIX(multiply) (U, p->x, t0);
00428        group->ecfp_reduce(U, U, group);
00429        PREFIX(multiply) (t1, t0, q->z);
00430        group->ecfp_reduce(t1, t1, group);
00431        PREFIX(multiply) (t0, p->y, t1);
00432        group->ecfp_reduce(S, t0, group);
00433 
00434        /* H = qx*(pz)^2 - U , R = (qy * pz^3 - S) */
00435        PREFIX(square) (t0, p->z);
00436        group->ecfp_reduce(t0, t0, group);
00437        PREFIX(multiply) (H, q->x, t0);
00438        PREFIX(subtractShort) (H, H, U);
00439        group->ecfp_reduce(H, H, group);
00440        PREFIX(multiply) (t1, t0, p->z);   /* t1 = pz^3 */
00441        group->ecfp_reduce(t1, t1, group);
00442        PREFIX(multiply) (t0, t1, q->y);   /* t0 = qy * pz^3 */
00443        PREFIX(subtractShort) (t0, t0, S);
00444        group->ecfp_reduce(R, t0, group);
00445 
00446        /* U = U*H^2, H3 = H^3 */
00447        PREFIX(square) (t0, H);
00448        group->ecfp_reduce(t0, t0, group);
00449        PREFIX(multiply) (t1, U, t0);
00450        group->ecfp_reduce(U, t1, group);
00451        PREFIX(multiply) (H3, t0, H);
00452        group->ecfp_reduce(H3, H3, group);
00453 
00454        /* rz = pz * qz * H */
00455        PREFIX(multiply) (t0, q->z, H);
00456        group->ecfp_reduce(t0, t0, group);
00457        PREFIX(multiply) (t1, t0, p->z);
00458        group->ecfp_reduce(r->z, t1, group);
00459 
00460        /* rx = R^2 - H^3 - 2 * U */
00461        PREFIX(square) (t0, R);
00462        PREFIX(subtractShort) (t0, t0, H3);
00463        PREFIX(subtractShort) (t0, t0, U);
00464        PREFIX(subtractShort) (t0, t0, U);
00465        group->ecfp_reduce(r->x, t0, group);
00466 
00467        /* ry = R(U - rx) - S*H3 */
00468        PREFIX(subtractShort) (t1, U, r->x);
00469        PREFIX(multiply) (t0, t1, R);
00470        PREFIX(multiply) (t1, S, H3);
00471        PREFIX(subtractLong) (t1, t0, t1);
00472        group->ecfp_reduce(r->y, t1, group);
00473 
00474   CLEANUP:
00475        return;
00476 }
00477 
00478 /* Perform a point doubling in Modified Jacobian coordinates. Input and
00479  * output should be multi-precision floating point integers. */
00480 void PREFIX(pt_dbl_jm) (const ecfp_jm_pt * p, ecfp_jm_pt * r,
00481                                           const EC_group_fp * group) {
00482 
00483        /* Temporary storage */
00484        double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
00485               M[2 * ECFP_NUMDOUBLES], S[2 * ECFP_NUMDOUBLES],
00486               U[2 * ECFP_NUMDOUBLES], T[2 * ECFP_NUMDOUBLES];
00487 
00488        /* Check for point at infinity */
00489        if (PREFIX(pt_is_inf_jm) (p) == MP_YES) {
00490               /* Set r = pt at infinity by setting rz = 0 */
00491               PREFIX(set_pt_inf_jm) (r);
00492               goto CLEANUP;
00493        }
00494 
00495        /* M = 3 (px^2) + a*(pz^4) */
00496        PREFIX(square) (t0, p->x);
00497        PREFIX(addLong) (M, t0, t0);
00498        PREFIX(addLong) (t0, t0, M);       /* t0 = 3(px^2) */
00499        PREFIX(addShort) (t0, t0, p->az4);
00500        group->ecfp_reduce(M, t0, group);
00501 
00502        /* rz = 2 * py * pz */
00503        PREFIX(multiply) (t1, p->y, p->z);
00504        PREFIX(addLong) (t1, t1, t1);
00505        group->ecfp_reduce(r->z, t1, group);
00506 
00507        /* t0 = 2y^2, U = 8y^4 */
00508        PREFIX(square) (t0, p->y);
00509        group->ecfp_reduce(t0, t0, group);
00510        PREFIX(addShort) (t0, t0, t0);
00511        PREFIX(square) (U, t0);
00512        group->ecfp_reduce(U, U, group);
00513        PREFIX(addShort) (U, U, U);
00514 
00515        /* S = 4 * px * py^2 = 2 * px * t0 */
00516        PREFIX(multiply) (S, p->x, t0);
00517        group->ecfp_reduce(S, S, group);
00518        PREFIX(addShort) (S, S, S);
00519 
00520        /* rx = M^2 - 2S */
00521        PREFIX(square) (T, M);
00522        PREFIX(subtractShort) (T, T, S);
00523        PREFIX(subtractShort) (T, T, S);
00524        group->ecfp_reduce(r->x, T, group);
00525 
00526        /* ry = M * (S - rx) - U */
00527        PREFIX(subtractShort) (S, S, r->x);
00528        PREFIX(multiply) (t0, M, S);
00529        PREFIX(subtractShort) (t0, t0, U);
00530        group->ecfp_reduce(r->y, t0, group);
00531 
00532        /* ra*z^4 = 2*U*(apz4) */
00533        PREFIX(multiply) (t1, U, p->az4);
00534        PREFIX(addLong) (t1, t1, t1);
00535        group->ecfp_reduce(r->az4, t1, group);
00536 
00537   CLEANUP:
00538        return;
00539 }
00540 
00541 /* Perform a point doubling using coordinates Affine -> Chudnovsky
00542  * Jacobian. Input and output should be multi-precision floating point
00543  * integers. */
00544 void PREFIX(pt_dbl_aff2chud) (const ecfp_aff_pt * p, ecfp_chud_pt * r,
00545                                                    const EC_group_fp * group) {
00546        double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
00547               M[2 * ECFP_NUMDOUBLES], twoY2[2 * ECFP_NUMDOUBLES],
00548               S[2 * ECFP_NUMDOUBLES];
00549 
00550        /* Check for point at infinity for p, if so set r = O */
00551        if (PREFIX(pt_is_inf_aff) (p) == MP_YES) {
00552               PREFIX(set_pt_inf_chud) (r);
00553               goto CLEANUP;
00554        }
00555 
00556        /* M = 3(px)^2 + a */
00557        PREFIX(square) (t0, p->x);
00558        PREFIX(addLong) (t1, t0, t0);
00559        PREFIX(addLong) (t1, t1, t0);
00560        PREFIX(addShort) (t1, t1, group->curvea);
00561        group->ecfp_reduce(M, t1, group);
00562 
00563        /* twoY2 = 2*(py)^2, S = 4(px)(py)^2 */
00564        PREFIX(square) (twoY2, p->y);
00565        PREFIX(addLong) (twoY2, twoY2, twoY2);
00566        group->ecfp_reduce(twoY2, twoY2, group);
00567        PREFIX(multiply) (S, p->x, twoY2);
00568        PREFIX(addLong) (S, S, S);
00569        group->ecfp_reduce(S, S, group);
00570 
00571        /* rx = M^2 - 2S */
00572        PREFIX(square) (t0, M);
00573        PREFIX(subtractShort) (t0, t0, S);
00574        PREFIX(subtractShort) (t0, t0, S);
00575        group->ecfp_reduce(r->x, t0, group);
00576 
00577        /* ry = M(S-rx) - 8y^4 */
00578        PREFIX(subtractShort) (t0, S, r->x);
00579        PREFIX(multiply) (t1, t0, M);
00580        PREFIX(square) (t0, twoY2);
00581        PREFIX(subtractLong) (t1, t1, t0);
00582        PREFIX(subtractLong) (t1, t1, t0);
00583        group->ecfp_reduce(r->y, t1, group);
00584 
00585        /* rz = 2py */
00586        PREFIX(addShort) (r->z, p->y, p->y);
00587 
00588        /* rz2 = rz^2 */
00589        PREFIX(square) (t0, r->z);
00590        group->ecfp_reduce(r->z2, t0, group);
00591 
00592        /* rz3 = rz^3 */
00593        PREFIX(multiply) (t0, r->z, r->z2);
00594        group->ecfp_reduce(r->z3, t0, group);
00595 
00596   CLEANUP:
00597        return;
00598 }
00599 
00600 /* Perform a point addition using coordinates: Modified Jacobian +
00601  * Chudnovsky Jacobian -> Modified Jacobian. Input and output should be
00602  * multi-precision floating point integers. */
00603 void PREFIX(pt_add_jm_chud) (ecfp_jm_pt * p, ecfp_chud_pt * q,
00604                                                   ecfp_jm_pt * r, const EC_group_fp * group) {
00605 
00606        double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
00607               U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
00608               S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
00609               H3[2 * ECFP_NUMDOUBLES], pz2[2 * ECFP_NUMDOUBLES];
00610 
00611        /* Check for point at infinity for p, if so set r = q need to convert
00612         * from Chudnovsky form to Modified Jacobian form */
00613        if (PREFIX(pt_is_inf_jm) (p) == MP_YES) {
00614               PREFIX(copy) (r->x, q->x);
00615               PREFIX(copy) (r->y, q->y);
00616               PREFIX(copy) (r->z, q->z);
00617               PREFIX(square) (t0, q->z2);
00618               group->ecfp_reduce(t0, t0, group);
00619               PREFIX(multiply) (t1, t0, group->curvea);
00620               group->ecfp_reduce(r->az4, t1, group);
00621               goto CLEANUP;
00622        }
00623        /* Check for point at infinity for q, if so set r = p */
00624        if (PREFIX(pt_is_inf_chud) (q) == MP_YES) {
00625               PREFIX(copy) (r->x, p->x);
00626               PREFIX(copy) (r->y, p->y);
00627               PREFIX(copy) (r->z, p->z);
00628               PREFIX(copy) (r->az4, p->az4);
00629               goto CLEANUP;
00630        }
00631 
00632        /* U = px * qz^2 */
00633        PREFIX(multiply) (U, p->x, q->z2);
00634        group->ecfp_reduce(U, U, group);
00635 
00636        /* H = qx*(pz)^2 - U */
00637        PREFIX(square) (t0, p->z);
00638        group->ecfp_reduce(pz2, t0, group);
00639        PREFIX(multiply) (H, pz2, q->x);
00640        group->ecfp_reduce(H, H, group);
00641        PREFIX(subtractShort) (H, H, U);
00642 
00643        /* U = U*H^2, H3 = H^3 */
00644        PREFIX(square) (t0, H);
00645        group->ecfp_reduce(t0, t0, group);
00646        PREFIX(multiply) (t1, U, t0);
00647        group->ecfp_reduce(U, t1, group);
00648        PREFIX(multiply) (H3, t0, H);
00649        group->ecfp_reduce(H3, H3, group);
00650 
00651        /* S = py * qz^3 */
00652        PREFIX(multiply) (S, p->y, q->z3);
00653        group->ecfp_reduce(S, S, group);
00654 
00655        /* R = (qy * z1^3 - s) */
00656        PREFIX(multiply) (t0, pz2, p->z);
00657        group->ecfp_reduce(t0, t0, group);
00658        PREFIX(multiply) (R, t0, q->y);
00659        PREFIX(subtractShort) (R, R, S);
00660        group->ecfp_reduce(R, R, group);
00661 
00662        /* rz = pz * qz * H */
00663        PREFIX(multiply) (t1, q->z, H);
00664        group->ecfp_reduce(t1, t1, group);
00665        PREFIX(multiply) (t0, p->z, t1);
00666        group->ecfp_reduce(r->z, t0, group);
00667 
00668        /* rx = R^2 - H^3 - 2 * U */
00669        PREFIX(square) (t0, R);
00670        PREFIX(subtractShort) (t0, t0, H3);
00671        PREFIX(subtractShort) (t0, t0, U);
00672        PREFIX(subtractShort) (t0, t0, U);
00673        group->ecfp_reduce(r->x, t0, group);
00674 
00675        /* ry = R(U - rx) - S*H3 */
00676        PREFIX(subtractShort) (t1, U, r->x);
00677        PREFIX(multiply) (t0, t1, R);
00678        PREFIX(multiply) (t1, S, H3);
00679        PREFIX(subtractLong) (t1, t0, t1);
00680        group->ecfp_reduce(r->y, t1, group);
00681 
00682        if (group->aIsM3) {                /* a == -3 */
00683               /* a(rz^4) = -3 * ((rz^2)^2) */
00684               PREFIX(square) (t0, r->z);
00685               group->ecfp_reduce(t0, t0, group);
00686               PREFIX(square) (t1, t0);
00687               PREFIX(addLong) (t0, t1, t1);
00688               PREFIX(addLong) (t0, t0, t1);
00689               PREFIX(negLong) (t0, t0);
00690               group->ecfp_reduce(r->az4, t0, group);
00691        } else {                                  /* Generic case */
00692               /* a(rz^4) = a * ((rz^2)^2) */
00693               PREFIX(square) (t0, r->z);
00694               group->ecfp_reduce(t0, t0, group);
00695               PREFIX(square) (t1, t0);
00696               group->ecfp_reduce(t1, t1, group);
00697               PREFIX(multiply) (t0, group->curvea, t1);
00698               group->ecfp_reduce(r->az4, t0, group);
00699        }
00700   CLEANUP:
00701        return;
00702 }
00703 
00704 /* Perform a point addition using Chudnovsky Jacobian coordinates. Input
00705  * and output should be multi-precision floating point integers. */
00706 void PREFIX(pt_add_chud) (const ecfp_chud_pt * p, const ecfp_chud_pt * q,
00707                                             ecfp_chud_pt * r, const EC_group_fp * group) {
00708 
00709        /* Temporary Storage */
00710        double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
00711               U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
00712               S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
00713               H3[2 * ECFP_NUMDOUBLES];
00714 
00715        /* Check for point at infinity for p, if so set r = q */
00716        if (PREFIX(pt_is_inf_chud) (p) == MP_YES) {
00717               PREFIX(copy) (r->x, q->x);
00718               PREFIX(copy) (r->y, q->y);
00719               PREFIX(copy) (r->z, q->z);
00720               PREFIX(copy) (r->z2, q->z2);
00721               PREFIX(copy) (r->z3, q->z3);
00722               goto CLEANUP;
00723        }
00724 
00725        /* Check for point at infinity for p, if so set r = q */
00726        if (PREFIX(pt_is_inf_chud) (q) == MP_YES) {
00727               PREFIX(copy) (r->x, p->x);
00728               PREFIX(copy) (r->y, p->y);
00729               PREFIX(copy) (r->z, p->z);
00730               PREFIX(copy) (r->z2, p->z2);
00731               PREFIX(copy) (r->z3, p->z3);
00732               goto CLEANUP;
00733        }
00734 
00735        /* U = px * qz^2 */
00736        PREFIX(multiply) (U, p->x, q->z2);
00737        group->ecfp_reduce(U, U, group);
00738 
00739        /* H = qx*(pz)^2 - U */
00740        PREFIX(multiply) (H, q->x, p->z2);
00741        PREFIX(subtractShort) (H, H, U);
00742        group->ecfp_reduce(H, H, group);
00743 
00744        /* U = U*H^2, H3 = H^3 */
00745        PREFIX(square) (t0, H);
00746        group->ecfp_reduce(t0, t0, group);
00747        PREFIX(multiply) (t1, U, t0);
00748        group->ecfp_reduce(U, t1, group);
00749        PREFIX(multiply) (H3, t0, H);
00750        group->ecfp_reduce(H3, H3, group);
00751 
00752        /* S = py * qz^3 */
00753        PREFIX(multiply) (S, p->y, q->z3);
00754        group->ecfp_reduce(S, S, group);
00755 
00756        /* rz = pz * qz * H */
00757        PREFIX(multiply) (t0, q->z, H);
00758        group->ecfp_reduce(t0, t0, group);
00759        PREFIX(multiply) (t1, t0, p->z);
00760        group->ecfp_reduce(r->z, t1, group);
00761 
00762        /* R = (qy * z1^3 - s) */
00763        PREFIX(multiply) (t0, q->y, p->z3);
00764        PREFIX(subtractShort) (t0, t0, S);
00765        group->ecfp_reduce(R, t0, group);
00766 
00767        /* rx = R^2 - H^3 - 2 * U */
00768        PREFIX(square) (t0, R);
00769        PREFIX(subtractShort) (t0, t0, H3);
00770        PREFIX(subtractShort) (t0, t0, U);
00771        PREFIX(subtractShort) (t0, t0, U);
00772        group->ecfp_reduce(r->x, t0, group);
00773 
00774        /* ry = R(U - rx) - S*H3 */
00775        PREFIX(subtractShort) (t1, U, r->x);
00776        PREFIX(multiply) (t0, t1, R);
00777        PREFIX(multiply) (t1, S, H3);
00778        PREFIX(subtractLong) (t1, t0, t1);
00779        group->ecfp_reduce(r->y, t1, group);
00780 
00781        /* rz2 = rz^2 */
00782        PREFIX(square) (t0, r->z);
00783        group->ecfp_reduce(r->z2, t0, group);
00784 
00785        /* rz3 = rz^3 */
00786        PREFIX(multiply) (t0, r->z, r->z2);
00787        group->ecfp_reduce(r->z3, t0, group);
00788 
00789   CLEANUP:
00790        return;
00791 }
00792 
00793 /* Expects out to be an array of size 16 of Chudnovsky Jacobian points.
00794  * Fills in Chudnovsky Jacobian form (x, y, z, z^2, z^3), for -15P, -13P,
00795  * -11P, -9P, -7P, -5P, -3P, -P, P, 3P, 5P, 7P, 9P, 11P, 13P, 15P */
00796 void PREFIX(precompute_chud) (ecfp_chud_pt * out, const ecfp_aff_pt * p,
00797                                                    const EC_group_fp * group) {
00798 
00799        ecfp_chud_pt p2;
00800 
00801        /* Set out[8] = P */
00802        PREFIX(copy) (out[8].x, p->x);
00803        PREFIX(copy) (out[8].y, p->y);
00804        PREFIX(one) (out[8].z);
00805        PREFIX(one) (out[8].z2);
00806        PREFIX(one) (out[8].z3);
00807 
00808        /* Set p2 = 2P */
00809        PREFIX(pt_dbl_aff2chud) (p, &p2, group);
00810 
00811        /* Set 3P, 5P, ..., 15P */
00812        PREFIX(pt_add_chud) (&out[8], &p2, &out[9], group);
00813        PREFIX(pt_add_chud) (&out[9], &p2, &out[10], group);
00814        PREFIX(pt_add_chud) (&out[10], &p2, &out[11], group);
00815        PREFIX(pt_add_chud) (&out[11], &p2, &out[12], group);
00816        PREFIX(pt_add_chud) (&out[12], &p2, &out[13], group);
00817        PREFIX(pt_add_chud) (&out[13], &p2, &out[14], group);
00818        PREFIX(pt_add_chud) (&out[14], &p2, &out[15], group);
00819 
00820        /* Set -15P, -13P, ..., -P */
00821        PREFIX(pt_neg_chud) (&out[8], &out[7]);
00822        PREFIX(pt_neg_chud) (&out[9], &out[6]);
00823        PREFIX(pt_neg_chud) (&out[10], &out[5]);
00824        PREFIX(pt_neg_chud) (&out[11], &out[4]);
00825        PREFIX(pt_neg_chud) (&out[12], &out[3]);
00826        PREFIX(pt_neg_chud) (&out[13], &out[2]);
00827        PREFIX(pt_neg_chud) (&out[14], &out[1]);
00828        PREFIX(pt_neg_chud) (&out[15], &out[0]);
00829 }
00830 
00831 /* Expects out to be an array of size 16 of Jacobian points. Fills in
00832  * Jacobian form (x, y, z), for O, P, 2P, ... 15P */
00833 void PREFIX(precompute_jac) (ecfp_jac_pt * precomp, const ecfp_aff_pt * p,
00834                                                   const EC_group_fp * group) {
00835        int i;
00836 
00837        /* fill precomputation table */
00838        /* set precomp[0] */
00839        PREFIX(set_pt_inf_jac) (&precomp[0]);
00840        /* set precomp[1] */
00841        PREFIX(copy) (precomp[1].x, p->x);
00842        PREFIX(copy) (precomp[1].y, p->y);
00843        if (PREFIX(pt_is_inf_aff) (p) == MP_YES) {
00844               PREFIX(zero) (precomp[1].z);
00845        } else {
00846               PREFIX(one) (precomp[1].z);
00847        }
00848        /* set precomp[2] */
00849        group->pt_dbl_jac(&precomp[1], &precomp[2], group);
00850 
00851        /* set rest of precomp */
00852        for (i = 3; i < 16; i++) {
00853               group->pt_add_jac_aff(&precomp[i - 1], p, &precomp[i], group);
00854        }
00855 }