|
1 /* This Source Code Form is subject to the terms of the Mozilla Public |
|
2 * License, v. 2.0. If a copy of the MPL was not distributed with this |
|
3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ |
|
4 |
|
5 #ifndef __ecp_fp_h_ |
|
6 #define __ecp_fp_h_ |
|
7 |
|
8 #include "mpi.h" |
|
9 #include "ecl.h" |
|
10 #include "ecp.h" |
|
11 |
|
12 #include <sys/types.h> |
|
13 #include "mpi-priv.h" |
|
14 |
|
15 #ifdef ECL_DEBUG |
|
16 #include <assert.h> |
|
17 #endif |
|
18 |
|
19 /* Largest number of doubles to store one reduced number in floating |
|
20 * point. Used for memory allocation on the stack. */ |
|
21 #define ECFP_MAXDOUBLES 10 |
|
22 |
|
23 /* For debugging purposes */ |
|
24 #ifndef ECL_DEBUG |
|
25 #define ECFP_ASSERT(x) |
|
26 #else |
|
27 #define ECFP_ASSERT(x) assert(x) |
|
28 #endif |
|
29 |
|
30 /* ECFP_Ti = 2^(i*24) Define as preprocessor constants so we can use in |
|
31 * multiple static constants */ |
|
32 #define ECFP_T0 1.0 |
|
33 #define ECFP_T1 16777216.0 |
|
34 #define ECFP_T2 281474976710656.0 |
|
35 #define ECFP_T3 4722366482869645213696.0 |
|
36 #define ECFP_T4 79228162514264337593543950336.0 |
|
37 #define ECFP_T5 1329227995784915872903807060280344576.0 |
|
38 #define ECFP_T6 22300745198530623141535718272648361505980416.0 |
|
39 #define ECFP_T7 374144419156711147060143317175368453031918731001856.0 |
|
40 #define ECFP_T8 6277101735386680763835789423207666416102355444464034512896.0 |
|
41 #define ECFP_T9 105312291668557186697918027683670432318895095400549111254310977536.0 |
|
42 #define ECFP_T10 1766847064778384329583297500742918515827483896875618958121606201292619776.0 |
|
43 #define ECFP_T11 29642774844752946028434172162224104410437116074403984394101141506025761187823616.0 |
|
44 #define ECFP_T12 497323236409786642155382248146820840100456150797347717440463976893159497012533375533056.0 |
|
45 #define ECFP_T13 8343699359066055009355553539724812947666814540455674882605631280555545803830627148527195652096.0 |
|
46 #define ECFP_T14 139984046386112763159840142535527767382602843577165595931249318810236991948760059086304843329475444736.0 |
|
47 #define ECFP_T15 2348542582773833227889480596789337027375682548908319870707290971532209025114608443463698998384768703031934976.0 |
|
48 #define ECFP_T16 39402006196394479212279040100143613805079739270465446667948293404245\ |
|
49 721771497210611414266254884915640806627990306816.0 |
|
50 #define ECFP_T17 66105596879024859895191530803277103982840468296428121928464879527440\ |
|
51 5791236311345825189210439715284847591212025023358304256.0 |
|
52 #define ECFP_T18 11090678776483259438313656736572334813745748301503266300681918322458\ |
|
53 485231222502492159897624416558312389564843845614287315896631296.0 |
|
54 #define ECFP_T19 18607071341967536398062689481932916079453218833595342343206149099024\ |
|
55 36577570298683715049089827234727835552055312041415509848580169253519\ |
|
56 36.0 |
|
57 |
|
58 #define ECFP_TWO160 1461501637330902918203684832716283019655932542976.0 |
|
59 #define ECFP_TWO192 6277101735386680763835789423207666416102355444464034512896.0 |
|
60 #define ECFP_TWO224 26959946667150639794667015087019630673637144422540572481103610249216.0 |
|
61 |
|
62 /* Multiplicative constants */ |
|
63 static const double ecfp_two32 = 4294967296.0; |
|
64 static const double ecfp_two64 = 18446744073709551616.0; |
|
65 static const double ecfp_twom16 = .0000152587890625; |
|
66 static const double ecfp_twom128 = |
|
67 .00000000000000000000000000000000000000293873587705571876992184134305561419454666389193021880377187926569604314863681793212890625; |
|
68 static const double ecfp_twom129 = |
|
69 .000000000000000000000000000000000000001469367938527859384960920671527807097273331945965109401885939632848021574318408966064453125; |
|
70 static const double ecfp_twom160 = |
|
71 .0000000000000000000000000000000000000000000000006842277657836020854119773355907793609766904013068924666782559979930620520927053718196475529111921787261962890625; |
|
72 static const double ecfp_twom192 = |
|
73 .000000000000000000000000000000000000000000000000000000000159309191113245227702888039776771180559110455519261878607388585338616290151305816094308987472018268594098344692611135542392730712890625; |
|
74 static const double ecfp_twom224 = |
|
75 .00000000000000000000000000000000000000000000000000000000000000000003709206150687421385731735261547639513367564778757791002453039058917581340095629358997312082723208437536338919136001159027049567384892725385725498199462890625; |
|
76 |
|
77 /* ecfp_exp[i] = 2^(i*ECFP_DSIZE) */ |
|
78 static const double ecfp_exp[2 * ECFP_MAXDOUBLES] = { |
|
79 ECFP_T0, ECFP_T1, ECFP_T2, ECFP_T3, ECFP_T4, ECFP_T5, |
|
80 ECFP_T6, ECFP_T7, ECFP_T8, ECFP_T9, ECFP_T10, ECFP_T11, |
|
81 ECFP_T12, ECFP_T13, ECFP_T14, ECFP_T15, ECFP_T16, ECFP_T17, ECFP_T18, |
|
82 ECFP_T19 |
|
83 }; |
|
84 |
|
85 /* 1.1 * 2^52 Uses 2^52 to truncate, the .1 is an extra 2^51 to protect |
|
86 * the 2^52 bit, so that adding alphas to a negative number won't borrow |
|
87 * and empty the important 2^52 bit */ |
|
88 #define ECFP_ALPHABASE_53 6755399441055744.0 |
|
89 /* Special case: On some platforms, notably x86 Linux, there is an |
|
90 * extended-precision floating point representation with 64-bits of |
|
91 * precision in the mantissa. These extra bits of precision require a |
|
92 * larger value of alpha to truncate, i.e. 1.1 * 2^63. */ |
|
93 #define ECFP_ALPHABASE_64 13835058055282163712.0 |
|
94 |
|
95 /* |
|
96 * ecfp_alpha[i] = 1.5 * 2^(52 + i*ECFP_DSIZE) we add and subtract alpha |
|
97 * to truncate floating point numbers to a certain number of bits for |
|
98 * tidying */ |
|
99 static const double ecfp_alpha_53[2 * ECFP_MAXDOUBLES] = { |
|
100 ECFP_ALPHABASE_53 * ECFP_T0, |
|
101 ECFP_ALPHABASE_53 * ECFP_T1, |
|
102 ECFP_ALPHABASE_53 * ECFP_T2, |
|
103 ECFP_ALPHABASE_53 * ECFP_T3, |
|
104 ECFP_ALPHABASE_53 * ECFP_T4, |
|
105 ECFP_ALPHABASE_53 * ECFP_T5, |
|
106 ECFP_ALPHABASE_53 * ECFP_T6, |
|
107 ECFP_ALPHABASE_53 * ECFP_T7, |
|
108 ECFP_ALPHABASE_53 * ECFP_T8, |
|
109 ECFP_ALPHABASE_53 * ECFP_T9, |
|
110 ECFP_ALPHABASE_53 * ECFP_T10, |
|
111 ECFP_ALPHABASE_53 * ECFP_T11, |
|
112 ECFP_ALPHABASE_53 * ECFP_T12, |
|
113 ECFP_ALPHABASE_53 * ECFP_T13, |
|
114 ECFP_ALPHABASE_53 * ECFP_T14, |
|
115 ECFP_ALPHABASE_53 * ECFP_T15, |
|
116 ECFP_ALPHABASE_53 * ECFP_T16, |
|
117 ECFP_ALPHABASE_53 * ECFP_T17, |
|
118 ECFP_ALPHABASE_53 * ECFP_T18, |
|
119 ECFP_ALPHABASE_53 * ECFP_T19 |
|
120 }; |
|
121 |
|
122 /* |
|
123 * ecfp_alpha[i] = 1.5 * 2^(63 + i*ECFP_DSIZE) we add and subtract alpha |
|
124 * to truncate floating point numbers to a certain number of bits for |
|
125 * tidying */ |
|
126 static const double ecfp_alpha_64[2 * ECFP_MAXDOUBLES] = { |
|
127 ECFP_ALPHABASE_64 * ECFP_T0, |
|
128 ECFP_ALPHABASE_64 * ECFP_T1, |
|
129 ECFP_ALPHABASE_64 * ECFP_T2, |
|
130 ECFP_ALPHABASE_64 * ECFP_T3, |
|
131 ECFP_ALPHABASE_64 * ECFP_T4, |
|
132 ECFP_ALPHABASE_64 * ECFP_T5, |
|
133 ECFP_ALPHABASE_64 * ECFP_T6, |
|
134 ECFP_ALPHABASE_64 * ECFP_T7, |
|
135 ECFP_ALPHABASE_64 * ECFP_T8, |
|
136 ECFP_ALPHABASE_64 * ECFP_T9, |
|
137 ECFP_ALPHABASE_64 * ECFP_T10, |
|
138 ECFP_ALPHABASE_64 * ECFP_T11, |
|
139 ECFP_ALPHABASE_64 * ECFP_T12, |
|
140 ECFP_ALPHABASE_64 * ECFP_T13, |
|
141 ECFP_ALPHABASE_64 * ECFP_T14, |
|
142 ECFP_ALPHABASE_64 * ECFP_T15, |
|
143 ECFP_ALPHABASE_64 * ECFP_T16, |
|
144 ECFP_ALPHABASE_64 * ECFP_T17, |
|
145 ECFP_ALPHABASE_64 * ECFP_T18, |
|
146 ECFP_ALPHABASE_64 * ECFP_T19 |
|
147 }; |
|
148 |
|
149 /* 0.011111111111111111111111 (binary) = 0.5 - 2^25 (24 ones) */ |
|
150 #define ECFP_BETABASE 0.4999999701976776123046875 |
|
151 |
|
152 /* |
|
153 * We subtract beta prior to using alpha to simulate rounding down. We |
|
154 * make this close to 0.5 to round almost everything down, but exactly 0.5 |
|
155 * would cause some incorrect rounding. */ |
|
156 static const double ecfp_beta[2 * ECFP_MAXDOUBLES] = { |
|
157 ECFP_BETABASE * ECFP_T0, |
|
158 ECFP_BETABASE * ECFP_T1, |
|
159 ECFP_BETABASE * ECFP_T2, |
|
160 ECFP_BETABASE * ECFP_T3, |
|
161 ECFP_BETABASE * ECFP_T4, |
|
162 ECFP_BETABASE * ECFP_T5, |
|
163 ECFP_BETABASE * ECFP_T6, |
|
164 ECFP_BETABASE * ECFP_T7, |
|
165 ECFP_BETABASE * ECFP_T8, |
|
166 ECFP_BETABASE * ECFP_T9, |
|
167 ECFP_BETABASE * ECFP_T10, |
|
168 ECFP_BETABASE * ECFP_T11, |
|
169 ECFP_BETABASE * ECFP_T12, |
|
170 ECFP_BETABASE * ECFP_T13, |
|
171 ECFP_BETABASE * ECFP_T14, |
|
172 ECFP_BETABASE * ECFP_T15, |
|
173 ECFP_BETABASE * ECFP_T16, |
|
174 ECFP_BETABASE * ECFP_T17, |
|
175 ECFP_BETABASE * ECFP_T18, |
|
176 ECFP_BETABASE * ECFP_T19 |
|
177 }; |
|
178 |
|
179 static const double ecfp_beta_160 = ECFP_BETABASE * ECFP_TWO160; |
|
180 static const double ecfp_beta_192 = ECFP_BETABASE * ECFP_TWO192; |
|
181 static const double ecfp_beta_224 = ECFP_BETABASE * ECFP_TWO224; |
|
182 |
|
183 /* Affine EC Point. This is the basic representation (x, y) of an elliptic |
|
184 * curve point. */ |
|
185 typedef struct { |
|
186 double x[ECFP_MAXDOUBLES]; |
|
187 double y[ECFP_MAXDOUBLES]; |
|
188 } ecfp_aff_pt; |
|
189 |
|
190 /* Jacobian EC Point. This coordinate system uses X = x/z^2, Y = y/z^3, |
|
191 * which enables calculations with fewer inversions than affine |
|
192 * coordinates. */ |
|
193 typedef struct { |
|
194 double x[ECFP_MAXDOUBLES]; |
|
195 double y[ECFP_MAXDOUBLES]; |
|
196 double z[ECFP_MAXDOUBLES]; |
|
197 } ecfp_jac_pt; |
|
198 |
|
199 /* Chudnovsky Jacobian EC Point. This coordinate system is the same as |
|
200 * Jacobian, except it keeps z^2, z^3 for faster additions. */ |
|
201 typedef struct { |
|
202 double x[ECFP_MAXDOUBLES]; |
|
203 double y[ECFP_MAXDOUBLES]; |
|
204 double z[ECFP_MAXDOUBLES]; |
|
205 double z2[ECFP_MAXDOUBLES]; |
|
206 double z3[ECFP_MAXDOUBLES]; |
|
207 } ecfp_chud_pt; |
|
208 |
|
209 /* Modified Jacobian EC Point. This coordinate system is the same as |
|
210 * Jacobian, except it keeps a*z^4 for faster doublings. */ |
|
211 typedef struct { |
|
212 double x[ECFP_MAXDOUBLES]; |
|
213 double y[ECFP_MAXDOUBLES]; |
|
214 double z[ECFP_MAXDOUBLES]; |
|
215 double az4[ECFP_MAXDOUBLES]; |
|
216 } ecfp_jm_pt; |
|
217 |
|
218 struct EC_group_fp_str; |
|
219 |
|
220 typedef struct EC_group_fp_str EC_group_fp; |
|
221 struct EC_group_fp_str { |
|
222 int fpPrecision; /* Set to number of bits in mantissa, 53 |
|
223 * or 64 */ |
|
224 int numDoubles; |
|
225 int primeBitSize; |
|
226 int orderBitSize; |
|
227 int doubleBitSize; |
|
228 int numInts; |
|
229 int aIsM3; /* True if curvea == -3 (mod p), then we |
|
230 * can optimize doubling */ |
|
231 double curvea[ECFP_MAXDOUBLES]; |
|
232 /* Used to truncate a double to the number of bits in the curve */ |
|
233 double bitSize_alpha; |
|
234 /* Pointer to either ecfp_alpha_53 or ecfp_alpha_64 */ |
|
235 const double *alpha; |
|
236 |
|
237 void (*ecfp_singleReduce) (double *r, const EC_group_fp * group); |
|
238 void (*ecfp_reduce) (double *r, double *x, const EC_group_fp * group); |
|
239 /* Performs a "tidy" operation, which performs carrying, moving excess |
|
240 * bits from one double to the next double, so that the precision of |
|
241 * the doubles is reduced to the regular precision ECFP_DSIZE. This |
|
242 * might result in some float digits being negative. */ |
|
243 void (*ecfp_tidy) (double *t, const double *alpha, |
|
244 const EC_group_fp * group); |
|
245 /* Perform a point addition using coordinate system Jacobian + Affine |
|
246 * -> Jacobian. Input and output should be multi-precision floating |
|
247 * point integers. */ |
|
248 void (*pt_add_jac_aff) (const ecfp_jac_pt * p, const ecfp_aff_pt * q, |
|
249 ecfp_jac_pt * r, const EC_group_fp * group); |
|
250 /* Perform a point doubling in Jacobian coordinates. Input and output |
|
251 * should be multi-precision floating point integers. */ |
|
252 void (*pt_dbl_jac) (const ecfp_jac_pt * dp, ecfp_jac_pt * dr, |
|
253 const EC_group_fp * group); |
|
254 /* Perform a point addition using Jacobian coordinate system. Input |
|
255 * and output should be multi-precision floating point integers. */ |
|
256 void (*pt_add_jac) (const ecfp_jac_pt * p, const ecfp_jac_pt * q, |
|
257 ecfp_jac_pt * r, const EC_group_fp * group); |
|
258 /* Perform a point doubling in Modified Jacobian coordinates. Input |
|
259 * and output should be multi-precision floating point integers. */ |
|
260 void (*pt_dbl_jm) (const ecfp_jm_pt * p, ecfp_jm_pt * r, |
|
261 const EC_group_fp * group); |
|
262 /* Perform a point doubling using coordinates Affine -> Chudnovsky |
|
263 * Jacobian. Input and output should be multi-precision floating point |
|
264 * integers. */ |
|
265 void (*pt_dbl_aff2chud) (const ecfp_aff_pt * p, ecfp_chud_pt * r, |
|
266 const EC_group_fp * group); |
|
267 /* Perform a point addition using coordinates: Modified Jacobian + |
|
268 * Chudnovsky Jacobian -> Modified Jacobian. Input and output should |
|
269 * be multi-precision floating point integers. */ |
|
270 void (*pt_add_jm_chud) (ecfp_jm_pt * p, ecfp_chud_pt * q, |
|
271 ecfp_jm_pt * r, const EC_group_fp * group); |
|
272 /* Perform a point addition using Chudnovsky Jacobian coordinates. |
|
273 * Input and output should be multi-precision floating point integers. |
|
274 */ |
|
275 void (*pt_add_chud) (const ecfp_chud_pt * p, const ecfp_chud_pt * q, |
|
276 ecfp_chud_pt * r, const EC_group_fp * group); |
|
277 /* Expects out to be an array of size 16 of Chudnovsky Jacobian |
|
278 * points. Fills in Chudnovsky Jacobian form (x, y, z, z^2, z^3), for |
|
279 * -15P, -13P, -11P, -9P, -7P, -5P, -3P, -P, P, 3P, 5P, 7P, 9P, 11P, |
|
280 * 13P, 15P */ |
|
281 void (*precompute_chud) (ecfp_chud_pt * out, const ecfp_aff_pt * p, |
|
282 const EC_group_fp * group); |
|
283 /* Expects out to be an array of size 16 of Jacobian points. Fills in |
|
284 * Chudnovsky Jacobian form (x, y, z), for O, P, 2P, ... 15P */ |
|
285 void (*precompute_jac) (ecfp_jac_pt * out, const ecfp_aff_pt * p, |
|
286 const EC_group_fp * group); |
|
287 |
|
288 }; |
|
289 |
|
290 /* Computes r = x*y. |
|
291 * r must be different (point to different memory) than x and y. |
|
292 * Does not tidy or reduce. */ |
|
293 void ecfp_multiply(double *r, const double *x, const double *y); |
|
294 |
|
295 /* Performs a "tidy" operation, which performs carrying, moving excess |
|
296 * bits from one double to the next double, so that the precision of the |
|
297 * doubles is reduced to the regular precision group->doubleBitSize. This |
|
298 * might result in some float digits being negative. */ |
|
299 void ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group); |
|
300 |
|
301 /* Performs tidying on only the upper float digits of a multi-precision |
|
302 * floating point integer, i.e. the digits beyond the regular length which |
|
303 * are removed in the reduction step. */ |
|
304 void ecfp_tidyUpper(double *t, const EC_group_fp * group); |
|
305 |
|
306 /* Performs tidying on a short multi-precision floating point integer (the |
|
307 * lower group->numDoubles floats). */ |
|
308 void ecfp_tidyShort(double *t, const EC_group_fp * group); |
|
309 |
|
310 /* Performs a more mathematically precise "tidying" so that each term is |
|
311 * positive. This is slower than the regular tidying, and is used for |
|
312 * conversion from floating point to integer. */ |
|
313 void ecfp_positiveTidy(double *t, const EC_group_fp * group); |
|
314 |
|
315 /* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters |
|
316 * a, b and p are the elliptic curve coefficients and the prime that |
|
317 * determines the field GFp. Elliptic curve points P and R can be |
|
318 * identical. Uses mixed Jacobian-affine coordinates. Uses 4-bit window |
|
319 * method. */ |
|
320 mp_err |
|
321 ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px, |
|
322 const mp_int *py, mp_int *rx, mp_int *ry, |
|
323 const ECGroup *ecgroup); |
|
324 |
|
325 /* Computes R = nP where R is (rx, ry) and P is the base point. The |
|
326 * parameters a, b and p are the elliptic curve coefficients and the prime |
|
327 * that determines the field GFp. Elliptic curve points P and R can be |
|
328 * identical. Uses mixed Jacobian-affine coordinates (Jacobian |
|
329 * coordinates for doubles and affine coordinates for additions; based on |
|
330 * recommendation from Brown et al.). Uses window NAF method (algorithm |
|
331 * 11) for scalar-point multiplication from Brown, Hankerson, Lopez, |
|
332 * Menezes. Software Implementation of the NIST Elliptic Curves Over Prime |
|
333 * Fields. */ |
|
334 mp_err ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px, |
|
335 const mp_int *py, mp_int *rx, mp_int *ry, |
|
336 const ECGroup *ecgroup); |
|
337 |
|
338 /* Uses mixed Jacobian-affine coordinates to perform a point |
|
339 * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine |
|
340 * coordinates (Jacobian coordinates for doubles and affine coordinates |
|
341 * for additions; based on recommendation from Brown et al.). Not very |
|
342 * time efficient but quite space efficient, no precomputation needed. |
|
343 * group contains the elliptic curve coefficients and the prime that |
|
344 * determines the field GFp. Elliptic curve points P and R can be |
|
345 * identical. Performs calculations in floating point number format, since |
|
346 * this is faster than the integer operations on the ULTRASPARC III. |
|
347 * Uses left-to-right binary method (double & add) (algorithm 9) for |
|
348 * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes. |
|
349 * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */ |
|
350 mp_err |
|
351 ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py, |
|
352 mp_int *rx, mp_int *ry, const ECGroup *ecgroup); |
|
353 |
|
354 /* Cleans up extra memory allocated in ECGroup for this implementation. */ |
|
355 void ec_GFp_extra_free_fp(ECGroup *group); |
|
356 |
|
357 /* Converts from a floating point representation into an mp_int. Expects |
|
358 * that d is already reduced. */ |
|
359 void |
|
360 ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup); |
|
361 |
|
362 /* Converts from an mpint into a floating point representation. */ |
|
363 void |
|
364 ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup); |
|
365 |
|
366 /* Tests what precision floating point arithmetic is set to. This should |
|
367 * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa |
|
368 * (extended precision on x86) and sets it into the EC_group_fp. Returns |
|
369 * either 53 or 64 accordingly. */ |
|
370 int ec_set_fp_precision(EC_group_fp * group); |
|
371 |
|
372 #endif |