Back to index

lightning-sunbird  0.9+nobinonly
ecp_fp.h
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 #ifndef __ecp_fp_h_
00040 #define __ecp_fp_h_
00041 
00042 #include "mpi.h"
00043 #include "ecl.h"
00044 #include "ecp.h"
00045 
00046 #include <sys/types.h>
00047 #include "mpi-priv.h"
00048 
00049 #ifdef ECL_DEBUG
00050 #include <assert.h>
00051 #endif
00052 
00053 /* Largest number of doubles to store one reduced number in floating
00054  * point. Used for memory allocation on the stack. */
00055 #define ECFP_MAXDOUBLES 10
00056 
00057 /* For debugging purposes */
00058 #ifndef ECL_DEBUG
00059 #define ECFP_ASSERT(x)
00060 #else
00061 #define ECFP_ASSERT(x) assert(x)
00062 #endif
00063 
00064 /* ECFP_Ti = 2^(i*24) Define as preprocessor constants so we can use in
00065  * multiple static constants */
00066 #define ECFP_T0 1.0
00067 #define ECFP_T1 16777216.0
00068 #define ECFP_T2 281474976710656.0
00069 #define ECFP_T3 4722366482869645213696.0
00070 #define ECFP_T4 79228162514264337593543950336.0
00071 #define ECFP_T5 1329227995784915872903807060280344576.0
00072 #define ECFP_T6 22300745198530623141535718272648361505980416.0
00073 #define ECFP_T7 374144419156711147060143317175368453031918731001856.0
00074 #define ECFP_T8 6277101735386680763835789423207666416102355444464034512896.0
00075 #define ECFP_T9 105312291668557186697918027683670432318895095400549111254310977536.0
00076 #define ECFP_T10  1766847064778384329583297500742918515827483896875618958121606201292619776.0
00077 #define ECFP_T11 29642774844752946028434172162224104410437116074403984394101141506025761187823616.0
00078 #define ECFP_T12 497323236409786642155382248146820840100456150797347717440463976893159497012533375533056.0
00079 #define ECFP_T13  8343699359066055009355553539724812947666814540455674882605631280555545803830627148527195652096.0
00080 #define ECFP_T14 139984046386112763159840142535527767382602843577165595931249318810236991948760059086304843329475444736.0
00081 #define ECFP_T15 2348542582773833227889480596789337027375682548908319870707290971532209025114608443463698998384768703031934976.0
00082 #define ECFP_T16 39402006196394479212279040100143613805079739270465446667948293404245\
00083 721771497210611414266254884915640806627990306816.0
00084 #define ECFP_T17 66105596879024859895191530803277103982840468296428121928464879527440\
00085 5791236311345825189210439715284847591212025023358304256.0
00086 #define ECFP_T18 11090678776483259438313656736572334813745748301503266300681918322458\
00087 485231222502492159897624416558312389564843845614287315896631296.0
00088 #define ECFP_T19 18607071341967536398062689481932916079453218833595342343206149099024\
00089 36577570298683715049089827234727835552055312041415509848580169253519\
00090 36.0
00091 
00092 #define ECFP_TWO160 1461501637330902918203684832716283019655932542976.0
00093 #define ECFP_TWO192 6277101735386680763835789423207666416102355444464034512896.0
00094 #define ECFP_TWO224 26959946667150639794667015087019630673637144422540572481103610249216.0
00095 
00096 /* Multiplicative constants */
00097 static const double ecfp_two32 = 4294967296.0;
00098 static const double ecfp_two64 = 18446744073709551616.0;
00099 static const double ecfp_twom16 = .0000152587890625;
00100 static const double ecfp_twom128 =
00101        .00000000000000000000000000000000000000293873587705571876992184134305561419454666389193021880377187926569604314863681793212890625;
00102 static const double ecfp_twom129 =
00103        .000000000000000000000000000000000000001469367938527859384960920671527807097273331945965109401885939632848021574318408966064453125;
00104 static const double ecfp_twom160 =
00105        .0000000000000000000000000000000000000000000000006842277657836020854119773355907793609766904013068924666782559979930620520927053718196475529111921787261962890625;
00106 static const double ecfp_twom192 =
00107        .000000000000000000000000000000000000000000000000000000000159309191113245227702888039776771180559110455519261878607388585338616290151305816094308987472018268594098344692611135542392730712890625;
00108 static const double ecfp_twom224 =
00109        .00000000000000000000000000000000000000000000000000000000000000000003709206150687421385731735261547639513367564778757791002453039058917581340095629358997312082723208437536338919136001159027049567384892725385725498199462890625;
00110 
00111 /* ecfp_exp[i] = 2^(i*ECFP_DSIZE) */
00112 static const double ecfp_exp[2 * ECFP_MAXDOUBLES] = {
00113        ECFP_T0, ECFP_T1, ECFP_T2, ECFP_T3, ECFP_T4, ECFP_T5,
00114        ECFP_T6, ECFP_T7, ECFP_T8, ECFP_T9, ECFP_T10, ECFP_T11,
00115        ECFP_T12, ECFP_T13, ECFP_T14, ECFP_T15, ECFP_T16, ECFP_T17, ECFP_T18,
00116        ECFP_T19
00117 };
00118 
00119 /* 1.1 * 2^52 Uses 2^52 to truncate, the .1 is an extra 2^51 to protect
00120  * the 2^52 bit, so that adding alphas to a negative number won't borrow
00121  * and empty the important 2^52 bit */
00122 #define ECFP_ALPHABASE_53 6755399441055744.0
00123 /* Special case: On some platforms, notably x86 Linux, there is an
00124  * extended-precision floating point representation with 64-bits of
00125  * precision in the mantissa.  These extra bits of precision require a
00126  * larger value of alpha to truncate, i.e. 1.1 * 2^63. */
00127 #define ECFP_ALPHABASE_64 13835058055282163712.0
00128 
00129 /* 
00130  * ecfp_alpha[i] = 1.5 * 2^(52 + i*ECFP_DSIZE) we add and subtract alpha
00131  * to truncate floating point numbers to a certain number of bits for
00132  * tidying */
00133 static const double ecfp_alpha_53[2 * ECFP_MAXDOUBLES] = {
00134        ECFP_ALPHABASE_53 * ECFP_T0,
00135        ECFP_ALPHABASE_53 * ECFP_T1,
00136        ECFP_ALPHABASE_53 * ECFP_T2,
00137        ECFP_ALPHABASE_53 * ECFP_T3,
00138        ECFP_ALPHABASE_53 * ECFP_T4,
00139        ECFP_ALPHABASE_53 * ECFP_T5,
00140        ECFP_ALPHABASE_53 * ECFP_T6,
00141        ECFP_ALPHABASE_53 * ECFP_T7,
00142        ECFP_ALPHABASE_53 * ECFP_T8,
00143        ECFP_ALPHABASE_53 * ECFP_T9,
00144        ECFP_ALPHABASE_53 * ECFP_T10,
00145        ECFP_ALPHABASE_53 * ECFP_T11,
00146        ECFP_ALPHABASE_53 * ECFP_T12,
00147        ECFP_ALPHABASE_53 * ECFP_T13,
00148        ECFP_ALPHABASE_53 * ECFP_T14,
00149        ECFP_ALPHABASE_53 * ECFP_T15,
00150        ECFP_ALPHABASE_53 * ECFP_T16,
00151        ECFP_ALPHABASE_53 * ECFP_T17,
00152        ECFP_ALPHABASE_53 * ECFP_T18,
00153        ECFP_ALPHABASE_53 * ECFP_T19
00154 };
00155 
00156 /* 
00157  * ecfp_alpha[i] = 1.5 * 2^(63 + i*ECFP_DSIZE) we add and subtract alpha
00158  * to truncate floating point numbers to a certain number of bits for
00159  * tidying */
00160 static const double ecfp_alpha_64[2 * ECFP_MAXDOUBLES] = {
00161        ECFP_ALPHABASE_64 * ECFP_T0,
00162        ECFP_ALPHABASE_64 * ECFP_T1,
00163        ECFP_ALPHABASE_64 * ECFP_T2,
00164        ECFP_ALPHABASE_64 * ECFP_T3,
00165        ECFP_ALPHABASE_64 * ECFP_T4,
00166        ECFP_ALPHABASE_64 * ECFP_T5,
00167        ECFP_ALPHABASE_64 * ECFP_T6,
00168        ECFP_ALPHABASE_64 * ECFP_T7,
00169        ECFP_ALPHABASE_64 * ECFP_T8,
00170        ECFP_ALPHABASE_64 * ECFP_T9,
00171        ECFP_ALPHABASE_64 * ECFP_T10,
00172        ECFP_ALPHABASE_64 * ECFP_T11,
00173        ECFP_ALPHABASE_64 * ECFP_T12,
00174        ECFP_ALPHABASE_64 * ECFP_T13,
00175        ECFP_ALPHABASE_64 * ECFP_T14,
00176        ECFP_ALPHABASE_64 * ECFP_T15,
00177        ECFP_ALPHABASE_64 * ECFP_T16,
00178        ECFP_ALPHABASE_64 * ECFP_T17,
00179        ECFP_ALPHABASE_64 * ECFP_T18,
00180        ECFP_ALPHABASE_64 * ECFP_T19
00181 };
00182 
00183 /* 0.011111111111111111111111 (binary) = 0.5 - 2^25 (24 ones) */
00184 #define ECFP_BETABASE 0.4999999701976776123046875
00185 
00186 /* 
00187  * We subtract beta prior to using alpha to simulate rounding down. We
00188  * make this close to 0.5 to round almost everything down, but exactly 0.5 
00189  * would cause some incorrect rounding. */
00190 static const double ecfp_beta[2 * ECFP_MAXDOUBLES] = {
00191        ECFP_BETABASE * ECFP_T0,
00192        ECFP_BETABASE * ECFP_T1,
00193        ECFP_BETABASE * ECFP_T2,
00194        ECFP_BETABASE * ECFP_T3,
00195        ECFP_BETABASE * ECFP_T4,
00196        ECFP_BETABASE * ECFP_T5,
00197        ECFP_BETABASE * ECFP_T6,
00198        ECFP_BETABASE * ECFP_T7,
00199        ECFP_BETABASE * ECFP_T8,
00200        ECFP_BETABASE * ECFP_T9,
00201        ECFP_BETABASE * ECFP_T10,
00202        ECFP_BETABASE * ECFP_T11,
00203        ECFP_BETABASE * ECFP_T12,
00204        ECFP_BETABASE * ECFP_T13,
00205        ECFP_BETABASE * ECFP_T14,
00206        ECFP_BETABASE * ECFP_T15,
00207        ECFP_BETABASE * ECFP_T16,
00208        ECFP_BETABASE * ECFP_T17,
00209        ECFP_BETABASE * ECFP_T18,
00210        ECFP_BETABASE * ECFP_T19
00211 };
00212 
00213 static const double ecfp_beta_160 = ECFP_BETABASE * ECFP_TWO160;
00214 static const double ecfp_beta_192 = ECFP_BETABASE * ECFP_TWO192;
00215 static const double ecfp_beta_224 = ECFP_BETABASE * ECFP_TWO224;
00216 
00217 /* Affine EC Point. This is the basic representation (x, y) of an elliptic 
00218  * curve point. */
00219 typedef struct {
00220        double x[ECFP_MAXDOUBLES];
00221        double y[ECFP_MAXDOUBLES];
00222 } ecfp_aff_pt;
00223 
00224 /* Jacobian EC Point. This coordinate system uses X = x/z^2, Y = y/z^3,
00225  * which enables calculations with fewer inversions than affine
00226  * coordinates. */
00227 typedef struct {
00228        double x[ECFP_MAXDOUBLES];
00229        double y[ECFP_MAXDOUBLES];
00230        double z[ECFP_MAXDOUBLES];
00231 } ecfp_jac_pt;
00232 
00233 /* Chudnovsky Jacobian EC Point. This coordinate system is the same as
00234  * Jacobian, except it keeps z^2, z^3 for faster additions. */
00235 typedef struct {
00236        double x[ECFP_MAXDOUBLES];
00237        double y[ECFP_MAXDOUBLES];
00238        double z[ECFP_MAXDOUBLES];
00239        double z2[ECFP_MAXDOUBLES];
00240        double z3[ECFP_MAXDOUBLES];
00241 } ecfp_chud_pt;
00242 
00243 /* Modified Jacobian EC Point. This coordinate system is the same as
00244  * Jacobian, except it keeps a*z^4 for faster doublings. */
00245 typedef struct {
00246        double x[ECFP_MAXDOUBLES];
00247        double y[ECFP_MAXDOUBLES];
00248        double z[ECFP_MAXDOUBLES];
00249        double az4[ECFP_MAXDOUBLES];
00250 } ecfp_jm_pt;
00251 
00252 struct EC_group_fp_str;
00253 
00254 typedef struct EC_group_fp_str EC_group_fp;
00255 struct EC_group_fp_str {
00256        int fpPrecision;                   /* Set to number of bits in mantissa, 53
00257                                                          * or 64 */
00258        int numDoubles;
00259        int primeBitSize;
00260        int orderBitSize;
00261        int doubleBitSize;
00262        int numInts;
00263        int aIsM3;                                /* True if curvea == -3 (mod p), then we
00264                                                          * can optimize doubling */
00265        double curvea[ECFP_MAXDOUBLES];
00266        /* Used to truncate a double to the number of bits in the curve */
00267        double bitSize_alpha;
00268        /* Pointer to either ecfp_alpha_53 or ecfp_alpha_64 */
00269        const double *alpha;
00270 
00271        void (*ecfp_singleReduce) (double *r, const EC_group_fp * group);
00272        void (*ecfp_reduce) (double *r, double *x, const EC_group_fp * group);
00273        /* Performs a "tidy" operation, which performs carrying, moving excess 
00274         * bits from one double to the next double, so that the precision of
00275         * the doubles is reduced to the regular precision ECFP_DSIZE. This
00276         * might result in some float digits being negative. */
00277        void (*ecfp_tidy) (double *t, const double *alpha,
00278                                       const EC_group_fp * group);
00279        /* Perform a point addition using coordinate system Jacobian + Affine
00280         * -> Jacobian. Input and output should be multi-precision floating
00281         * point integers. */
00282        void (*pt_add_jac_aff) (const ecfp_jac_pt * p, const ecfp_aff_pt * q,
00283                                                  ecfp_jac_pt * r, const EC_group_fp * group);
00284        /* Perform a point doubling in Jacobian coordinates. Input and output
00285         * should be multi-precision floating point integers. */
00286        void (*pt_dbl_jac) (const ecfp_jac_pt * dp, ecfp_jac_pt * dr,
00287                                           const EC_group_fp * group);
00288        /* Perform a point addition using Jacobian coordinate system. Input
00289         * and output should be multi-precision floating point integers. */
00290        void (*pt_add_jac) (const ecfp_jac_pt * p, const ecfp_jac_pt * q,
00291                                           ecfp_jac_pt * r, const EC_group_fp * group);
00292        /* Perform a point doubling in Modified Jacobian coordinates. Input
00293         * and output should be multi-precision floating point integers. */
00294        void (*pt_dbl_jm) (const ecfp_jm_pt * p, ecfp_jm_pt * r,
00295                                       const EC_group_fp * group);
00296        /* Perform a point doubling using coordinates Affine -> Chudnovsky
00297         * Jacobian. Input and output should be multi-precision floating point 
00298         * integers. */
00299        void (*pt_dbl_aff2chud) (const ecfp_aff_pt * p, ecfp_chud_pt * r,
00300                                                   const EC_group_fp * group);
00301        /* Perform a point addition using coordinates: Modified Jacobian +
00302         * Chudnovsky Jacobian -> Modified Jacobian. Input and output should
00303         * be multi-precision floating point integers. */
00304        void (*pt_add_jm_chud) (ecfp_jm_pt * p, ecfp_chud_pt * q,
00305                                                  ecfp_jm_pt * r, const EC_group_fp * group);
00306        /* Perform a point addition using Chudnovsky Jacobian coordinates.
00307         * Input and output should be multi-precision floating point integers. 
00308         */
00309        void (*pt_add_chud) (const ecfp_chud_pt * p, const ecfp_chud_pt * q,
00310                                            ecfp_chud_pt * r, const EC_group_fp * group);
00311        /* Expects out to be an array of size 16 of Chudnovsky Jacobian
00312         * points. Fills in Chudnovsky Jacobian form (x, y, z, z^2, z^3), for
00313         * -15P, -13P, -11P, -9P, -7P, -5P, -3P, -P, P, 3P, 5P, 7P, 9P, 11P,
00314         * 13P, 15P */
00315        void (*precompute_chud) (ecfp_chud_pt * out, const ecfp_aff_pt * p,
00316                                                   const EC_group_fp * group);
00317        /* Expects out to be an array of size 16 of Jacobian points. Fills in
00318         * Chudnovsky Jacobian form (x, y, z), for O, P, 2P, ... 15P */
00319        void (*precompute_jac) (ecfp_jac_pt * out, const ecfp_aff_pt * p,
00320                                                  const EC_group_fp * group);
00321 
00322 };
00323 
00324 /* Computes r = x*y.
00325  * r must be different (point to different memory) than x and y.
00326  * Does not tidy or reduce. */
00327 void ecfp_multiply(double *r, const double *x, const double *y);
00328 
00329 /* Performs a "tidy" operation, which performs carrying, moving excess
00330  * bits from one double to the next double, so that the precision of the
00331  * doubles is reduced to the regular precision group->doubleBitSize. This
00332  * might result in some float digits being negative. */
00333 void ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group);
00334 
00335 /* Performs tidying on only the upper float digits of a multi-precision
00336  * floating point integer, i.e. the digits beyond the regular length which 
00337  * are removed in the reduction step. */
00338 void ecfp_tidyUpper(double *t, const EC_group_fp * group);
00339 
00340 /* Performs tidying on a short multi-precision floating point integer (the 
00341  * lower group->numDoubles floats). */
00342 void ecfp_tidyShort(double *t, const EC_group_fp * group);
00343 
00344 /* Performs a more mathematically precise "tidying" so that each term is
00345  * positive.  This is slower than the regular tidying, and is used for
00346  * conversion from floating point to integer. */
00347 void ecfp_positiveTidy(double *t, const EC_group_fp * group);
00348 
00349 /* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters
00350  * a, b and p are the elliptic curve coefficients and the prime that
00351  * determines the field GFp.  Elliptic curve points P and R can be
00352  * identical.  Uses mixed Jacobian-affine coordinates. Uses 4-bit window
00353  * method. */
00354 mp_err
00355  ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px,
00356                                                  const mp_int *py, mp_int *rx, mp_int *ry,
00357                                                  const ECGroup *ecgroup);
00358 
00359 /* Computes R = nP where R is (rx, ry) and P is the base point. The
00360  * parameters a, b and p are the elliptic curve coefficients and the prime 
00361  * that determines the field GFp.  Elliptic curve points P and R can be
00362  * identical.  Uses mixed Jacobian-affine coordinates (Jacobian
00363  * coordinates for doubles and affine coordinates for additions; based on
00364  * recommendation from Brown et al.). Uses window NAF method (algorithm
00365  * 11) for scalar-point multiplication from Brown, Hankerson, Lopez,
00366  * Menezes. Software Implementation of the NIST Elliptic Curves Over Prime
00367  * Fields. */
00368 mp_err ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px,
00369                                                         const mp_int *py, mp_int *rx, mp_int *ry,
00370                                                         const ECGroup *ecgroup);
00371 
00372 /* Uses mixed Jacobian-affine coordinates to perform a point
00373  * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine
00374  * coordinates (Jacobian coordinates for doubles and affine coordinates
00375  * for additions; based on recommendation from Brown et al.). Not very
00376  * time efficient but quite space efficient, no precomputation needed.
00377  * group contains the elliptic curve coefficients and the prime that
00378  * determines the field GFp.  Elliptic curve points P and R can be
00379  * identical. Performs calculations in floating point number format, since 
00380  * this is faster than the integer operations on the ULTRASPARC III.
00381  * Uses left-to-right binary method (double & add) (algorithm 9) for
00382  * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes.
00383  * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */
00384 mp_err
00385  ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py,
00386                                      mp_int *rx, mp_int *ry, const ECGroup *ecgroup);
00387 
00388 /* Cleans up extra memory allocated in ECGroup for this implementation. */
00389 void ec_GFp_extra_free_fp(ECGroup *group);
00390 
00391 /* Converts from a floating point representation into an mp_int. Expects
00392  * that d is already reduced. */
00393 void
00394  ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup);
00395 
00396 /* Converts from an mpint into a floating point representation. */
00397 void
00398  ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup);
00399 
00400 /* Tests what precision floating point arithmetic is set to. This should
00401  * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa
00402  * (extended precision on x86) and sets it into the EC_group_fp. Returns
00403  * either 53 or 64 accordingly. */
00404 int ec_set_fp_precision(EC_group_fp * group);
00405 
00406 #endif