1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/security/nss/lib/freebl/ecl/ecp_fp.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,531 @@
1.4 +/* This Source Code Form is subject to the terms of the Mozilla Public
1.5 + * License, v. 2.0. If a copy of the MPL was not distributed with this
1.6 + * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
1.7 +
1.8 +#include "ecp_fp.h"
1.9 +#include "ecl-priv.h"
1.10 +#include <stdlib.h>
1.11 +
1.12 +/* Performs tidying on a short multi-precision floating point integer (the
1.13 + * lower group->numDoubles floats). */
1.14 +void
1.15 +ecfp_tidyShort(double *t, const EC_group_fp * group)
1.16 +{
1.17 + group->ecfp_tidy(t, group->alpha, group);
1.18 +}
1.19 +
1.20 +/* Performs tidying on only the upper float digits of a multi-precision
1.21 + * floating point integer, i.e. the digits beyond the regular length which
1.22 + * are removed in the reduction step. */
1.23 +void
1.24 +ecfp_tidyUpper(double *t, const EC_group_fp * group)
1.25 +{
1.26 + group->ecfp_tidy(t + group->numDoubles,
1.27 + group->alpha + group->numDoubles, group);
1.28 +}
1.29 +
1.30 +/* Performs a "tidy" operation, which performs carrying, moving excess
1.31 + * bits from one double to the next double, so that the precision of the
1.32 + * doubles is reduced to the regular precision group->doubleBitSize. This
1.33 + * might result in some float digits being negative. Alternative C version
1.34 + * for portability. */
1.35 +void
1.36 +ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group)
1.37 +{
1.38 + double q;
1.39 + int i;
1.40 +
1.41 + /* Do carrying */
1.42 + for (i = 0; i < group->numDoubles - 1; i++) {
1.43 + q = t[i] + alpha[i + 1];
1.44 + q -= alpha[i + 1];
1.45 + t[i] -= q;
1.46 + t[i + 1] += q;
1.47 +
1.48 + /* If we don't assume that truncation rounding is used, then q
1.49 + * might be 2^n bigger than expected (if it rounds up), then t[0]
1.50 + * could be negative and t[1] 2^n larger than expected. */
1.51 + }
1.52 +}
1.53 +
1.54 +/* Performs a more mathematically precise "tidying" so that each term is
1.55 + * positive. This is slower than the regular tidying, and is used for
1.56 + * conversion from floating point to integer. */
1.57 +void
1.58 +ecfp_positiveTidy(double *t, const EC_group_fp * group)
1.59 +{
1.60 + double q;
1.61 + int i;
1.62 +
1.63 + /* Do carrying */
1.64 + for (i = 0; i < group->numDoubles - 1; i++) {
1.65 + /* Subtract beta to force rounding down */
1.66 + q = t[i] - ecfp_beta[i + 1];
1.67 + q += group->alpha[i + 1];
1.68 + q -= group->alpha[i + 1];
1.69 + t[i] -= q;
1.70 + t[i + 1] += q;
1.71 +
1.72 + /* Due to subtracting ecfp_beta, we should have each term a
1.73 + * non-negative int */
1.74 + ECFP_ASSERT(t[i] / ecfp_exp[i] ==
1.75 + (unsigned long long) (t[i] / ecfp_exp[i]));
1.76 + ECFP_ASSERT(t[i] >= 0);
1.77 + }
1.78 +}
1.79 +
1.80 +/* Converts from a floating point representation into an mp_int. Expects
1.81 + * that d is already reduced. */
1.82 +void
1.83 +ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup)
1.84 +{
1.85 + EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
1.86 + unsigned short i16[(group->primeBitSize + 15) / 16];
1.87 + double q = 1;
1.88 +
1.89 +#ifdef ECL_THIRTY_TWO_BIT
1.90 + /* TEST uint32_t z = 0; */
1.91 + unsigned int z = 0;
1.92 +#else
1.93 + uint64_t z = 0;
1.94 +#endif
1.95 + int zBits = 0;
1.96 + int copiedBits = 0;
1.97 + int i = 0;
1.98 + int j = 0;
1.99 +
1.100 + mp_digit *out;
1.101 +
1.102 + /* Result should always be >= 0, so set sign accordingly */
1.103 + MP_SIGN(mpout) = MP_ZPOS;
1.104 +
1.105 + /* Tidy up so we're just dealing with positive numbers */
1.106 + ecfp_positiveTidy(d, group);
1.107 +
1.108 + /* We might need to do this reduction step more than once if the
1.109 + * reduction adds smaller terms which carry-over to cause another
1.110 + * reduction. However, this should happen very rarely, if ever,
1.111 + * depending on the elliptic curve. */
1.112 + do {
1.113 + /* Init loop data */
1.114 + z = 0;
1.115 + zBits = 0;
1.116 + q = 1;
1.117 + i = 0;
1.118 + j = 0;
1.119 + copiedBits = 0;
1.120 +
1.121 + /* Might have to do a bit more reduction */
1.122 + group->ecfp_singleReduce(d, group);
1.123 +
1.124 + /* Grow the size of the mpint if it's too small */
1.125 + s_mp_grow(mpout, group->numInts);
1.126 + MP_USED(mpout) = group->numInts;
1.127 + out = MP_DIGITS(mpout);
1.128 +
1.129 + /* Convert double to 16 bit integers */
1.130 + while (copiedBits < group->primeBitSize) {
1.131 + if (zBits < 16) {
1.132 + z += d[i] * q;
1.133 + i++;
1.134 + ECFP_ASSERT(i < (group->primeBitSize + 15) / 16);
1.135 + zBits += group->doubleBitSize;
1.136 + }
1.137 + i16[j] = z;
1.138 + j++;
1.139 + z >>= 16;
1.140 + zBits -= 16;
1.141 + q *= ecfp_twom16;
1.142 + copiedBits += 16;
1.143 + }
1.144 + } while (z != 0);
1.145 +
1.146 + /* Convert 16 bit integers to mp_digit */
1.147 +#ifdef ECL_THIRTY_TWO_BIT
1.148 + for (i = 0; i < (group->primeBitSize + 15) / 16; i += 2) {
1.149 + *out = 0;
1.150 + if (i + 1 < (group->primeBitSize + 15) / 16) {
1.151 + *out = i16[i + 1];
1.152 + *out <<= 16;
1.153 + }
1.154 + *out++ += i16[i];
1.155 + }
1.156 +#else /* 64 bit */
1.157 + for (i = 0; i < (group->primeBitSize + 15) / 16; i += 4) {
1.158 + *out = 0;
1.159 + if (i + 3 < (group->primeBitSize + 15) / 16) {
1.160 + *out = i16[i + 3];
1.161 + *out <<= 16;
1.162 + }
1.163 + if (i + 2 < (group->primeBitSize + 15) / 16) {
1.164 + *out += i16[i + 2];
1.165 + *out <<= 16;
1.166 + }
1.167 + if (i + 1 < (group->primeBitSize + 15) / 16) {
1.168 + *out += i16[i + 1];
1.169 + *out <<= 16;
1.170 + }
1.171 + *out++ += i16[i];
1.172 + }
1.173 +#endif
1.174 +
1.175 + /* Perform final reduction. mpout should already be the same number
1.176 + * of bits as p, but might not be less than p. Make it so. Since
1.177 + * mpout has the same number of bits as p, and 2p has a larger bit
1.178 + * size, then mpout < 2p, so a single subtraction of p will suffice. */
1.179 + if (mp_cmp(mpout, &ecgroup->meth->irr) >= 0) {
1.180 + mp_sub(mpout, &ecgroup->meth->irr, mpout);
1.181 + }
1.182 +
1.183 + /* Shrink the size of the mp_int to the actual used size (required for
1.184 + * mp_cmp_z == 0) */
1.185 + out = MP_DIGITS(mpout);
1.186 + for (i = group->numInts - 1; i > 0; i--) {
1.187 + if (out[i] != 0)
1.188 + break;
1.189 + }
1.190 + MP_USED(mpout) = i + 1;
1.191 +
1.192 + /* Should be between 0 and p-1 */
1.193 + ECFP_ASSERT(mp_cmp(mpout, &ecgroup->meth->irr) < 0);
1.194 + ECFP_ASSERT(mp_cmp_z(mpout) >= 0);
1.195 +}
1.196 +
1.197 +/* Converts from an mpint into a floating point representation. */
1.198 +void
1.199 +ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup)
1.200 +{
1.201 + int i;
1.202 + int j = 0;
1.203 + int size;
1.204 + double shift = 1;
1.205 + mp_digit *in;
1.206 + EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
1.207 +
1.208 +#ifdef ECL_DEBUG
1.209 + /* if debug mode, convert result back using ecfp_fp2i into cmp, then
1.210 + * compare to x. */
1.211 + mp_int cmp;
1.212 +
1.213 + MP_DIGITS(&cmp) = NULL;
1.214 + mp_init(&cmp);
1.215 +#endif
1.216 +
1.217 + ECFP_ASSERT(group != NULL);
1.218 +
1.219 + /* init output to 0 (since we skip over some terms) */
1.220 + for (i = 0; i < group->numDoubles; i++)
1.221 + out[i] = 0;
1.222 + i = 0;
1.223 +
1.224 + size = MP_USED(x);
1.225 + in = MP_DIGITS(x);
1.226 +
1.227 + /* Copy from int into doubles */
1.228 +#ifdef ECL_THIRTY_TWO_BIT
1.229 + while (j < size) {
1.230 + while (group->doubleBitSize * (i + 1) <= 32 * j) {
1.231 + i++;
1.232 + }
1.233 + ECFP_ASSERT(group->doubleBitSize * i <= 32 * j);
1.234 + out[i] = in[j];
1.235 + out[i] *= shift;
1.236 + shift *= ecfp_two32;
1.237 + j++;
1.238 + }
1.239 +#else
1.240 + while (j < size) {
1.241 + while (group->doubleBitSize * (i + 1) <= 64 * j) {
1.242 + i++;
1.243 + }
1.244 + ECFP_ASSERT(group->doubleBitSize * i <= 64 * j);
1.245 + out[i] = (in[j] & 0x00000000FFFFFFFF) * shift;
1.246 +
1.247 + while (group->doubleBitSize * (i + 1) <= 64 * j + 32) {
1.248 + i++;
1.249 + }
1.250 + ECFP_ASSERT(24 * i <= 64 * j + 32);
1.251 + out[i] = (in[j] & 0xFFFFFFFF00000000) * shift;
1.252 +
1.253 + shift *= ecfp_two64;
1.254 + j++;
1.255 + }
1.256 +#endif
1.257 + /* Realign bits to match double boundaries */
1.258 + ecfp_tidyShort(out, group);
1.259 +
1.260 +#ifdef ECL_DEBUG
1.261 + /* Convert result back to mp_int, compare to original */
1.262 + ecfp_fp2i(&cmp, out, ecgroup);
1.263 + ECFP_ASSERT(mp_cmp(&cmp, x) == 0);
1.264 + mp_clear(&cmp);
1.265 +#endif
1.266 +}
1.267 +
1.268 +/* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters
1.269 + * a, b and p are the elliptic curve coefficients and the prime that
1.270 + * determines the field GFp. Elliptic curve points P and R can be
1.271 + * identical. Uses Jacobian coordinates. Uses 4-bit window method. */
1.272 +mp_err
1.273 +ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px,
1.274 + const mp_int *py, mp_int *rx, mp_int *ry,
1.275 + const ECGroup *ecgroup)
1.276 +{
1.277 + mp_err res = MP_OKAY;
1.278 + ecfp_jac_pt precomp[16], r;
1.279 + ecfp_aff_pt p;
1.280 + EC_group_fp *group;
1.281 +
1.282 + mp_int rz;
1.283 + int i, ni, d;
1.284 +
1.285 + ARGCHK(ecgroup != NULL, MP_BADARG);
1.286 + ARGCHK((n != NULL) && (px != NULL) && (py != NULL), MP_BADARG);
1.287 +
1.288 + group = (EC_group_fp *) ecgroup->extra1;
1.289 + MP_DIGITS(&rz) = 0;
1.290 + MP_CHECKOK(mp_init(&rz));
1.291 +
1.292 + /* init p, da */
1.293 + ecfp_i2fp(p.x, px, ecgroup);
1.294 + ecfp_i2fp(p.y, py, ecgroup);
1.295 + ecfp_i2fp(group->curvea, &ecgroup->curvea, ecgroup);
1.296 +
1.297 + /* Do precomputation */
1.298 + group->precompute_jac(precomp, &p, group);
1.299 +
1.300 + /* Do main body of calculations */
1.301 + d = (mpl_significant_bits(n) + 3) / 4;
1.302 +
1.303 + /* R = inf */
1.304 + for (i = 0; i < group->numDoubles; i++) {
1.305 + r.z[i] = 0;
1.306 + }
1.307 +
1.308 + for (i = d - 1; i >= 0; i--) {
1.309 + /* compute window ni */
1.310 + ni = MP_GET_BIT(n, 4 * i + 3);
1.311 + ni <<= 1;
1.312 + ni |= MP_GET_BIT(n, 4 * i + 2);
1.313 + ni <<= 1;
1.314 + ni |= MP_GET_BIT(n, 4 * i + 1);
1.315 + ni <<= 1;
1.316 + ni |= MP_GET_BIT(n, 4 * i);
1.317 +
1.318 + /* R = 2^4 * R */
1.319 + group->pt_dbl_jac(&r, &r, group);
1.320 + group->pt_dbl_jac(&r, &r, group);
1.321 + group->pt_dbl_jac(&r, &r, group);
1.322 + group->pt_dbl_jac(&r, &r, group);
1.323 +
1.324 + /* R = R + (ni * P) */
1.325 + group->pt_add_jac(&r, &precomp[ni], &r, group);
1.326 + }
1.327 +
1.328 + /* Convert back to integer */
1.329 + ecfp_fp2i(rx, r.x, ecgroup);
1.330 + ecfp_fp2i(ry, r.y, ecgroup);
1.331 + ecfp_fp2i(&rz, r.z, ecgroup);
1.332 +
1.333 + /* convert result S to affine coordinates */
1.334 + MP_CHECKOK(ec_GFp_pt_jac2aff(rx, ry, &rz, rx, ry, ecgroup));
1.335 +
1.336 + CLEANUP:
1.337 + mp_clear(&rz);
1.338 + return res;
1.339 +}
1.340 +
1.341 +/* Uses mixed Jacobian-affine coordinates to perform a point
1.342 + * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine
1.343 + * coordinates (Jacobian coordinates for doubles and affine coordinates
1.344 + * for additions; based on recommendation from Brown et al.). Not very
1.345 + * time efficient but quite space efficient, no precomputation needed.
1.346 + * group contains the elliptic curve coefficients and the prime that
1.347 + * determines the field GFp. Elliptic curve points P and R can be
1.348 + * identical. Performs calculations in floating point number format, since
1.349 + * this is faster than the integer operations on the ULTRASPARC III.
1.350 + * Uses left-to-right binary method (double & add) (algorithm 9) for
1.351 + * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes.
1.352 + * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */
1.353 +mp_err
1.354 +ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py,
1.355 + mp_int *rx, mp_int *ry, const ECGroup *ecgroup)
1.356 +{
1.357 + mp_err res;
1.358 + mp_int sx, sy, sz;
1.359 +
1.360 + ecfp_aff_pt p;
1.361 + ecfp_jac_pt r;
1.362 + EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
1.363 +
1.364 + int i, l;
1.365 +
1.366 + MP_DIGITS(&sx) = 0;
1.367 + MP_DIGITS(&sy) = 0;
1.368 + MP_DIGITS(&sz) = 0;
1.369 + MP_CHECKOK(mp_init(&sx));
1.370 + MP_CHECKOK(mp_init(&sy));
1.371 + MP_CHECKOK(mp_init(&sz));
1.372 +
1.373 + /* if n = 0 then r = inf */
1.374 + if (mp_cmp_z(n) == 0) {
1.375 + mp_zero(rx);
1.376 + mp_zero(ry);
1.377 + res = MP_OKAY;
1.378 + goto CLEANUP;
1.379 + /* if n < 0 then out of range error */
1.380 + } else if (mp_cmp_z(n) < 0) {
1.381 + res = MP_RANGE;
1.382 + goto CLEANUP;
1.383 + }
1.384 +
1.385 + /* Convert from integer to floating point */
1.386 + ecfp_i2fp(p.x, px, ecgroup);
1.387 + ecfp_i2fp(p.y, py, ecgroup);
1.388 + ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
1.389 +
1.390 + /* Init r to point at infinity */
1.391 + for (i = 0; i < group->numDoubles; i++) {
1.392 + r.z[i] = 0;
1.393 + }
1.394 +
1.395 + /* double and add method */
1.396 + l = mpl_significant_bits(n) - 1;
1.397 +
1.398 + for (i = l; i >= 0; i--) {
1.399 + /* R = 2R */
1.400 + group->pt_dbl_jac(&r, &r, group);
1.401 +
1.402 + /* if n_i = 1, then R = R + Q */
1.403 + if (MP_GET_BIT(n, i) != 0) {
1.404 + group->pt_add_jac_aff(&r, &p, &r, group);
1.405 + }
1.406 + }
1.407 +
1.408 + /* Convert from floating point to integer */
1.409 + ecfp_fp2i(&sx, r.x, ecgroup);
1.410 + ecfp_fp2i(&sy, r.y, ecgroup);
1.411 + ecfp_fp2i(&sz, r.z, ecgroup);
1.412 +
1.413 + /* convert result R to affine coordinates */
1.414 + MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
1.415 +
1.416 + CLEANUP:
1.417 + mp_clear(&sx);
1.418 + mp_clear(&sy);
1.419 + mp_clear(&sz);
1.420 + return res;
1.421 +}
1.422 +
1.423 +/* Computes R = nP where R is (rx, ry) and P is the base point. Elliptic
1.424 + * curve points P and R can be identical. Uses mixed Modified-Jacobian
1.425 + * co-ordinates for doubling and Chudnovsky Jacobian coordinates for
1.426 + * additions. Uses 5-bit window NAF method (algorithm 11) for scalar-point
1.427 + * multiplication from Brown, Hankerson, Lopez, Menezes. Software
1.428 + * Implementation of the NIST Elliptic Curves Over Prime Fields. */
1.429 +mp_err
1.430 +ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px,
1.431 + const mp_int *py, mp_int *rx, mp_int *ry,
1.432 + const ECGroup *ecgroup)
1.433 +{
1.434 + mp_err res = MP_OKAY;
1.435 + mp_int sx, sy, sz;
1.436 + EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
1.437 + ecfp_chud_pt precomp[16];
1.438 +
1.439 + ecfp_aff_pt p;
1.440 + ecfp_jm_pt r;
1.441 +
1.442 + signed char naf[group->orderBitSize + 1];
1.443 + int i;
1.444 +
1.445 + MP_DIGITS(&sx) = 0;
1.446 + MP_DIGITS(&sy) = 0;
1.447 + MP_DIGITS(&sz) = 0;
1.448 + MP_CHECKOK(mp_init(&sx));
1.449 + MP_CHECKOK(mp_init(&sy));
1.450 + MP_CHECKOK(mp_init(&sz));
1.451 +
1.452 + /* if n = 0 then r = inf */
1.453 + if (mp_cmp_z(n) == 0) {
1.454 + mp_zero(rx);
1.455 + mp_zero(ry);
1.456 + res = MP_OKAY;
1.457 + goto CLEANUP;
1.458 + /* if n < 0 then out of range error */
1.459 + } else if (mp_cmp_z(n) < 0) {
1.460 + res = MP_RANGE;
1.461 + goto CLEANUP;
1.462 + }
1.463 +
1.464 + /* Convert from integer to floating point */
1.465 + ecfp_i2fp(p.x, px, ecgroup);
1.466 + ecfp_i2fp(p.y, py, ecgroup);
1.467 + ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
1.468 +
1.469 + /* Perform precomputation */
1.470 + group->precompute_chud(precomp, &p, group);
1.471 +
1.472 + /* Compute 5NAF */
1.473 + ec_compute_wNAF(naf, group->orderBitSize, n, 5);
1.474 +
1.475 + /* Init R = pt at infinity */
1.476 + for (i = 0; i < group->numDoubles; i++) {
1.477 + r.z[i] = 0;
1.478 + }
1.479 +
1.480 + /* wNAF method */
1.481 + for (i = group->orderBitSize; i >= 0; i--) {
1.482 + /* R = 2R */
1.483 + group->pt_dbl_jm(&r, &r, group);
1.484 +
1.485 + if (naf[i] != 0) {
1.486 + group->pt_add_jm_chud(&r, &precomp[(naf[i] + 15) / 2], &r,
1.487 + group);
1.488 + }
1.489 + }
1.490 +
1.491 + /* Convert from floating point to integer */
1.492 + ecfp_fp2i(&sx, r.x, ecgroup);
1.493 + ecfp_fp2i(&sy, r.y, ecgroup);
1.494 + ecfp_fp2i(&sz, r.z, ecgroup);
1.495 +
1.496 + /* convert result R to affine coordinates */
1.497 + MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
1.498 +
1.499 + CLEANUP:
1.500 + mp_clear(&sx);
1.501 + mp_clear(&sy);
1.502 + mp_clear(&sz);
1.503 + return res;
1.504 +}
1.505 +
1.506 +/* Cleans up extra memory allocated in ECGroup for this implementation. */
1.507 +void
1.508 +ec_GFp_extra_free_fp(ECGroup *group)
1.509 +{
1.510 + if (group->extra1 != NULL) {
1.511 + free(group->extra1);
1.512 + group->extra1 = NULL;
1.513 + }
1.514 +}
1.515 +
1.516 +/* Tests what precision floating point arithmetic is set to. This should
1.517 + * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa
1.518 + * (extended precision on x86) and sets it into the EC_group_fp. Returns
1.519 + * either 53 or 64 accordingly. */
1.520 +int
1.521 +ec_set_fp_precision(EC_group_fp * group)
1.522 +{
1.523 + double a = 9007199254740992.0; /* 2^53 */
1.524 + double b = a + 1;
1.525 +
1.526 + if (a == b) {
1.527 + group->fpPrecision = 53;
1.528 + group->alpha = ecfp_alpha_53;
1.529 + return 53;
1.530 + }
1.531 + group->fpPrecision = 64;
1.532 + group->alpha = ecfp_alpha_64;
1.533 + return 64;
1.534 +}
