/security/nss/lib/freebl/ecl/ecp_fp.c
C | 568 lines | 348 code | 76 blank | 144 comment | 52 complexity | a5686f450c1342809032f159bd28da30 MD5 | raw file
1/* 2 * ***** BEGIN LICENSE BLOCK ***** 3 * Version: MPL 1.1/GPL 2.0/LGPL 2.1 4 * 5 * The contents of this file are subject to the Mozilla Public License Version 6 * 1.1 (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * http://www.mozilla.org/MPL/ 9 * 10 * Software distributed under the License is distributed on an "AS IS" basis, 11 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 12 * for the specific language governing rights and limitations under the 13 * License. 14 * 15 * The Original Code is the elliptic curve math library for prime field curves 16 * using floating point operations. 17 * 18 * The Initial Developer of the Original Code is 19 * Sun Microsystems, Inc. 20 * Portions created by the Initial Developer are Copyright (C) 2003 21 * the Initial Developer. All Rights Reserved. 22 * 23 * Contributor(s): 24 * Sheueling Chang-Shantz <sheueling.chang@sun.com>, 25 * Stephen Fung <fungstep@hotmail.com>, and 26 * Douglas Stebila <douglas@stebila.ca>, Sun Microsystems Laboratories. 27 * 28 * Alternatively, the contents of this file may be used under the terms of 29 * either the GNU General Public License Version 2 or later (the "GPL"), or 30 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 31 * in which case the provisions of the GPL or the LGPL are applicable instead 32 * of those above. If you wish to allow use of your version of this file only 33 * under the terms of either the GPL or the LGPL, and not to allow others to 34 * use your version of this file under the terms of the MPL, indicate your 35 * decision by deleting the provisions above and replace them with the notice 36 * and other provisions required by the GPL or the LGPL. If you do not delete 37 * the provisions above, a recipient may use your version of this file under 38 * the terms of any one of the MPL, the GPL or the LGPL. 39 * 40 * ***** END LICENSE BLOCK ***** */ 41 42#include "ecp_fp.h" 43#include "ecl-priv.h" 44#include <stdlib.h> 45 46/* Performs tidying on a short multi-precision floating point integer (the 47 * lower group->numDoubles floats). */ 48void 49ecfp_tidyShort(double *t, const EC_group_fp * group) 50{ 51 group->ecfp_tidy(t, group->alpha, group); 52} 53 54/* Performs tidying on only the upper float digits of a multi-precision 55 * floating point integer, i.e. the digits beyond the regular length which 56 * are removed in the reduction step. */ 57void 58ecfp_tidyUpper(double *t, const EC_group_fp * group) 59{ 60 group->ecfp_tidy(t + group->numDoubles, 61 group->alpha + group->numDoubles, group); 62} 63 64/* Performs a "tidy" operation, which performs carrying, moving excess 65 * bits from one double to the next double, so that the precision of the 66 * doubles is reduced to the regular precision group->doubleBitSize. This 67 * might result in some float digits being negative. Alternative C version 68 * for portability. */ 69void 70ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group) 71{ 72 double q; 73 int i; 74 75 /* Do carrying */ 76 for (i = 0; i < group->numDoubles - 1; i++) { 77 q = t[i] + alpha[i + 1]; 78 q -= alpha[i + 1]; 79 t[i] -= q; 80 t[i + 1] += q; 81 82 /* If we don't assume that truncation rounding is used, then q 83 * might be 2^n bigger than expected (if it rounds up), then t[0] 84 * could be negative and t[1] 2^n larger than expected. */ 85 } 86} 87 88/* Performs a more mathematically precise "tidying" so that each term is 89 * positive. This is slower than the regular tidying, and is used for 90 * conversion from floating point to integer. */ 91void 92ecfp_positiveTidy(double *t, const EC_group_fp * group) 93{ 94 double q; 95 int i; 96 97 /* Do carrying */ 98 for (i = 0; i < group->numDoubles - 1; i++) { 99 /* Subtract beta to force rounding down */ 100 q = t[i] - ecfp_beta[i + 1]; 101 q += group->alpha[i + 1]; 102 q -= group->alpha[i + 1]; 103 t[i] -= q; 104 t[i + 1] += q; 105 106 /* Due to subtracting ecfp_beta, we should have each term a 107 * non-negative int */ 108 ECFP_ASSERT(t[i] / ecfp_exp[i] == 109 (unsigned long long) (t[i] / ecfp_exp[i])); 110 ECFP_ASSERT(t[i] >= 0); 111 } 112} 113 114/* Converts from a floating point representation into an mp_int. Expects 115 * that d is already reduced. */ 116void 117ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup) 118{ 119 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1; 120 unsigned short i16[(group->primeBitSize + 15) / 16]; 121 double q = 1; 122 123#ifdef ECL_THIRTY_TWO_BIT 124 /* TEST uint32_t z = 0; */ 125 unsigned int z = 0; 126#else 127 uint64_t z = 0; 128#endif 129 int zBits = 0; 130 int copiedBits = 0; 131 int i = 0; 132 int j = 0; 133 134 mp_digit *out; 135 136 /* Result should always be >= 0, so set sign accordingly */ 137 MP_SIGN(mpout) = MP_ZPOS; 138 139 /* Tidy up so we're just dealing with positive numbers */ 140 ecfp_positiveTidy(d, group); 141 142 /* We might need to do this reduction step more than once if the 143 * reduction adds smaller terms which carry-over to cause another 144 * reduction. However, this should happen very rarely, if ever, 145 * depending on the elliptic curve. */ 146 do { 147 /* Init loop data */ 148 z = 0; 149 zBits = 0; 150 q = 1; 151 i = 0; 152 j = 0; 153 copiedBits = 0; 154 155 /* Might have to do a bit more reduction */ 156 group->ecfp_singleReduce(d, group); 157 158 /* Grow the size of the mpint if it's too small */ 159 s_mp_grow(mpout, group->numInts); 160 MP_USED(mpout) = group->numInts; 161 out = MP_DIGITS(mpout); 162 163 /* Convert double to 16 bit integers */ 164 while (copiedBits < group->primeBitSize) { 165 if (zBits < 16) { 166 z += d[i] * q; 167 i++; 168 ECFP_ASSERT(i < (group->primeBitSize + 15) / 16); 169 zBits += group->doubleBitSize; 170 } 171 i16[j] = z; 172 j++; 173 z >>= 16; 174 zBits -= 16; 175 q *= ecfp_twom16; 176 copiedBits += 16; 177 } 178 } while (z != 0); 179 180 /* Convert 16 bit integers to mp_digit */ 181#ifdef ECL_THIRTY_TWO_BIT 182 for (i = 0; i < (group->primeBitSize + 15) / 16; i += 2) { 183 *out = 0; 184 if (i + 1 < (group->primeBitSize + 15) / 16) { 185 *out = i16[i + 1]; 186 *out <<= 16; 187 } 188 *out++ += i16[i]; 189 } 190#else /* 64 bit */ 191 for (i = 0; i < (group->primeBitSize + 15) / 16; i += 4) { 192 *out = 0; 193 if (i + 3 < (group->primeBitSize + 15) / 16) { 194 *out = i16[i + 3]; 195 *out <<= 16; 196 } 197 if (i + 2 < (group->primeBitSize + 15) / 16) { 198 *out += i16[i + 2]; 199 *out <<= 16; 200 } 201 if (i + 1 < (group->primeBitSize + 15) / 16) { 202 *out += i16[i + 1]; 203 *out <<= 16; 204 } 205 *out++ += i16[i]; 206 } 207#endif 208 209 /* Perform final reduction. mpout should already be the same number 210 * of bits as p, but might not be less than p. Make it so. Since 211 * mpout has the same number of bits as p, and 2p has a larger bit 212 * size, then mpout < 2p, so a single subtraction of p will suffice. */ 213 if (mp_cmp(mpout, &ecgroup->meth->irr) >= 0) { 214 mp_sub(mpout, &ecgroup->meth->irr, mpout); 215 } 216 217 /* Shrink the size of the mp_int to the actual used size (required for 218 * mp_cmp_z == 0) */ 219 out = MP_DIGITS(mpout); 220 for (i = group->numInts - 1; i > 0; i--) { 221 if (out[i] != 0) 222 break; 223 } 224 MP_USED(mpout) = i + 1; 225 226 /* Should be between 0 and p-1 */ 227 ECFP_ASSERT(mp_cmp(mpout, &ecgroup->meth->irr) < 0); 228 ECFP_ASSERT(mp_cmp_z(mpout) >= 0); 229} 230 231/* Converts from an mpint into a floating point representation. */ 232void 233ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup) 234{ 235 int i; 236 int j = 0; 237 int size; 238 double shift = 1; 239 mp_digit *in; 240 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1; 241 242#ifdef ECL_DEBUG 243 /* if debug mode, convert result back using ecfp_fp2i into cmp, then 244 * compare to x. */ 245 mp_int cmp; 246 247 MP_DIGITS(&cmp) = NULL; 248 mp_init(&cmp); 249#endif 250 251 ECFP_ASSERT(group != NULL); 252 253 /* init output to 0 (since we skip over some terms) */ 254 for (i = 0; i < group->numDoubles; i++) 255 out[i] = 0; 256 i = 0; 257 258 size = MP_USED(x); 259 in = MP_DIGITS(x); 260 261 /* Copy from int into doubles */ 262#ifdef ECL_THIRTY_TWO_BIT 263 while (j < size) { 264 while (group->doubleBitSize * (i + 1) <= 32 * j) { 265 i++; 266 } 267 ECFP_ASSERT(group->doubleBitSize * i <= 32 * j); 268 out[i] = in[j]; 269 out[i] *= shift; 270 shift *= ecfp_two32; 271 j++; 272 } 273#else 274 while (j < size) { 275 while (group->doubleBitSize * (i + 1) <= 64 * j) { 276 i++; 277 } 278 ECFP_ASSERT(group->doubleBitSize * i <= 64 * j); 279 out[i] = (in[j] & 0x00000000FFFFFFFF) * shift; 280 281 while (group->doubleBitSize * (i + 1) <= 64 * j + 32) { 282 i++; 283 } 284 ECFP_ASSERT(24 * i <= 64 * j + 32); 285 out[i] = (in[j] & 0xFFFFFFFF00000000) * shift; 286 287 shift *= ecfp_two64; 288 j++; 289 } 290#endif 291 /* Realign bits to match double boundaries */ 292 ecfp_tidyShort(out, group); 293 294#ifdef ECL_DEBUG 295 /* Convert result back to mp_int, compare to original */ 296 ecfp_fp2i(&cmp, out, ecgroup); 297 ECFP_ASSERT(mp_cmp(&cmp, x) == 0); 298 mp_clear(&cmp); 299#endif 300} 301 302/* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters 303 * a, b and p are the elliptic curve coefficients and the prime that 304 * determines the field GFp. Elliptic curve points P and R can be 305 * identical. Uses Jacobian coordinates. Uses 4-bit window method. */ 306mp_err 307ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px, 308 const mp_int *py, mp_int *rx, mp_int *ry, 309 const ECGroup *ecgroup) 310{ 311 mp_err res = MP_OKAY; 312 ecfp_jac_pt precomp[16], r; 313 ecfp_aff_pt p; 314 EC_group_fp *group; 315 316 mp_int rz; 317 int i, ni, d; 318 319 ARGCHK(ecgroup != NULL, MP_BADARG); 320 ARGCHK((n != NULL) && (px != NULL) && (py != NULL), MP_BADARG); 321 322 group = (EC_group_fp *) ecgroup->extra1; 323 MP_DIGITS(&rz) = 0; 324 MP_CHECKOK(mp_init(&rz)); 325 326 /* init p, da */ 327 ecfp_i2fp(p.x, px, ecgroup); 328 ecfp_i2fp(p.y, py, ecgroup); 329 ecfp_i2fp(group->curvea, &ecgroup->curvea, ecgroup); 330 331 /* Do precomputation */ 332 group->precompute_jac(precomp, &p, group); 333 334 /* Do main body of calculations */ 335 d = (mpl_significant_bits(n) + 3) / 4; 336 337 /* R = inf */ 338 for (i = 0; i < group->numDoubles; i++) { 339 r.z[i] = 0; 340 } 341 342 for (i = d - 1; i >= 0; i--) { 343 /* compute window ni */ 344 ni = MP_GET_BIT(n, 4 * i + 3); 345 ni <<= 1; 346 ni |= MP_GET_BIT(n, 4 * i + 2); 347 ni <<= 1; 348 ni |= MP_GET_BIT(n, 4 * i + 1); 349 ni <<= 1; 350 ni |= MP_GET_BIT(n, 4 * i); 351 352 /* R = 2^4 * R */ 353 group->pt_dbl_jac(&r, &r, group); 354 group->pt_dbl_jac(&r, &r, group); 355 group->pt_dbl_jac(&r, &r, group); 356 group->pt_dbl_jac(&r, &r, group); 357 358 /* R = R + (ni * P) */ 359 group->pt_add_jac(&r, &precomp[ni], &r, group); 360 } 361 362 /* Convert back to integer */ 363 ecfp_fp2i(rx, r.x, ecgroup); 364 ecfp_fp2i(ry, r.y, ecgroup); 365 ecfp_fp2i(&rz, r.z, ecgroup); 366 367 /* convert result S to affine coordinates */ 368 MP_CHECKOK(ec_GFp_pt_jac2aff(rx, ry, &rz, rx, ry, ecgroup)); 369 370 CLEANUP: 371 mp_clear(&rz); 372 return res; 373} 374 375/* Uses mixed Jacobian-affine coordinates to perform a point 376 * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine 377 * coordinates (Jacobian coordinates for doubles and affine coordinates 378 * for additions; based on recommendation from Brown et al.). Not very 379 * time efficient but quite space efficient, no precomputation needed. 380 * group contains the elliptic curve coefficients and the prime that 381 * determines the field GFp. Elliptic curve points P and R can be 382 * identical. Performs calculations in floating point number format, since 383 * this is faster than the integer operations on the ULTRASPARC III. 384 * Uses left-to-right binary method (double & add) (algorithm 9) for 385 * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes. 386 * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */ 387mp_err 388ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py, 389 mp_int *rx, mp_int *ry, const ECGroup *ecgroup) 390{ 391 mp_err res; 392 mp_int sx, sy, sz; 393 394 ecfp_aff_pt p; 395 ecfp_jac_pt r; 396 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1; 397 398 int i, l; 399 400 MP_DIGITS(&sx) = 0; 401 MP_DIGITS(&sy) = 0; 402 MP_DIGITS(&sz) = 0; 403 MP_CHECKOK(mp_init(&sx)); 404 MP_CHECKOK(mp_init(&sy)); 405 MP_CHECKOK(mp_init(&sz)); 406 407 /* if n = 0 then r = inf */ 408 if (mp_cmp_z(n) == 0) { 409 mp_zero(rx); 410 mp_zero(ry); 411 res = MP_OKAY; 412 goto CLEANUP; 413 /* if n < 0 then out of range error */ 414 } else if (mp_cmp_z(n) < 0) { 415 res = MP_RANGE; 416 goto CLEANUP; 417 } 418 419 /* Convert from integer to floating point */ 420 ecfp_i2fp(p.x, px, ecgroup); 421 ecfp_i2fp(p.y, py, ecgroup); 422 ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup); 423 424 /* Init r to point at infinity */ 425 for (i = 0; i < group->numDoubles; i++) { 426 r.z[i] = 0; 427 } 428 429 /* double and add method */ 430 l = mpl_significant_bits(n) - 1; 431 432 for (i = l; i >= 0; i--) { 433 /* R = 2R */ 434 group->pt_dbl_jac(&r, &r, group); 435 436 /* if n_i = 1, then R = R + Q */ 437 if (MP_GET_BIT(n, i) != 0) { 438 group->pt_add_jac_aff(&r, &p, &r, group); 439 } 440 } 441 442 /* Convert from floating point to integer */ 443 ecfp_fp2i(&sx, r.x, ecgroup); 444 ecfp_fp2i(&sy, r.y, ecgroup); 445 ecfp_fp2i(&sz, r.z, ecgroup); 446 447 /* convert result R to affine coordinates */ 448 MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup)); 449 450 CLEANUP: 451 mp_clear(&sx); 452 mp_clear(&sy); 453 mp_clear(&sz); 454 return res; 455} 456 457/* Computes R = nP where R is (rx, ry) and P is the base point. Elliptic 458 * curve points P and R can be identical. Uses mixed Modified-Jacobian 459 * co-ordinates for doubling and Chudnovsky Jacobian coordinates for 460 * additions. Uses 5-bit window NAF method (algorithm 11) for scalar-point 461 * multiplication from Brown, Hankerson, Lopez, Menezes. Software 462 * Implementation of the NIST Elliptic Curves Over Prime Fields. */ 463mp_err 464ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px, 465 const mp_int *py, mp_int *rx, mp_int *ry, 466 const ECGroup *ecgroup) 467{ 468 mp_err res = MP_OKAY; 469 mp_int sx, sy, sz; 470 EC_group_fp *group = (EC_group_fp *) ecgroup->extra1; 471 ecfp_chud_pt precomp[16]; 472 473 ecfp_aff_pt p; 474 ecfp_jm_pt r; 475 476 signed char naf[group->orderBitSize + 1]; 477 int i; 478 479 MP_DIGITS(&sx) = 0; 480 MP_DIGITS(&sy) = 0; 481 MP_DIGITS(&sz) = 0; 482 MP_CHECKOK(mp_init(&sx)); 483 MP_CHECKOK(mp_init(&sy)); 484 MP_CHECKOK(mp_init(&sz)); 485 486 /* if n = 0 then r = inf */ 487 if (mp_cmp_z(n) == 0) { 488 mp_zero(rx); 489 mp_zero(ry); 490 res = MP_OKAY; 491 goto CLEANUP; 492 /* if n < 0 then out of range error */ 493 } else if (mp_cmp_z(n) < 0) { 494 res = MP_RANGE; 495 goto CLEANUP; 496 } 497 498 /* Convert from integer to floating point */ 499 ecfp_i2fp(p.x, px, ecgroup); 500 ecfp_i2fp(p.y, py, ecgroup); 501 ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup); 502 503 /* Perform precomputation */ 504 group->precompute_chud(precomp, &p, group); 505 506 /* Compute 5NAF */ 507 ec_compute_wNAF(naf, group->orderBitSize, n, 5); 508 509 /* Init R = pt at infinity */ 510 for (i = 0; i < group->numDoubles; i++) { 511 r.z[i] = 0; 512 } 513 514 /* wNAF method */ 515 for (i = group->orderBitSize; i >= 0; i--) { 516 /* R = 2R */ 517 group->pt_dbl_jm(&r, &r, group); 518 519 if (naf[i] != 0) { 520 group->pt_add_jm_chud(&r, &precomp[(naf[i] + 15) / 2], &r, 521 group); 522 } 523 } 524 525 /* Convert from floating point to integer */ 526 ecfp_fp2i(&sx, r.x, ecgroup); 527 ecfp_fp2i(&sy, r.y, ecgroup); 528 ecfp_fp2i(&sz, r.z, ecgroup); 529 530 /* convert result R to affine coordinates */ 531 MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup)); 532 533 CLEANUP: 534 mp_clear(&sx); 535 mp_clear(&sy); 536 mp_clear(&sz); 537 return res; 538} 539 540/* Cleans up extra memory allocated in ECGroup for this implementation. */ 541void 542ec_GFp_extra_free_fp(ECGroup *group) 543{ 544 if (group->extra1 != NULL) { 545 free(group->extra1); 546 group->extra1 = NULL; 547 } 548} 549 550/* Tests what precision floating point arithmetic is set to. This should 551 * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa 552 * (extended precision on x86) and sets it into the EC_group_fp. Returns 553 * either 53 or 64 accordingly. */ 554int 555ec_set_fp_precision(EC_group_fp * group) 556{ 557 double a = 9007199254740992.0; /* 2^53 */ 558 double b = a + 1; 559 560 if (a == b) { 561 group->fpPrecision = 53; 562 group->alpha = ecfp_alpha_53; 563 return 53; 564 } 565 group->fpPrecision = 64; 566 group->alpha = ecfp_alpha_64; 567 return 64; 568}