Back to index

lightning-sunbird  0.9+nobinonly
ecp_fp.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
00016  * using floating point operations.
00017  *
00018  * The Initial Developer of the Original Code is
00019  * Sun Microsystems, Inc.
00020  * Portions created by the Initial Developer are Copyright (C) 2003
00021  * the Initial Developer. All Rights Reserved.
00022  *
00023  * Contributor(s):
00024  *   Sheueling Chang-Shantz <sheueling.chang@sun.com>,
00025  *   Stephen Fung <fungstep@hotmail.com>, and
00026  *   Douglas Stebila <douglas@stebila.ca>, Sun Microsystems Laboratories.
00027  *
00028  * Alternatively, the contents of this file may be used under the terms of
00029  * either the GNU General Public License Version 2 or later (the "GPL"), or
00030  * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00031  * in which case the provisions of the GPL or the LGPL are applicable instead
00032  * of those above. If you wish to allow use of your version of this file only
00033  * under the terms of either the GPL or the LGPL, and not to allow others to
00034  * use your version of this file under the terms of the MPL, indicate your
00035  * decision by deleting the provisions above and replace them with the notice
00036  * and other provisions required by the GPL or the LGPL. If you do not delete
00037  * the provisions above, a recipient may use your version of this file under
00038  * the terms of any one of the MPL, the GPL or the LGPL.
00039  *
00040  * ***** END LICENSE BLOCK ***** */
00041 
00042 #include "ecp_fp.h"
00043 #include "ecl-priv.h"
00044 #include <stdlib.h>
00045 
00046 /* Performs tidying on a short multi-precision floating point integer (the 
00047  * lower group->numDoubles floats). */
00048 void
00049 ecfp_tidyShort(double *t, const EC_group_fp * group)
00050 {
00051        group->ecfp_tidy(t, group->alpha, group);
00052 }
00053 
00054 /* Performs tidying on only the upper float digits of a multi-precision
00055  * floating point integer, i.e. the digits beyond the regular length which 
00056  * are removed in the reduction step. */
00057 void
00058 ecfp_tidyUpper(double *t, const EC_group_fp * group)
00059 {
00060        group->ecfp_tidy(t + group->numDoubles,
00061                                     group->alpha + group->numDoubles, group);
00062 }
00063 
00064 /* Performs a "tidy" operation, which performs carrying, moving excess
00065  * bits from one double to the next double, so that the precision of the
00066  * doubles is reduced to the regular precision group->doubleBitSize. This
00067  * might result in some float digits being negative. Alternative C version 
00068  * for portability. */
00069 void
00070 ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group)
00071 {
00072        double q;
00073        int i;
00074 
00075        /* Do carrying */
00076        for (i = 0; i < group->numDoubles - 1; i++) {
00077               q = t[i] + alpha[i + 1];
00078               q -= alpha[i + 1];
00079               t[i] -= q;
00080               t[i + 1] += q;
00081 
00082               /* If we don't assume that truncation rounding is used, then q
00083                * might be 2^n bigger than expected (if it rounds up), then t[0]
00084                * could be negative and t[1] 2^n larger than expected. */
00085        }
00086 }
00087 
00088 /* Performs a more mathematically precise "tidying" so that each term is
00089  * positive.  This is slower than the regular tidying, and is used for
00090  * conversion from floating point to integer. */
00091 void
00092 ecfp_positiveTidy(double *t, const EC_group_fp * group)
00093 {
00094        double q;
00095        int i;
00096 
00097        /* Do carrying */
00098        for (i = 0; i < group->numDoubles - 1; i++) {
00099               /* Subtract beta to force rounding down */
00100               q = t[i] - ecfp_beta[i + 1];
00101               q += group->alpha[i + 1];
00102               q -= group->alpha[i + 1];
00103               t[i] -= q;
00104               t[i + 1] += q;
00105 
00106               /* Due to subtracting ecfp_beta, we should have each term a
00107                * non-negative int */
00108               ECFP_ASSERT(t[i] / ecfp_exp[i] ==
00109                                    (unsigned long long) (t[i] / ecfp_exp[i]));
00110               ECFP_ASSERT(t[i] >= 0);
00111        }
00112 }
00113 
00114 /* Converts from a floating point representation into an mp_int. Expects
00115  * that d is already reduced. */
00116 void
00117 ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup)
00118 {
00119        EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
00120        unsigned short i16[(group->primeBitSize + 15) / 16];
00121        double q = 1;
00122 
00123 #ifdef ECL_THIRTY_TWO_BIT
00124        /* TEST uint32_t z = 0; */
00125        unsigned int z = 0;
00126 #else
00127        uint64_t z = 0;
00128 #endif
00129        int zBits = 0;
00130        int copiedBits = 0;
00131        int i = 0;
00132        int j = 0;
00133 
00134        mp_digit *out;
00135 
00136        /* Result should always be >= 0, so set sign accordingly */
00137        MP_SIGN(mpout) = MP_ZPOS;
00138 
00139        /* Tidy up so we're just dealing with positive numbers */
00140        ecfp_positiveTidy(d, group);
00141 
00142        /* We might need to do this reduction step more than once if the
00143         * reduction adds smaller terms which carry-over to cause another
00144         * reduction. However, this should happen very rarely, if ever,
00145         * depending on the elliptic curve. */
00146        do {
00147               /* Init loop data */
00148               z = 0;
00149               zBits = 0;
00150               q = 1;
00151               i = 0;
00152               j = 0;
00153               copiedBits = 0;
00154 
00155               /* Might have to do a bit more reduction */
00156               group->ecfp_singleReduce(d, group);
00157 
00158               /* Grow the size of the mpint if it's too small */
00159               s_mp_grow(mpout, group->numInts);
00160               MP_USED(mpout) = group->numInts;
00161               out = MP_DIGITS(mpout);
00162 
00163               /* Convert double to 16 bit integers */
00164               while (copiedBits < group->primeBitSize) {
00165                      if (zBits < 16) {
00166                             z += d[i] * q;
00167                             i++;
00168                             ECFP_ASSERT(i < (group->primeBitSize + 15) / 16);
00169                             zBits += group->doubleBitSize;
00170                      }
00171                      i16[j] = z;
00172                      j++;
00173                      z >>= 16;
00174                      zBits -= 16;
00175                      q *= ecfp_twom16;
00176                      copiedBits += 16;
00177               }
00178        } while (z != 0);
00179 
00180        /* Convert 16 bit integers to mp_digit */
00181 #ifdef ECL_THIRTY_TWO_BIT
00182        for (i = 0; i < (group->primeBitSize + 15) / 16; i += 2) {
00183               *out = 0;
00184               if (i + 1 < (group->primeBitSize + 15) / 16) {
00185                      *out = i16[i + 1];
00186                      *out <<= 16;
00187               }
00188               *out++ += i16[i];
00189        }
00190 #else                                            /* 64 bit */
00191        for (i = 0; i < (group->primeBitSize + 15) / 16; i += 4) {
00192               *out = 0;
00193               if (i + 3 < (group->primeBitSize + 15) / 16) {
00194                      *out = i16[i + 3];
00195                      *out <<= 16;
00196               }
00197               if (i + 2 < (group->primeBitSize + 15) / 16) {
00198                      *out += i16[i + 2];
00199                      *out <<= 16;
00200               }
00201               if (i + 1 < (group->primeBitSize + 15) / 16) {
00202                      *out += i16[i + 1];
00203                      *out <<= 16;
00204               }
00205               *out++ += i16[i];
00206        }
00207 #endif
00208 
00209        /* Perform final reduction.  mpout should already be the same number
00210         * of bits as p, but might not be less than p.  Make it so. Since
00211         * mpout has the same number of bits as p, and 2p has a larger bit
00212         * size, then mpout < 2p, so a single subtraction of p will suffice. */
00213        if (mp_cmp(mpout, &ecgroup->meth->irr) >= 0) {
00214               mp_sub(mpout, &ecgroup->meth->irr, mpout);
00215        }
00216 
00217        /* Shrink the size of the mp_int to the actual used size (required for 
00218         * mp_cmp_z == 0) */
00219        out = MP_DIGITS(mpout);
00220        for (i = group->numInts - 1; i > 0; i--) {
00221               if (out[i] != 0)
00222                      break;
00223        }
00224        MP_USED(mpout) = i + 1;
00225 
00226        /* Should be between 0 and p-1 */
00227        ECFP_ASSERT(mp_cmp(mpout, &ecgroup->meth->irr) < 0);
00228        ECFP_ASSERT(mp_cmp_z(mpout) >= 0);
00229 }
00230 
00231 /* Converts from an mpint into a floating point representation. */
00232 void
00233 ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup)
00234 {
00235        int i;
00236        int j = 0;
00237        int size;
00238        double shift = 1;
00239        mp_digit *in;
00240        EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
00241 
00242 #ifdef ECL_DEBUG
00243        /* if debug mode, convert result back using ecfp_fp2i into cmp, then
00244         * compare to x. */
00245        mp_int cmp;
00246 
00247        MP_DIGITS(&cmp) = NULL;
00248        mp_init(&cmp);
00249 #endif
00250 
00251        ECFP_ASSERT(group != NULL);
00252 
00253        /* init output to 0 (since we skip over some terms) */
00254        for (i = 0; i < group->numDoubles; i++)
00255               out[i] = 0;
00256        i = 0;
00257 
00258        size = MP_USED(x);
00259        in = MP_DIGITS(x);
00260 
00261        /* Copy from int into doubles */
00262 #ifdef ECL_THIRTY_TWO_BIT
00263        while (j < size) {
00264               while (group->doubleBitSize * (i + 1) <= 32 * j) {
00265                      i++;
00266               }
00267               ECFP_ASSERT(group->doubleBitSize * i <= 32 * j);
00268               out[i] = in[j];
00269               out[i] *= shift;
00270               shift *= ecfp_two32;
00271               j++;
00272        }
00273 #else
00274        while (j < size) {
00275               while (group->doubleBitSize * (i + 1) <= 64 * j) {
00276                      i++;
00277               }
00278               ECFP_ASSERT(group->doubleBitSize * i <= 64 * j);
00279               out[i] = (in[j] & 0x00000000FFFFFFFF) * shift;
00280 
00281               while (group->doubleBitSize * (i + 1) <= 64 * j + 32) {
00282                      i++;
00283               }
00284               ECFP_ASSERT(24 * i <= 64 * j + 32);
00285               out[i] = (in[j] & 0xFFFFFFFF00000000) * shift;
00286 
00287               shift *= ecfp_two64;
00288               j++;
00289        }
00290 #endif
00291        /* Realign bits to match double boundaries */
00292        ecfp_tidyShort(out, group);
00293 
00294 #ifdef ECL_DEBUG
00295        /* Convert result back to mp_int, compare to original */
00296        ecfp_fp2i(&cmp, out, ecgroup);
00297        ECFP_ASSERT(mp_cmp(&cmp, x) == 0);
00298        mp_clear(&cmp);
00299 #endif
00300 }
00301 
00302 /* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters
00303  * a, b and p are the elliptic curve coefficients and the prime that
00304  * determines the field GFp.  Elliptic curve points P and R can be
00305  * identical.  Uses Jacobian coordinates. Uses 4-bit window method. */
00306 mp_err
00307 ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px,
00308                                              const mp_int *py, mp_int *rx, mp_int *ry,
00309                                              const ECGroup *ecgroup)
00310 {
00311        mp_err res = MP_OKAY;
00312        ecfp_jac_pt precomp[16], r;
00313        ecfp_aff_pt p;
00314        EC_group_fp *group;
00315 
00316        mp_int rz;
00317        int i, ni, d;
00318 
00319        ARGCHK(ecgroup != NULL, MP_BADARG);
00320        ARGCHK((n != NULL) && (px != NULL) && (py != NULL), MP_BADARG);
00321 
00322        group = (EC_group_fp *) ecgroup->extra1;
00323        MP_DIGITS(&rz) = 0;
00324        MP_CHECKOK(mp_init(&rz));
00325 
00326        /* init p, da */
00327        ecfp_i2fp(p.x, px, ecgroup);
00328        ecfp_i2fp(p.y, py, ecgroup);
00329        ecfp_i2fp(group->curvea, &ecgroup->curvea, ecgroup);
00330 
00331        /* Do precomputation */
00332        group->precompute_jac(precomp, &p, group);
00333 
00334        /* Do main body of calculations */
00335        d = (mpl_significant_bits(n) + 3) / 4;
00336 
00337        /* R = inf */
00338        for (i = 0; i < group->numDoubles; i++) {
00339               r.z[i] = 0;
00340        }
00341 
00342        for (i = d - 1; i >= 0; i--) {
00343               /* compute window ni */
00344               ni = MP_GET_BIT(n, 4 * i + 3);
00345               ni <<= 1;
00346               ni |= MP_GET_BIT(n, 4 * i + 2);
00347               ni <<= 1;
00348               ni |= MP_GET_BIT(n, 4 * i + 1);
00349               ni <<= 1;
00350               ni |= MP_GET_BIT(n, 4 * i);
00351 
00352               /* R = 2^4 * R */
00353               group->pt_dbl_jac(&r, &r, group);
00354               group->pt_dbl_jac(&r, &r, group);
00355               group->pt_dbl_jac(&r, &r, group);
00356               group->pt_dbl_jac(&r, &r, group);
00357 
00358               /* R = R + (ni * P) */
00359               group->pt_add_jac(&r, &precomp[ni], &r, group);
00360        }
00361 
00362        /* Convert back to integer */
00363        ecfp_fp2i(rx, r.x, ecgroup);
00364        ecfp_fp2i(ry, r.y, ecgroup);
00365        ecfp_fp2i(&rz, r.z, ecgroup);
00366 
00367        /* convert result S to affine coordinates */
00368        MP_CHECKOK(ec_GFp_pt_jac2aff(rx, ry, &rz, rx, ry, ecgroup));
00369 
00370   CLEANUP:
00371        mp_clear(&rz);
00372        return res;
00373 }
00374 
00375 /* Uses mixed Jacobian-affine coordinates to perform a point
00376  * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine
00377  * coordinates (Jacobian coordinates for doubles and affine coordinates
00378  * for additions; based on recommendation from Brown et al.). Not very
00379  * time efficient but quite space efficient, no precomputation needed.
00380  * group contains the elliptic curve coefficients and the prime that
00381  * determines the field GFp.  Elliptic curve points P and R can be
00382  * identical. Performs calculations in floating point number format, since 
00383  * this is faster than the integer operations on the ULTRASPARC III.
00384  * Uses left-to-right binary method (double & add) (algorithm 9) for
00385  * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes.
00386  * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */
00387 mp_err
00388 ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py,
00389                                     mp_int *rx, mp_int *ry, const ECGroup *ecgroup)
00390 {
00391        mp_err res;
00392        mp_int sx, sy, sz;
00393 
00394        ecfp_aff_pt p;
00395        ecfp_jac_pt r;
00396        EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
00397 
00398        int i, l;
00399 
00400        MP_DIGITS(&sx) = 0;
00401        MP_DIGITS(&sy) = 0;
00402        MP_DIGITS(&sz) = 0;
00403        MP_CHECKOK(mp_init(&sx));
00404        MP_CHECKOK(mp_init(&sy));
00405        MP_CHECKOK(mp_init(&sz));
00406 
00407        /* if n = 0 then r = inf */
00408        if (mp_cmp_z(n) == 0) {
00409               mp_zero(rx);
00410               mp_zero(ry);
00411               res = MP_OKAY;
00412               goto CLEANUP;
00413               /* if n < 0 then out of range error */
00414        } else if (mp_cmp_z(n) < 0) {
00415               res = MP_RANGE;
00416               goto CLEANUP;
00417        }
00418 
00419        /* Convert from integer to floating point */
00420        ecfp_i2fp(p.x, px, ecgroup);
00421        ecfp_i2fp(p.y, py, ecgroup);
00422        ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
00423 
00424        /* Init r to point at infinity */
00425        for (i = 0; i < group->numDoubles; i++) {
00426               r.z[i] = 0;
00427        }
00428 
00429        /* double and add method */
00430        l = mpl_significant_bits(n) - 1;
00431 
00432        for (i = l; i >= 0; i--) {
00433               /* R = 2R */
00434               group->pt_dbl_jac(&r, &r, group);
00435 
00436               /* if n_i = 1, then R = R + Q */
00437               if (MP_GET_BIT(n, i) != 0) {
00438                      group->pt_add_jac_aff(&r, &p, &r, group);
00439               }
00440        }
00441 
00442        /* Convert from floating point to integer */
00443        ecfp_fp2i(&sx, r.x, ecgroup);
00444        ecfp_fp2i(&sy, r.y, ecgroup);
00445        ecfp_fp2i(&sz, r.z, ecgroup);
00446 
00447        /* convert result R to affine coordinates */
00448        MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
00449 
00450   CLEANUP:
00451        mp_clear(&sx);
00452        mp_clear(&sy);
00453        mp_clear(&sz);
00454        return res;
00455 }
00456 
00457 /* Computes R = nP where R is (rx, ry) and P is the base point. Elliptic
00458  * curve points P and R can be identical. Uses mixed Modified-Jacobian
00459  * co-ordinates for doubling and Chudnovsky Jacobian coordinates for
00460  * additions. Uses 5-bit window NAF method (algorithm 11) for scalar-point 
00461  * multiplication from Brown, Hankerson, Lopez, Menezes. Software
00462  * Implementation of the NIST Elliptic Curves Over Prime Fields. */
00463 mp_err
00464 ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px,
00465                                            const mp_int *py, mp_int *rx, mp_int *ry,
00466                                            const ECGroup *ecgroup)
00467 {
00468        mp_err res = MP_OKAY;
00469        mp_int sx, sy, sz;
00470        EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
00471        ecfp_chud_pt precomp[16];
00472 
00473        ecfp_aff_pt p;
00474        ecfp_jm_pt r;
00475 
00476        signed char naf[group->orderBitSize + 1];
00477        int i;
00478 
00479        MP_DIGITS(&sx) = 0;
00480        MP_DIGITS(&sy) = 0;
00481        MP_DIGITS(&sz) = 0;
00482        MP_CHECKOK(mp_init(&sx));
00483        MP_CHECKOK(mp_init(&sy));
00484        MP_CHECKOK(mp_init(&sz));
00485 
00486        /* if n = 0 then r = inf */
00487        if (mp_cmp_z(n) == 0) {
00488               mp_zero(rx);
00489               mp_zero(ry);
00490               res = MP_OKAY;
00491               goto CLEANUP;
00492               /* if n < 0 then out of range error */
00493        } else if (mp_cmp_z(n) < 0) {
00494               res = MP_RANGE;
00495               goto CLEANUP;
00496        }
00497 
00498        /* Convert from integer to floating point */
00499        ecfp_i2fp(p.x, px, ecgroup);
00500        ecfp_i2fp(p.y, py, ecgroup);
00501        ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
00502 
00503        /* Perform precomputation */
00504        group->precompute_chud(precomp, &p, group);
00505 
00506        /* Compute 5NAF */
00507        ec_compute_wNAF(naf, group->orderBitSize, n, 5);
00508 
00509        /* Init R = pt at infinity */
00510        for (i = 0; i < group->numDoubles; i++) {
00511               r.z[i] = 0;
00512        }
00513 
00514        /* wNAF method */
00515        for (i = group->orderBitSize; i >= 0; i--) {
00516               /* R = 2R */
00517               group->pt_dbl_jm(&r, &r, group);
00518 
00519               if (naf[i] != 0) {
00520                      group->pt_add_jm_chud(&r, &precomp[(naf[i] + 15) / 2], &r,
00521                                                           group);
00522               }
00523        }
00524 
00525        /* Convert from floating point to integer */
00526        ecfp_fp2i(&sx, r.x, ecgroup);
00527        ecfp_fp2i(&sy, r.y, ecgroup);
00528        ecfp_fp2i(&sz, r.z, ecgroup);
00529 
00530        /* convert result R to affine coordinates */
00531        MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
00532 
00533   CLEANUP:
00534        mp_clear(&sx);
00535        mp_clear(&sy);
00536        mp_clear(&sz);
00537        return res;
00538 }
00539 
00540 /* Cleans up extra memory allocated in ECGroup for this implementation. */
00541 void
00542 ec_GFp_extra_free_fp(ECGroup *group)
00543 {
00544        if (group->extra1 != NULL) {
00545               free(group->extra1);
00546               group->extra1 = NULL;
00547        }
00548 }
00549 
00550 /* Tests what precision floating point arithmetic is set to. This should
00551  * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa
00552  * (extended precision on x86) and sets it into the EC_group_fp. Returns
00553  * either 53 or 64 accordingly. */
00554 int
00555 ec_set_fp_precision(EC_group_fp * group)
00556 {
00557        double a = 9007199254740992.0;     /* 2^53 */
00558        double b = a + 1;
00559 
00560        if (a == b) {
00561               group->fpPrecision = 53;
00562               group->alpha = ecfp_alpha_53;
00563               return 53;
00564        }
00565        group->fpPrecision = 64;
00566        group->alpha = ecfp_alpha_64;
00567        return 64;
00568 }