/security/nss/lib/freebl/ecl/ecp_fp.c

http://github.com/zpao/v8monkey · 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. #include "ecp_fp.h"
  42. #include "ecl-priv.h"
  43. #include <stdlib.h>
  44. /* Performs tidying on a short multi-precision floating point integer (the
  45. * lower group->numDoubles floats). */
  46. void
  47. ecfp_tidyShort(double *t, const EC_group_fp * group)
  48. {
  49. group->ecfp_tidy(t, group->alpha, group);
  50. }
  51. /* Performs tidying on only the upper float digits of a multi-precision
  52. * floating point integer, i.e. the digits beyond the regular length which
  53. * are removed in the reduction step. */
  54. void
  55. ecfp_tidyUpper(double *t, const EC_group_fp * group)
  56. {
  57. group->ecfp_tidy(t + group->numDoubles,
  58. group->alpha + group->numDoubles, group);
  59. }
  60. /* Performs a "tidy" operation, which performs carrying, moving excess
  61. * bits from one double to the next double, so that the precision of the
  62. * doubles is reduced to the regular precision group->doubleBitSize. This
  63. * might result in some float digits being negative. Alternative C version
  64. * for portability. */
  65. void
  66. ecfp_tidy(double *t, const double *alpha, const EC_group_fp * group)
  67. {
  68. double q;
  69. int i;
  70. /* Do carrying */
  71. for (i = 0; i < group->numDoubles - 1; i++) {
  72. q = t[i] + alpha[i + 1];
  73. q -= alpha[i + 1];
  74. t[i] -= q;
  75. t[i + 1] += q;
  76. /* If we don't assume that truncation rounding is used, then q
  77. * might be 2^n bigger than expected (if it rounds up), then t[0]
  78. * could be negative and t[1] 2^n larger than expected. */
  79. }
  80. }
  81. /* Performs a more mathematically precise "tidying" so that each term is
  82. * positive. This is slower than the regular tidying, and is used for
  83. * conversion from floating point to integer. */
  84. void
  85. ecfp_positiveTidy(double *t, const EC_group_fp * group)
  86. {
  87. double q;
  88. int i;
  89. /* Do carrying */
  90. for (i = 0; i < group->numDoubles - 1; i++) {
  91. /* Subtract beta to force rounding down */
  92. q = t[i] - ecfp_beta[i + 1];
  93. q += group->alpha[i + 1];
  94. q -= group->alpha[i + 1];
  95. t[i] -= q;
  96. t[i + 1] += q;
  97. /* Due to subtracting ecfp_beta, we should have each term a
  98. * non-negative int */
  99. ECFP_ASSERT(t[i] / ecfp_exp[i] ==
  100. (unsigned long long) (t[i] / ecfp_exp[i]));
  101. ECFP_ASSERT(t[i] >= 0);
  102. }
  103. }
  104. /* Converts from a floating point representation into an mp_int. Expects
  105. * that d is already reduced. */
  106. void
  107. ecfp_fp2i(mp_int *mpout, double *d, const ECGroup *ecgroup)
  108. {
  109. EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
  110. unsigned short i16[(group->primeBitSize + 15) / 16];
  111. double q = 1;
  112. #ifdef ECL_THIRTY_TWO_BIT
  113. /* TEST uint32_t z = 0; */
  114. unsigned int z = 0;
  115. #else
  116. uint64_t z = 0;
  117. #endif
  118. int zBits = 0;
  119. int copiedBits = 0;
  120. int i = 0;
  121. int j = 0;
  122. mp_digit *out;
  123. /* Result should always be >= 0, so set sign accordingly */
  124. MP_SIGN(mpout) = MP_ZPOS;
  125. /* Tidy up so we're just dealing with positive numbers */
  126. ecfp_positiveTidy(d, group);
  127. /* We might need to do this reduction step more than once if the
  128. * reduction adds smaller terms which carry-over to cause another
  129. * reduction. However, this should happen very rarely, if ever,
  130. * depending on the elliptic curve. */
  131. do {
  132. /* Init loop data */
  133. z = 0;
  134. zBits = 0;
  135. q = 1;
  136. i = 0;
  137. j = 0;
  138. copiedBits = 0;
  139. /* Might have to do a bit more reduction */
  140. group->ecfp_singleReduce(d, group);
  141. /* Grow the size of the mpint if it's too small */
  142. s_mp_grow(mpout, group->numInts);
  143. MP_USED(mpout) = group->numInts;
  144. out = MP_DIGITS(mpout);
  145. /* Convert double to 16 bit integers */
  146. while (copiedBits < group->primeBitSize) {
  147. if (zBits < 16) {
  148. z += d[i] * q;
  149. i++;
  150. ECFP_ASSERT(i < (group->primeBitSize + 15) / 16);
  151. zBits += group->doubleBitSize;
  152. }
  153. i16[j] = z;
  154. j++;
  155. z >>= 16;
  156. zBits -= 16;
  157. q *= ecfp_twom16;
  158. copiedBits += 16;
  159. }
  160. } while (z != 0);
  161. /* Convert 16 bit integers to mp_digit */
  162. #ifdef ECL_THIRTY_TWO_BIT
  163. for (i = 0; i < (group->primeBitSize + 15) / 16; i += 2) {
  164. *out = 0;
  165. if (i + 1 < (group->primeBitSize + 15) / 16) {
  166. *out = i16[i + 1];
  167. *out <<= 16;
  168. }
  169. *out++ += i16[i];
  170. }
  171. #else /* 64 bit */
  172. for (i = 0; i < (group->primeBitSize + 15) / 16; i += 4) {
  173. *out = 0;
  174. if (i + 3 < (group->primeBitSize + 15) / 16) {
  175. *out = i16[i + 3];
  176. *out <<= 16;
  177. }
  178. if (i + 2 < (group->primeBitSize + 15) / 16) {
  179. *out += i16[i + 2];
  180. *out <<= 16;
  181. }
  182. if (i + 1 < (group->primeBitSize + 15) / 16) {
  183. *out += i16[i + 1];
  184. *out <<= 16;
  185. }
  186. *out++ += i16[i];
  187. }
  188. #endif
  189. /* Perform final reduction. mpout should already be the same number
  190. * of bits as p, but might not be less than p. Make it so. Since
  191. * mpout has the same number of bits as p, and 2p has a larger bit
  192. * size, then mpout < 2p, so a single subtraction of p will suffice. */
  193. if (mp_cmp(mpout, &ecgroup->meth->irr) >= 0) {
  194. mp_sub(mpout, &ecgroup->meth->irr, mpout);
  195. }
  196. /* Shrink the size of the mp_int to the actual used size (required for
  197. * mp_cmp_z == 0) */
  198. out = MP_DIGITS(mpout);
  199. for (i = group->numInts - 1; i > 0; i--) {
  200. if (out[i] != 0)
  201. break;
  202. }
  203. MP_USED(mpout) = i + 1;
  204. /* Should be between 0 and p-1 */
  205. ECFP_ASSERT(mp_cmp(mpout, &ecgroup->meth->irr) < 0);
  206. ECFP_ASSERT(mp_cmp_z(mpout) >= 0);
  207. }
  208. /* Converts from an mpint into a floating point representation. */
  209. void
  210. ecfp_i2fp(double *out, const mp_int *x, const ECGroup *ecgroup)
  211. {
  212. int i;
  213. int j = 0;
  214. int size;
  215. double shift = 1;
  216. mp_digit *in;
  217. EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
  218. #ifdef ECL_DEBUG
  219. /* if debug mode, convert result back using ecfp_fp2i into cmp, then
  220. * compare to x. */
  221. mp_int cmp;
  222. MP_DIGITS(&cmp) = NULL;
  223. mp_init(&cmp);
  224. #endif
  225. ECFP_ASSERT(group != NULL);
  226. /* init output to 0 (since we skip over some terms) */
  227. for (i = 0; i < group->numDoubles; i++)
  228. out[i] = 0;
  229. i = 0;
  230. size = MP_USED(x);
  231. in = MP_DIGITS(x);
  232. /* Copy from int into doubles */
  233. #ifdef ECL_THIRTY_TWO_BIT
  234. while (j < size) {
  235. while (group->doubleBitSize * (i + 1) <= 32 * j) {
  236. i++;
  237. }
  238. ECFP_ASSERT(group->doubleBitSize * i <= 32 * j);
  239. out[i] = in[j];
  240. out[i] *= shift;
  241. shift *= ecfp_two32;
  242. j++;
  243. }
  244. #else
  245. while (j < size) {
  246. while (group->doubleBitSize * (i + 1) <= 64 * j) {
  247. i++;
  248. }
  249. ECFP_ASSERT(group->doubleBitSize * i <= 64 * j);
  250. out[i] = (in[j] & 0x00000000FFFFFFFF) * shift;
  251. while (group->doubleBitSize * (i + 1) <= 64 * j + 32) {
  252. i++;
  253. }
  254. ECFP_ASSERT(24 * i <= 64 * j + 32);
  255. out[i] = (in[j] & 0xFFFFFFFF00000000) * shift;
  256. shift *= ecfp_two64;
  257. j++;
  258. }
  259. #endif
  260. /* Realign bits to match double boundaries */
  261. ecfp_tidyShort(out, group);
  262. #ifdef ECL_DEBUG
  263. /* Convert result back to mp_int, compare to original */
  264. ecfp_fp2i(&cmp, out, ecgroup);
  265. ECFP_ASSERT(mp_cmp(&cmp, x) == 0);
  266. mp_clear(&cmp);
  267. #endif
  268. }
  269. /* Computes R = nP where R is (rx, ry) and P is (px, py). The parameters
  270. * a, b and p are the elliptic curve coefficients and the prime that
  271. * determines the field GFp. Elliptic curve points P and R can be
  272. * identical. Uses Jacobian coordinates. Uses 4-bit window method. */
  273. mp_err
  274. ec_GFp_point_mul_jac_4w_fp(const mp_int *n, const mp_int *px,
  275. const mp_int *py, mp_int *rx, mp_int *ry,
  276. const ECGroup *ecgroup)
  277. {
  278. mp_err res = MP_OKAY;
  279. ecfp_jac_pt precomp[16], r;
  280. ecfp_aff_pt p;
  281. EC_group_fp *group;
  282. mp_int rz;
  283. int i, ni, d;
  284. ARGCHK(ecgroup != NULL, MP_BADARG);
  285. ARGCHK((n != NULL) && (px != NULL) && (py != NULL), MP_BADARG);
  286. group = (EC_group_fp *) ecgroup->extra1;
  287. MP_DIGITS(&rz) = 0;
  288. MP_CHECKOK(mp_init(&rz));
  289. /* init p, da */
  290. ecfp_i2fp(p.x, px, ecgroup);
  291. ecfp_i2fp(p.y, py, ecgroup);
  292. ecfp_i2fp(group->curvea, &ecgroup->curvea, ecgroup);
  293. /* Do precomputation */
  294. group->precompute_jac(precomp, &p, group);
  295. /* Do main body of calculations */
  296. d = (mpl_significant_bits(n) + 3) / 4;
  297. /* R = inf */
  298. for (i = 0; i < group->numDoubles; i++) {
  299. r.z[i] = 0;
  300. }
  301. for (i = d - 1; i >= 0; i--) {
  302. /* compute window ni */
  303. ni = MP_GET_BIT(n, 4 * i + 3);
  304. ni <<= 1;
  305. ni |= MP_GET_BIT(n, 4 * i + 2);
  306. ni <<= 1;
  307. ni |= MP_GET_BIT(n, 4 * i + 1);
  308. ni <<= 1;
  309. ni |= MP_GET_BIT(n, 4 * i);
  310. /* R = 2^4 * R */
  311. group->pt_dbl_jac(&r, &r, group);
  312. group->pt_dbl_jac(&r, &r, group);
  313. group->pt_dbl_jac(&r, &r, group);
  314. group->pt_dbl_jac(&r, &r, group);
  315. /* R = R + (ni * P) */
  316. group->pt_add_jac(&r, &precomp[ni], &r, group);
  317. }
  318. /* Convert back to integer */
  319. ecfp_fp2i(rx, r.x, ecgroup);
  320. ecfp_fp2i(ry, r.y, ecgroup);
  321. ecfp_fp2i(&rz, r.z, ecgroup);
  322. /* convert result S to affine coordinates */
  323. MP_CHECKOK(ec_GFp_pt_jac2aff(rx, ry, &rz, rx, ry, ecgroup));
  324. CLEANUP:
  325. mp_clear(&rz);
  326. return res;
  327. }
  328. /* Uses mixed Jacobian-affine coordinates to perform a point
  329. * multiplication: R = n * P, n scalar. Uses mixed Jacobian-affine
  330. * coordinates (Jacobian coordinates for doubles and affine coordinates
  331. * for additions; based on recommendation from Brown et al.). Not very
  332. * time efficient but quite space efficient, no precomputation needed.
  333. * group contains the elliptic curve coefficients and the prime that
  334. * determines the field GFp. Elliptic curve points P and R can be
  335. * identical. Performs calculations in floating point number format, since
  336. * this is faster than the integer operations on the ULTRASPARC III.
  337. * Uses left-to-right binary method (double & add) (algorithm 9) for
  338. * scalar-point multiplication from Brown, Hankerson, Lopez, Menezes.
  339. * Software Implementation of the NIST Elliptic Curves Over Prime Fields. */
  340. mp_err
  341. ec_GFp_pt_mul_jac_fp(const mp_int *n, const mp_int *px, const mp_int *py,
  342. mp_int *rx, mp_int *ry, const ECGroup *ecgroup)
  343. {
  344. mp_err res;
  345. mp_int sx, sy, sz;
  346. ecfp_aff_pt p;
  347. ecfp_jac_pt r;
  348. EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
  349. int i, l;
  350. MP_DIGITS(&sx) = 0;
  351. MP_DIGITS(&sy) = 0;
  352. MP_DIGITS(&sz) = 0;
  353. MP_CHECKOK(mp_init(&sx));
  354. MP_CHECKOK(mp_init(&sy));
  355. MP_CHECKOK(mp_init(&sz));
  356. /* if n = 0 then r = inf */
  357. if (mp_cmp_z(n) == 0) {
  358. mp_zero(rx);
  359. mp_zero(ry);
  360. res = MP_OKAY;
  361. goto CLEANUP;
  362. /* if n < 0 then out of range error */
  363. } else if (mp_cmp_z(n) < 0) {
  364. res = MP_RANGE;
  365. goto CLEANUP;
  366. }
  367. /* Convert from integer to floating point */
  368. ecfp_i2fp(p.x, px, ecgroup);
  369. ecfp_i2fp(p.y, py, ecgroup);
  370. ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
  371. /* Init r to point at infinity */
  372. for (i = 0; i < group->numDoubles; i++) {
  373. r.z[i] = 0;
  374. }
  375. /* double and add method */
  376. l = mpl_significant_bits(n) - 1;
  377. for (i = l; i >= 0; i--) {
  378. /* R = 2R */
  379. group->pt_dbl_jac(&r, &r, group);
  380. /* if n_i = 1, then R = R + Q */
  381. if (MP_GET_BIT(n, i) != 0) {
  382. group->pt_add_jac_aff(&r, &p, &r, group);
  383. }
  384. }
  385. /* Convert from floating point to integer */
  386. ecfp_fp2i(&sx, r.x, ecgroup);
  387. ecfp_fp2i(&sy, r.y, ecgroup);
  388. ecfp_fp2i(&sz, r.z, ecgroup);
  389. /* convert result R to affine coordinates */
  390. MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
  391. CLEANUP:
  392. mp_clear(&sx);
  393. mp_clear(&sy);
  394. mp_clear(&sz);
  395. return res;
  396. }
  397. /* Computes R = nP where R is (rx, ry) and P is the base point. Elliptic
  398. * curve points P and R can be identical. Uses mixed Modified-Jacobian
  399. * co-ordinates for doubling and Chudnovsky Jacobian coordinates for
  400. * additions. Uses 5-bit window NAF method (algorithm 11) for scalar-point
  401. * multiplication from Brown, Hankerson, Lopez, Menezes. Software
  402. * Implementation of the NIST Elliptic Curves Over Prime Fields. */
  403. mp_err
  404. ec_GFp_point_mul_wNAF_fp(const mp_int *n, const mp_int *px,
  405. const mp_int *py, mp_int *rx, mp_int *ry,
  406. const ECGroup *ecgroup)
  407. {
  408. mp_err res = MP_OKAY;
  409. mp_int sx, sy, sz;
  410. EC_group_fp *group = (EC_group_fp *) ecgroup->extra1;
  411. ecfp_chud_pt precomp[16];
  412. ecfp_aff_pt p;
  413. ecfp_jm_pt r;
  414. signed char naf[group->orderBitSize + 1];
  415. int i;
  416. MP_DIGITS(&sx) = 0;
  417. MP_DIGITS(&sy) = 0;
  418. MP_DIGITS(&sz) = 0;
  419. MP_CHECKOK(mp_init(&sx));
  420. MP_CHECKOK(mp_init(&sy));
  421. MP_CHECKOK(mp_init(&sz));
  422. /* if n = 0 then r = inf */
  423. if (mp_cmp_z(n) == 0) {
  424. mp_zero(rx);
  425. mp_zero(ry);
  426. res = MP_OKAY;
  427. goto CLEANUP;
  428. /* if n < 0 then out of range error */
  429. } else if (mp_cmp_z(n) < 0) {
  430. res = MP_RANGE;
  431. goto CLEANUP;
  432. }
  433. /* Convert from integer to floating point */
  434. ecfp_i2fp(p.x, px, ecgroup);
  435. ecfp_i2fp(p.y, py, ecgroup);
  436. ecfp_i2fp(group->curvea, &(ecgroup->curvea), ecgroup);
  437. /* Perform precomputation */
  438. group->precompute_chud(precomp, &p, group);
  439. /* Compute 5NAF */
  440. ec_compute_wNAF(naf, group->orderBitSize, n, 5);
  441. /* Init R = pt at infinity */
  442. for (i = 0; i < group->numDoubles; i++) {
  443. r.z[i] = 0;
  444. }
  445. /* wNAF method */
  446. for (i = group->orderBitSize; i >= 0; i--) {
  447. /* R = 2R */
  448. group->pt_dbl_jm(&r, &r, group);
  449. if (naf[i] != 0) {
  450. group->pt_add_jm_chud(&r, &precomp[(naf[i] + 15) / 2], &r,
  451. group);
  452. }
  453. }
  454. /* Convert from floating point to integer */
  455. ecfp_fp2i(&sx, r.x, ecgroup);
  456. ecfp_fp2i(&sy, r.y, ecgroup);
  457. ecfp_fp2i(&sz, r.z, ecgroup);
  458. /* convert result R to affine coordinates */
  459. MP_CHECKOK(ec_GFp_pt_jac2aff(&sx, &sy, &sz, rx, ry, ecgroup));
  460. CLEANUP:
  461. mp_clear(&sx);
  462. mp_clear(&sy);
  463. mp_clear(&sz);
  464. return res;
  465. }
  466. /* Cleans up extra memory allocated in ECGroup for this implementation. */
  467. void
  468. ec_GFp_extra_free_fp(ECGroup *group)
  469. {
  470. if (group->extra1 != NULL) {
  471. free(group->extra1);
  472. group->extra1 = NULL;
  473. }
  474. }
  475. /* Tests what precision floating point arithmetic is set to. This should
  476. * be either a 53-bit mantissa (IEEE standard) or a 64-bit mantissa
  477. * (extended precision on x86) and sets it into the EC_group_fp. Returns
  478. * either 53 or 64 accordingly. */
  479. int
  480. ec_set_fp_precision(EC_group_fp * group)
  481. {
  482. double a = 9007199254740992.0; /* 2^53 */
  483. double b = a + 1;
  484. if (a == b) {
  485. group->fpPrecision = 53;
  486. group->alpha = ecfp_alpha_53;
  487. return 53;
  488. }
  489. group->fpPrecision = 64;
  490. group->alpha = ecfp_alpha_64;
  491. return 64;
  492. }