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

Wed, 31 Dec 2014 06:09:35 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Wed, 31 Dec 2014 06:09:35 +0100
changeset 0
6474c204b198
permissions
-rw-r--r--

Cloned upstream origin tor-browser at tor-browser-31.3.0esr-4.5-1-build1
revision ID fc1c9ff7c1b2defdbc039f12214767608f46423f for hacking purpose.

michael@0 1 /* This Source Code Form is subject to the terms of the Mozilla Public
michael@0 2 * License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0 3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
michael@0 4
michael@0 5 /* This source file is meant to be included by other source files
michael@0 6 * (ecp_fp###.c, where ### is one of 160, 192, 224) and should not
michael@0 7 * constitute an independent compilation unit. It requires the following
michael@0 8 * preprocessor definitions be made: ECFP_BSIZE - the number of bits in
michael@0 9 * the field's prime
michael@0 10 * ECFP_NUMDOUBLES - the number of doubles to store one
michael@0 11 * multi-precision integer in floating point
michael@0 12
michael@0 13 /* Adds a prefix to a given token to give a unique token name. Prefixes
michael@0 14 * with "ecfp" + ECFP_BSIZE + "_". e.g. if ECFP_BSIZE = 160, then
michael@0 15 * PREFIX(hello) = ecfp160_hello This optimization allows static function
michael@0 16 * linking and compiler loop unrolling without code duplication. */
michael@0 17 #ifndef PREFIX
michael@0 18 #define PREFIX(b) PREFIX1(ECFP_BSIZE, b)
michael@0 19 #define PREFIX1(bsize, b) PREFIX2(bsize, b)
michael@0 20 #define PREFIX2(bsize, b) ecfp ## bsize ## _ ## b
michael@0 21 #endif
michael@0 22
michael@0 23 /* Returns true iff every double in d is 0. (If d == 0 and it is tidied,
michael@0 24 * this will be true.) */
michael@0 25 mp_err PREFIX(isZero) (const double *d) {
michael@0 26 int i;
michael@0 27
michael@0 28 for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0 29 if (d[i] != 0)
michael@0 30 return MP_NO;
michael@0 31 }
michael@0 32 return MP_YES;
michael@0 33 }
michael@0 34
michael@0 35 /* Sets the multi-precision floating point number at t = 0 */
michael@0 36 void PREFIX(zero) (double *t) {
michael@0 37 int i;
michael@0 38
michael@0 39 for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0 40 t[i] = 0;
michael@0 41 }
michael@0 42 }
michael@0 43
michael@0 44 /* Sets the multi-precision floating point number at t = 1 */
michael@0 45 void PREFIX(one) (double *t) {
michael@0 46 int i;
michael@0 47
michael@0 48 t[0] = 1;
michael@0 49 for (i = 1; i < ECFP_NUMDOUBLES; i++) {
michael@0 50 t[i] = 0;
michael@0 51 }
michael@0 52 }
michael@0 53
michael@0 54 /* Checks if point P(x, y, z) is at infinity. Uses Jacobian coordinates. */
michael@0 55 mp_err PREFIX(pt_is_inf_jac) (const ecfp_jac_pt * p) {
michael@0 56 return PREFIX(isZero) (p->z);
michael@0 57 }
michael@0 58
michael@0 59 /* Sets the Jacobian point P to be at infinity. */
michael@0 60 void PREFIX(set_pt_inf_jac) (ecfp_jac_pt * p) {
michael@0 61 PREFIX(zero) (p->z);
michael@0 62 }
michael@0 63
michael@0 64 /* Checks if point P(x, y) is at infinity. Uses Affine coordinates. */
michael@0 65 mp_err PREFIX(pt_is_inf_aff) (const ecfp_aff_pt * p) {
michael@0 66 if (PREFIX(isZero) (p->x) == MP_YES && PREFIX(isZero) (p->y) == MP_YES)
michael@0 67 return MP_YES;
michael@0 68 return MP_NO;
michael@0 69 }
michael@0 70
michael@0 71 /* Sets the affine point P to be at infinity. */
michael@0 72 void PREFIX(set_pt_inf_aff) (ecfp_aff_pt * p) {
michael@0 73 PREFIX(zero) (p->x);
michael@0 74 PREFIX(zero) (p->y);
michael@0 75 }
michael@0 76
michael@0 77 /* Checks if point P(x, y, z, a*z^4) is at infinity. Uses Modified
michael@0 78 * Jacobian coordinates. */
michael@0 79 mp_err PREFIX(pt_is_inf_jm) (const ecfp_jm_pt * p) {
michael@0 80 return PREFIX(isZero) (p->z);
michael@0 81 }
michael@0 82
michael@0 83 /* Sets the Modified Jacobian point P to be at infinity. */
michael@0 84 void PREFIX(set_pt_inf_jm) (ecfp_jm_pt * p) {
michael@0 85 PREFIX(zero) (p->z);
michael@0 86 }
michael@0 87
michael@0 88 /* Checks if point P(x, y, z, z^2, z^3) is at infinity. Uses Chudnovsky
michael@0 89 * Jacobian coordinates */
michael@0 90 mp_err PREFIX(pt_is_inf_chud) (const ecfp_chud_pt * p) {
michael@0 91 return PREFIX(isZero) (p->z);
michael@0 92 }
michael@0 93
michael@0 94 /* Sets the Chudnovsky Jacobian point P to be at infinity. */
michael@0 95 void PREFIX(set_pt_inf_chud) (ecfp_chud_pt * p) {
michael@0 96 PREFIX(zero) (p->z);
michael@0 97 }
michael@0 98
michael@0 99 /* Copies a multi-precision floating point number, Setting dest = src */
michael@0 100 void PREFIX(copy) (double *dest, const double *src) {
michael@0 101 int i;
michael@0 102
michael@0 103 for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0 104 dest[i] = src[i];
michael@0 105 }
michael@0 106 }
michael@0 107
michael@0 108 /* Sets dest = -src */
michael@0 109 void PREFIX(negLong) (double *dest, const double *src) {
michael@0 110 int i;
michael@0 111
michael@0 112 for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
michael@0 113 dest[i] = -src[i];
michael@0 114 }
michael@0 115 }
michael@0 116
michael@0 117 /* Sets r = -p p = (x, y, z, z2, z3) r = (x, -y, z, z2, z3) Uses
michael@0 118 * Chudnovsky Jacobian coordinates. */
michael@0 119 /* TODO reverse order */
michael@0 120 void PREFIX(pt_neg_chud) (const ecfp_chud_pt * p, ecfp_chud_pt * r) {
michael@0 121 int i;
michael@0 122
michael@0 123 PREFIX(copy) (r->x, p->x);
michael@0 124 PREFIX(copy) (r->z, p->z);
michael@0 125 PREFIX(copy) (r->z2, p->z2);
michael@0 126 PREFIX(copy) (r->z3, p->z3);
michael@0 127 for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0 128 r->y[i] = -p->y[i];
michael@0 129 }
michael@0 130 }
michael@0 131
michael@0 132 /* Computes r = x + y. Does not tidy or reduce. Any combinations of r, x,
michael@0 133 * y can point to the same data. Componentwise adds first ECFP_NUMDOUBLES
michael@0 134 * doubles of x and y and stores the result in r. */
michael@0 135 void PREFIX(addShort) (double *r, const double *x, const double *y) {
michael@0 136 int i;
michael@0 137
michael@0 138 for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0 139 *r++ = *x++ + *y++;
michael@0 140 }
michael@0 141 }
michael@0 142
michael@0 143 /* Computes r = x + y. Does not tidy or reduce. Any combinations of r, x,
michael@0 144 * y can point to the same data. Componentwise adds first
michael@0 145 * 2*ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
michael@0 146 void PREFIX(addLong) (double *r, const double *x, const double *y) {
michael@0 147 int i;
michael@0 148
michael@0 149 for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
michael@0 150 *r++ = *x++ + *y++;
michael@0 151 }
michael@0 152 }
michael@0 153
michael@0 154 /* Computes r = x - y. Does not tidy or reduce. Any combinations of r, x,
michael@0 155 * y can point to the same data. Componentwise subtracts first
michael@0 156 * ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
michael@0 157 void PREFIX(subtractShort) (double *r, const double *x, const double *y) {
michael@0 158 int i;
michael@0 159
michael@0 160 for (i = 0; i < ECFP_NUMDOUBLES; i++) {
michael@0 161 *r++ = *x++ - *y++;
michael@0 162 }
michael@0 163 }
michael@0 164
michael@0 165 /* Computes r = x - y. Does not tidy or reduce. Any combinations of r, x,
michael@0 166 * y can point to the same data. Componentwise subtracts first
michael@0 167 * 2*ECFP_NUMDOUBLES doubles of x and y and stores the result in r. */
michael@0 168 void PREFIX(subtractLong) (double *r, const double *x, const double *y) {
michael@0 169 int i;
michael@0 170
michael@0 171 for (i = 0; i < 2 * ECFP_NUMDOUBLES; i++) {
michael@0 172 *r++ = *x++ - *y++;
michael@0 173 }
michael@0 174 }
michael@0 175
michael@0 176 /* Computes r = x*y. Both x and y should be tidied and reduced,
michael@0 177 * r must be different (point to different memory) than x and y.
michael@0 178 * Does not tidy or reduce. */
michael@0 179 void PREFIX(multiply)(double *r, const double *x, const double *y) {
michael@0 180 int i, j;
michael@0 181
michael@0 182 for(j=0;j<ECFP_NUMDOUBLES-1;j++) {
michael@0 183 r[j] = x[0] * y[j];
michael@0 184 r[j+(ECFP_NUMDOUBLES-1)] = x[ECFP_NUMDOUBLES-1] * y[j];
michael@0 185 }
michael@0 186 r[ECFP_NUMDOUBLES-1] = x[0] * y[ECFP_NUMDOUBLES-1];
michael@0 187 r[ECFP_NUMDOUBLES-1] += x[ECFP_NUMDOUBLES-1] * y[0];
michael@0 188 r[2*ECFP_NUMDOUBLES-2] = x[ECFP_NUMDOUBLES-1] * y[ECFP_NUMDOUBLES-1];
michael@0 189 r[2*ECFP_NUMDOUBLES-1] = 0;
michael@0 190
michael@0 191 for(i=1;i<ECFP_NUMDOUBLES-1;i++) {
michael@0 192 for(j=0;j<ECFP_NUMDOUBLES;j++) {
michael@0 193 r[i+j] += (x[i] * y[j]);
michael@0 194 }
michael@0 195 }
michael@0 196 }
michael@0 197
michael@0 198 /* Computes the square of x and stores the result in r. x should be
michael@0 199 * tidied & reduced, r will be neither tidied nor reduced.
michael@0 200 * r should point to different memory than x */
michael@0 201 void PREFIX(square) (double *r, const double *x) {
michael@0 202 PREFIX(multiply) (r, x, x);
michael@0 203 }
michael@0 204
michael@0 205 /* Perform a point doubling in Jacobian coordinates. Input and output
michael@0 206 * should be multi-precision floating point integers. */
michael@0 207 void PREFIX(pt_dbl_jac) (const ecfp_jac_pt * dp, ecfp_jac_pt * dr,
michael@0 208 const EC_group_fp * group) {
michael@0 209 double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0 210 M[2 * ECFP_NUMDOUBLES], S[2 * ECFP_NUMDOUBLES];
michael@0 211
michael@0 212 /* Check for point at infinity */
michael@0 213 if (PREFIX(pt_is_inf_jac) (dp) == MP_YES) {
michael@0 214 /* Set r = pt at infinity */
michael@0 215 PREFIX(set_pt_inf_jac) (dr);
michael@0 216 goto CLEANUP;
michael@0 217 }
michael@0 218
michael@0 219 /* Perform typical point doubling operations */
michael@0 220
michael@0 221 /* TODO? is it worthwhile to do optimizations for when pz = 1? */
michael@0 222
michael@0 223 if (group->aIsM3) {
michael@0 224 /* When a = -3, M = 3(px - pz^2)(px + pz^2) */
michael@0 225 PREFIX(square) (t1, dp->z);
michael@0 226 group->ecfp_reduce(t1, t1, group); /* 2^23 since the negative
michael@0 227 * rounding buys another bit */
michael@0 228 PREFIX(addShort) (t0, dp->x, t1); /* 2*2^23 */
michael@0 229 PREFIX(subtractShort) (t1, dp->x, t1); /* 2 * 2^23 */
michael@0 230 PREFIX(multiply) (M, t0, t1); /* 40 * 2^46 */
michael@0 231 PREFIX(addLong) (t0, M, M); /* 80 * 2^46 */
michael@0 232 PREFIX(addLong) (M, t0, M); /* 120 * 2^46 < 2^53 */
michael@0 233 group->ecfp_reduce(M, M, group);
michael@0 234 } else {
michael@0 235 /* Generic case */
michael@0 236 /* M = 3 (px^2) + a*(pz^4) */
michael@0 237 PREFIX(square) (t0, dp->x);
michael@0 238 PREFIX(addLong) (M, t0, t0);
michael@0 239 PREFIX(addLong) (t0, t0, M); /* t0 = 3(px^2) */
michael@0 240 PREFIX(square) (M, dp->z);
michael@0 241 group->ecfp_reduce(M, M, group);
michael@0 242 PREFIX(square) (t1, M);
michael@0 243 group->ecfp_reduce(t1, t1, group);
michael@0 244 PREFIX(multiply) (M, t1, group->curvea); /* M = a(pz^4) */
michael@0 245 PREFIX(addLong) (M, M, t0);
michael@0 246 group->ecfp_reduce(M, M, group);
michael@0 247 }
michael@0 248
michael@0 249 /* rz = 2 * py * pz */
michael@0 250 PREFIX(multiply) (t1, dp->y, dp->z);
michael@0 251 PREFIX(addLong) (t1, t1, t1);
michael@0 252 group->ecfp_reduce(dr->z, t1, group);
michael@0 253
michael@0 254 /* t0 = 2y^2 */
michael@0 255 PREFIX(square) (t0, dp->y);
michael@0 256 group->ecfp_reduce(t0, t0, group);
michael@0 257 PREFIX(addShort) (t0, t0, t0);
michael@0 258
michael@0 259 /* S = 4 * px * py^2 = 2 * px * t0 */
michael@0 260 PREFIX(multiply) (S, dp->x, t0);
michael@0 261 PREFIX(addLong) (S, S, S);
michael@0 262 group->ecfp_reduce(S, S, group);
michael@0 263
michael@0 264 /* rx = M^2 - 2 * S */
michael@0 265 PREFIX(square) (t1, M);
michael@0 266 PREFIX(subtractShort) (t1, t1, S);
michael@0 267 PREFIX(subtractShort) (t1, t1, S);
michael@0 268 group->ecfp_reduce(dr->x, t1, group);
michael@0 269
michael@0 270 /* ry = M * (S - rx) - 8 * py^4 */
michael@0 271 PREFIX(square) (t1, t0); /* t1 = 4y^4 */
michael@0 272 PREFIX(subtractShort) (S, S, dr->x);
michael@0 273 PREFIX(multiply) (t0, M, S);
michael@0 274 PREFIX(subtractLong) (t0, t0, t1);
michael@0 275 PREFIX(subtractLong) (t0, t0, t1);
michael@0 276 group->ecfp_reduce(dr->y, t0, group);
michael@0 277
michael@0 278 CLEANUP:
michael@0 279 return;
michael@0 280 }
michael@0 281
michael@0 282 /* Perform a point addition using coordinate system Jacobian + Affine ->
michael@0 283 * Jacobian. Input and output should be multi-precision floating point
michael@0 284 * integers. */
michael@0 285 void PREFIX(pt_add_jac_aff) (const ecfp_jac_pt * p, const ecfp_aff_pt * q,
michael@0 286 ecfp_jac_pt * r, const EC_group_fp * group) {
michael@0 287 /* Temporary storage */
michael@0 288 double A[2 * ECFP_NUMDOUBLES], B[2 * ECFP_NUMDOUBLES],
michael@0 289 C[2 * ECFP_NUMDOUBLES], C2[2 * ECFP_NUMDOUBLES],
michael@0 290 D[2 * ECFP_NUMDOUBLES], C3[2 * ECFP_NUMDOUBLES];
michael@0 291
michael@0 292 /* Check for point at infinity for p or q */
michael@0 293 if (PREFIX(pt_is_inf_aff) (q) == MP_YES) {
michael@0 294 PREFIX(copy) (r->x, p->x);
michael@0 295 PREFIX(copy) (r->y, p->y);
michael@0 296 PREFIX(copy) (r->z, p->z);
michael@0 297 goto CLEANUP;
michael@0 298 } else if (PREFIX(pt_is_inf_jac) (p) == MP_YES) {
michael@0 299 PREFIX(copy) (r->x, q->x);
michael@0 300 PREFIX(copy) (r->y, q->y);
michael@0 301 /* Since the affine point is not infinity, we can set r->z = 1 */
michael@0 302 PREFIX(one) (r->z);
michael@0 303 goto CLEANUP;
michael@0 304 }
michael@0 305
michael@0 306 /* Calculates c = qx * pz^2 - px d = (qy * b - py) rx = d^2 - c^3 + 2
michael@0 307 * (px * c^2) ry = d * (c-rx) - py*c^3 rz = c * pz */
michael@0 308
michael@0 309 /* A = pz^2, B = pz^3 */
michael@0 310 PREFIX(square) (A, p->z);
michael@0 311 group->ecfp_reduce(A, A, group);
michael@0 312 PREFIX(multiply) (B, A, p->z);
michael@0 313 group->ecfp_reduce(B, B, group);
michael@0 314
michael@0 315 /* C = qx * A - px */
michael@0 316 PREFIX(multiply) (C, q->x, A);
michael@0 317 PREFIX(subtractShort) (C, C, p->x);
michael@0 318 group->ecfp_reduce(C, C, group);
michael@0 319
michael@0 320 /* D = qy * B - py */
michael@0 321 PREFIX(multiply) (D, q->y, B);
michael@0 322 PREFIX(subtractShort) (D, D, p->y);
michael@0 323 group->ecfp_reduce(D, D, group);
michael@0 324
michael@0 325 /* C2 = C^2, C3 = C^3 */
michael@0 326 PREFIX(square) (C2, C);
michael@0 327 group->ecfp_reduce(C2, C2, group);
michael@0 328 PREFIX(multiply) (C3, C2, C);
michael@0 329 group->ecfp_reduce(C3, C3, group);
michael@0 330
michael@0 331 /* rz = A = pz * C */
michael@0 332 PREFIX(multiply) (A, p->z, C);
michael@0 333 group->ecfp_reduce(r->z, A, group);
michael@0 334
michael@0 335 /* C = px * C^2, untidied, unreduced */
michael@0 336 PREFIX(multiply) (C, p->x, C2);
michael@0 337
michael@0 338 /* A = D^2, untidied, unreduced */
michael@0 339 PREFIX(square) (A, D);
michael@0 340
michael@0 341 /* rx = B = A - C3 - C - C = D^2 - (C^3 + 2 * (px * C^2) */
michael@0 342 PREFIX(subtractShort) (A, A, C3);
michael@0 343 PREFIX(subtractLong) (A, A, C);
michael@0 344 PREFIX(subtractLong) (A, A, C);
michael@0 345 group->ecfp_reduce(r->x, A, group);
michael@0 346
michael@0 347 /* B = py * C3, untidied, unreduced */
michael@0 348 PREFIX(multiply) (B, p->y, C3);
michael@0 349
michael@0 350 /* C = px * C^2 - rx */
michael@0 351 PREFIX(subtractShort) (C, C, r->x);
michael@0 352 group->ecfp_reduce(C, C, group);
michael@0 353
michael@0 354 /* ry = A = D * C - py * C^3 */
michael@0 355 PREFIX(multiply) (A, D, C);
michael@0 356 PREFIX(subtractLong) (A, A, B);
michael@0 357 group->ecfp_reduce(r->y, A, group);
michael@0 358
michael@0 359 CLEANUP:
michael@0 360 return;
michael@0 361 }
michael@0 362
michael@0 363 /* Perform a point addition using Jacobian coordinate system. Input and
michael@0 364 * output should be multi-precision floating point integers. */
michael@0 365 void PREFIX(pt_add_jac) (const ecfp_jac_pt * p, const ecfp_jac_pt * q,
michael@0 366 ecfp_jac_pt * r, const EC_group_fp * group) {
michael@0 367
michael@0 368 /* Temporary Storage */
michael@0 369 double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0 370 U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
michael@0 371 S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
michael@0 372 H3[2 * ECFP_NUMDOUBLES];
michael@0 373
michael@0 374 /* Check for point at infinity for p, if so set r = q */
michael@0 375 if (PREFIX(pt_is_inf_jac) (p) == MP_YES) {
michael@0 376 PREFIX(copy) (r->x, q->x);
michael@0 377 PREFIX(copy) (r->y, q->y);
michael@0 378 PREFIX(copy) (r->z, q->z);
michael@0 379 goto CLEANUP;
michael@0 380 }
michael@0 381
michael@0 382 /* Check for point at infinity for p, if so set r = q */
michael@0 383 if (PREFIX(pt_is_inf_jac) (q) == MP_YES) {
michael@0 384 PREFIX(copy) (r->x, p->x);
michael@0 385 PREFIX(copy) (r->y, p->y);
michael@0 386 PREFIX(copy) (r->z, p->z);
michael@0 387 goto CLEANUP;
michael@0 388 }
michael@0 389
michael@0 390 /* U = px * qz^2 , S = py * qz^3 */
michael@0 391 PREFIX(square) (t0, q->z);
michael@0 392 group->ecfp_reduce(t0, t0, group);
michael@0 393 PREFIX(multiply) (U, p->x, t0);
michael@0 394 group->ecfp_reduce(U, U, group);
michael@0 395 PREFIX(multiply) (t1, t0, q->z);
michael@0 396 group->ecfp_reduce(t1, t1, group);
michael@0 397 PREFIX(multiply) (t0, p->y, t1);
michael@0 398 group->ecfp_reduce(S, t0, group);
michael@0 399
michael@0 400 /* H = qx*(pz)^2 - U , R = (qy * pz^3 - S) */
michael@0 401 PREFIX(square) (t0, p->z);
michael@0 402 group->ecfp_reduce(t0, t0, group);
michael@0 403 PREFIX(multiply) (H, q->x, t0);
michael@0 404 PREFIX(subtractShort) (H, H, U);
michael@0 405 group->ecfp_reduce(H, H, group);
michael@0 406 PREFIX(multiply) (t1, t0, p->z); /* t1 = pz^3 */
michael@0 407 group->ecfp_reduce(t1, t1, group);
michael@0 408 PREFIX(multiply) (t0, t1, q->y); /* t0 = qy * pz^3 */
michael@0 409 PREFIX(subtractShort) (t0, t0, S);
michael@0 410 group->ecfp_reduce(R, t0, group);
michael@0 411
michael@0 412 /* U = U*H^2, H3 = H^3 */
michael@0 413 PREFIX(square) (t0, H);
michael@0 414 group->ecfp_reduce(t0, t0, group);
michael@0 415 PREFIX(multiply) (t1, U, t0);
michael@0 416 group->ecfp_reduce(U, t1, group);
michael@0 417 PREFIX(multiply) (H3, t0, H);
michael@0 418 group->ecfp_reduce(H3, H3, group);
michael@0 419
michael@0 420 /* rz = pz * qz * H */
michael@0 421 PREFIX(multiply) (t0, q->z, H);
michael@0 422 group->ecfp_reduce(t0, t0, group);
michael@0 423 PREFIX(multiply) (t1, t0, p->z);
michael@0 424 group->ecfp_reduce(r->z, t1, group);
michael@0 425
michael@0 426 /* rx = R^2 - H^3 - 2 * U */
michael@0 427 PREFIX(square) (t0, R);
michael@0 428 PREFIX(subtractShort) (t0, t0, H3);
michael@0 429 PREFIX(subtractShort) (t0, t0, U);
michael@0 430 PREFIX(subtractShort) (t0, t0, U);
michael@0 431 group->ecfp_reduce(r->x, t0, group);
michael@0 432
michael@0 433 /* ry = R(U - rx) - S*H3 */
michael@0 434 PREFIX(subtractShort) (t1, U, r->x);
michael@0 435 PREFIX(multiply) (t0, t1, R);
michael@0 436 PREFIX(multiply) (t1, S, H3);
michael@0 437 PREFIX(subtractLong) (t1, t0, t1);
michael@0 438 group->ecfp_reduce(r->y, t1, group);
michael@0 439
michael@0 440 CLEANUP:
michael@0 441 return;
michael@0 442 }
michael@0 443
michael@0 444 /* Perform a point doubling in Modified Jacobian coordinates. Input and
michael@0 445 * output should be multi-precision floating point integers. */
michael@0 446 void PREFIX(pt_dbl_jm) (const ecfp_jm_pt * p, ecfp_jm_pt * r,
michael@0 447 const EC_group_fp * group) {
michael@0 448
michael@0 449 /* Temporary storage */
michael@0 450 double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0 451 M[2 * ECFP_NUMDOUBLES], S[2 * ECFP_NUMDOUBLES],
michael@0 452 U[2 * ECFP_NUMDOUBLES], T[2 * ECFP_NUMDOUBLES];
michael@0 453
michael@0 454 /* Check for point at infinity */
michael@0 455 if (PREFIX(pt_is_inf_jm) (p) == MP_YES) {
michael@0 456 /* Set r = pt at infinity by setting rz = 0 */
michael@0 457 PREFIX(set_pt_inf_jm) (r);
michael@0 458 goto CLEANUP;
michael@0 459 }
michael@0 460
michael@0 461 /* M = 3 (px^2) + a*(pz^4) */
michael@0 462 PREFIX(square) (t0, p->x);
michael@0 463 PREFIX(addLong) (M, t0, t0);
michael@0 464 PREFIX(addLong) (t0, t0, M); /* t0 = 3(px^2) */
michael@0 465 PREFIX(addShort) (t0, t0, p->az4);
michael@0 466 group->ecfp_reduce(M, t0, group);
michael@0 467
michael@0 468 /* rz = 2 * py * pz */
michael@0 469 PREFIX(multiply) (t1, p->y, p->z);
michael@0 470 PREFIX(addLong) (t1, t1, t1);
michael@0 471 group->ecfp_reduce(r->z, t1, group);
michael@0 472
michael@0 473 /* t0 = 2y^2, U = 8y^4 */
michael@0 474 PREFIX(square) (t0, p->y);
michael@0 475 group->ecfp_reduce(t0, t0, group);
michael@0 476 PREFIX(addShort) (t0, t0, t0);
michael@0 477 PREFIX(square) (U, t0);
michael@0 478 group->ecfp_reduce(U, U, group);
michael@0 479 PREFIX(addShort) (U, U, U);
michael@0 480
michael@0 481 /* S = 4 * px * py^2 = 2 * px * t0 */
michael@0 482 PREFIX(multiply) (S, p->x, t0);
michael@0 483 group->ecfp_reduce(S, S, group);
michael@0 484 PREFIX(addShort) (S, S, S);
michael@0 485
michael@0 486 /* rx = M^2 - 2S */
michael@0 487 PREFIX(square) (T, M);
michael@0 488 PREFIX(subtractShort) (T, T, S);
michael@0 489 PREFIX(subtractShort) (T, T, S);
michael@0 490 group->ecfp_reduce(r->x, T, group);
michael@0 491
michael@0 492 /* ry = M * (S - rx) - U */
michael@0 493 PREFIX(subtractShort) (S, S, r->x);
michael@0 494 PREFIX(multiply) (t0, M, S);
michael@0 495 PREFIX(subtractShort) (t0, t0, U);
michael@0 496 group->ecfp_reduce(r->y, t0, group);
michael@0 497
michael@0 498 /* ra*z^4 = 2*U*(apz4) */
michael@0 499 PREFIX(multiply) (t1, U, p->az4);
michael@0 500 PREFIX(addLong) (t1, t1, t1);
michael@0 501 group->ecfp_reduce(r->az4, t1, group);
michael@0 502
michael@0 503 CLEANUP:
michael@0 504 return;
michael@0 505 }
michael@0 506
michael@0 507 /* Perform a point doubling using coordinates Affine -> Chudnovsky
michael@0 508 * Jacobian. Input and output should be multi-precision floating point
michael@0 509 * integers. */
michael@0 510 void PREFIX(pt_dbl_aff2chud) (const ecfp_aff_pt * p, ecfp_chud_pt * r,
michael@0 511 const EC_group_fp * group) {
michael@0 512 double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0 513 M[2 * ECFP_NUMDOUBLES], twoY2[2 * ECFP_NUMDOUBLES],
michael@0 514 S[2 * ECFP_NUMDOUBLES];
michael@0 515
michael@0 516 /* Check for point at infinity for p, if so set r = O */
michael@0 517 if (PREFIX(pt_is_inf_aff) (p) == MP_YES) {
michael@0 518 PREFIX(set_pt_inf_chud) (r);
michael@0 519 goto CLEANUP;
michael@0 520 }
michael@0 521
michael@0 522 /* M = 3(px)^2 + a */
michael@0 523 PREFIX(square) (t0, p->x);
michael@0 524 PREFIX(addLong) (t1, t0, t0);
michael@0 525 PREFIX(addLong) (t1, t1, t0);
michael@0 526 PREFIX(addShort) (t1, t1, group->curvea);
michael@0 527 group->ecfp_reduce(M, t1, group);
michael@0 528
michael@0 529 /* twoY2 = 2*(py)^2, S = 4(px)(py)^2 */
michael@0 530 PREFIX(square) (twoY2, p->y);
michael@0 531 PREFIX(addLong) (twoY2, twoY2, twoY2);
michael@0 532 group->ecfp_reduce(twoY2, twoY2, group);
michael@0 533 PREFIX(multiply) (S, p->x, twoY2);
michael@0 534 PREFIX(addLong) (S, S, S);
michael@0 535 group->ecfp_reduce(S, S, group);
michael@0 536
michael@0 537 /* rx = M^2 - 2S */
michael@0 538 PREFIX(square) (t0, M);
michael@0 539 PREFIX(subtractShort) (t0, t0, S);
michael@0 540 PREFIX(subtractShort) (t0, t0, S);
michael@0 541 group->ecfp_reduce(r->x, t0, group);
michael@0 542
michael@0 543 /* ry = M(S-rx) - 8y^4 */
michael@0 544 PREFIX(subtractShort) (t0, S, r->x);
michael@0 545 PREFIX(multiply) (t1, t0, M);
michael@0 546 PREFIX(square) (t0, twoY2);
michael@0 547 PREFIX(subtractLong) (t1, t1, t0);
michael@0 548 PREFIX(subtractLong) (t1, t1, t0);
michael@0 549 group->ecfp_reduce(r->y, t1, group);
michael@0 550
michael@0 551 /* rz = 2py */
michael@0 552 PREFIX(addShort) (r->z, p->y, p->y);
michael@0 553
michael@0 554 /* rz2 = rz^2 */
michael@0 555 PREFIX(square) (t0, r->z);
michael@0 556 group->ecfp_reduce(r->z2, t0, group);
michael@0 557
michael@0 558 /* rz3 = rz^3 */
michael@0 559 PREFIX(multiply) (t0, r->z, r->z2);
michael@0 560 group->ecfp_reduce(r->z3, t0, group);
michael@0 561
michael@0 562 CLEANUP:
michael@0 563 return;
michael@0 564 }
michael@0 565
michael@0 566 /* Perform a point addition using coordinates: Modified Jacobian +
michael@0 567 * Chudnovsky Jacobian -> Modified Jacobian. Input and output should be
michael@0 568 * multi-precision floating point integers. */
michael@0 569 void PREFIX(pt_add_jm_chud) (ecfp_jm_pt * p, ecfp_chud_pt * q,
michael@0 570 ecfp_jm_pt * r, const EC_group_fp * group) {
michael@0 571
michael@0 572 double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0 573 U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
michael@0 574 S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
michael@0 575 H3[2 * ECFP_NUMDOUBLES], pz2[2 * ECFP_NUMDOUBLES];
michael@0 576
michael@0 577 /* Check for point at infinity for p, if so set r = q need to convert
michael@0 578 * from Chudnovsky form to Modified Jacobian form */
michael@0 579 if (PREFIX(pt_is_inf_jm) (p) == MP_YES) {
michael@0 580 PREFIX(copy) (r->x, q->x);
michael@0 581 PREFIX(copy) (r->y, q->y);
michael@0 582 PREFIX(copy) (r->z, q->z);
michael@0 583 PREFIX(square) (t0, q->z2);
michael@0 584 group->ecfp_reduce(t0, t0, group);
michael@0 585 PREFIX(multiply) (t1, t0, group->curvea);
michael@0 586 group->ecfp_reduce(r->az4, t1, group);
michael@0 587 goto CLEANUP;
michael@0 588 }
michael@0 589 /* Check for point at infinity for q, if so set r = p */
michael@0 590 if (PREFIX(pt_is_inf_chud) (q) == MP_YES) {
michael@0 591 PREFIX(copy) (r->x, p->x);
michael@0 592 PREFIX(copy) (r->y, p->y);
michael@0 593 PREFIX(copy) (r->z, p->z);
michael@0 594 PREFIX(copy) (r->az4, p->az4);
michael@0 595 goto CLEANUP;
michael@0 596 }
michael@0 597
michael@0 598 /* U = px * qz^2 */
michael@0 599 PREFIX(multiply) (U, p->x, q->z2);
michael@0 600 group->ecfp_reduce(U, U, group);
michael@0 601
michael@0 602 /* H = qx*(pz)^2 - U */
michael@0 603 PREFIX(square) (t0, p->z);
michael@0 604 group->ecfp_reduce(pz2, t0, group);
michael@0 605 PREFIX(multiply) (H, pz2, q->x);
michael@0 606 group->ecfp_reduce(H, H, group);
michael@0 607 PREFIX(subtractShort) (H, H, U);
michael@0 608
michael@0 609 /* U = U*H^2, H3 = H^3 */
michael@0 610 PREFIX(square) (t0, H);
michael@0 611 group->ecfp_reduce(t0, t0, group);
michael@0 612 PREFIX(multiply) (t1, U, t0);
michael@0 613 group->ecfp_reduce(U, t1, group);
michael@0 614 PREFIX(multiply) (H3, t0, H);
michael@0 615 group->ecfp_reduce(H3, H3, group);
michael@0 616
michael@0 617 /* S = py * qz^3 */
michael@0 618 PREFIX(multiply) (S, p->y, q->z3);
michael@0 619 group->ecfp_reduce(S, S, group);
michael@0 620
michael@0 621 /* R = (qy * z1^3 - s) */
michael@0 622 PREFIX(multiply) (t0, pz2, p->z);
michael@0 623 group->ecfp_reduce(t0, t0, group);
michael@0 624 PREFIX(multiply) (R, t0, q->y);
michael@0 625 PREFIX(subtractShort) (R, R, S);
michael@0 626 group->ecfp_reduce(R, R, group);
michael@0 627
michael@0 628 /* rz = pz * qz * H */
michael@0 629 PREFIX(multiply) (t1, q->z, H);
michael@0 630 group->ecfp_reduce(t1, t1, group);
michael@0 631 PREFIX(multiply) (t0, p->z, t1);
michael@0 632 group->ecfp_reduce(r->z, t0, group);
michael@0 633
michael@0 634 /* rx = R^2 - H^3 - 2 * U */
michael@0 635 PREFIX(square) (t0, R);
michael@0 636 PREFIX(subtractShort) (t0, t0, H3);
michael@0 637 PREFIX(subtractShort) (t0, t0, U);
michael@0 638 PREFIX(subtractShort) (t0, t0, U);
michael@0 639 group->ecfp_reduce(r->x, t0, group);
michael@0 640
michael@0 641 /* ry = R(U - rx) - S*H3 */
michael@0 642 PREFIX(subtractShort) (t1, U, r->x);
michael@0 643 PREFIX(multiply) (t0, t1, R);
michael@0 644 PREFIX(multiply) (t1, S, H3);
michael@0 645 PREFIX(subtractLong) (t1, t0, t1);
michael@0 646 group->ecfp_reduce(r->y, t1, group);
michael@0 647
michael@0 648 if (group->aIsM3) { /* a == -3 */
michael@0 649 /* a(rz^4) = -3 * ((rz^2)^2) */
michael@0 650 PREFIX(square) (t0, r->z);
michael@0 651 group->ecfp_reduce(t0, t0, group);
michael@0 652 PREFIX(square) (t1, t0);
michael@0 653 PREFIX(addLong) (t0, t1, t1);
michael@0 654 PREFIX(addLong) (t0, t0, t1);
michael@0 655 PREFIX(negLong) (t0, t0);
michael@0 656 group->ecfp_reduce(r->az4, t0, group);
michael@0 657 } else { /* Generic case */
michael@0 658 /* a(rz^4) = a * ((rz^2)^2) */
michael@0 659 PREFIX(square) (t0, r->z);
michael@0 660 group->ecfp_reduce(t0, t0, group);
michael@0 661 PREFIX(square) (t1, t0);
michael@0 662 group->ecfp_reduce(t1, t1, group);
michael@0 663 PREFIX(multiply) (t0, group->curvea, t1);
michael@0 664 group->ecfp_reduce(r->az4, t0, group);
michael@0 665 }
michael@0 666 CLEANUP:
michael@0 667 return;
michael@0 668 }
michael@0 669
michael@0 670 /* Perform a point addition using Chudnovsky Jacobian coordinates. Input
michael@0 671 * and output should be multi-precision floating point integers. */
michael@0 672 void PREFIX(pt_add_chud) (const ecfp_chud_pt * p, const ecfp_chud_pt * q,
michael@0 673 ecfp_chud_pt * r, const EC_group_fp * group) {
michael@0 674
michael@0 675 /* Temporary Storage */
michael@0 676 double t0[2 * ECFP_NUMDOUBLES], t1[2 * ECFP_NUMDOUBLES],
michael@0 677 U[2 * ECFP_NUMDOUBLES], R[2 * ECFP_NUMDOUBLES],
michael@0 678 S[2 * ECFP_NUMDOUBLES], H[2 * ECFP_NUMDOUBLES],
michael@0 679 H3[2 * ECFP_NUMDOUBLES];
michael@0 680
michael@0 681 /* Check for point at infinity for p, if so set r = q */
michael@0 682 if (PREFIX(pt_is_inf_chud) (p) == MP_YES) {
michael@0 683 PREFIX(copy) (r->x, q->x);
michael@0 684 PREFIX(copy) (r->y, q->y);
michael@0 685 PREFIX(copy) (r->z, q->z);
michael@0 686 PREFIX(copy) (r->z2, q->z2);
michael@0 687 PREFIX(copy) (r->z3, q->z3);
michael@0 688 goto CLEANUP;
michael@0 689 }
michael@0 690
michael@0 691 /* Check for point at infinity for p, if so set r = q */
michael@0 692 if (PREFIX(pt_is_inf_chud) (q) == MP_YES) {
michael@0 693 PREFIX(copy) (r->x, p->x);
michael@0 694 PREFIX(copy) (r->y, p->y);
michael@0 695 PREFIX(copy) (r->z, p->z);
michael@0 696 PREFIX(copy) (r->z2, p->z2);
michael@0 697 PREFIX(copy) (r->z3, p->z3);
michael@0 698 goto CLEANUP;
michael@0 699 }
michael@0 700
michael@0 701 /* U = px * qz^2 */
michael@0 702 PREFIX(multiply) (U, p->x, q->z2);
michael@0 703 group->ecfp_reduce(U, U, group);
michael@0 704
michael@0 705 /* H = qx*(pz)^2 - U */
michael@0 706 PREFIX(multiply) (H, q->x, p->z2);
michael@0 707 PREFIX(subtractShort) (H, H, U);
michael@0 708 group->ecfp_reduce(H, H, group);
michael@0 709
michael@0 710 /* U = U*H^2, H3 = H^3 */
michael@0 711 PREFIX(square) (t0, H);
michael@0 712 group->ecfp_reduce(t0, t0, group);
michael@0 713 PREFIX(multiply) (t1, U, t0);
michael@0 714 group->ecfp_reduce(U, t1, group);
michael@0 715 PREFIX(multiply) (H3, t0, H);
michael@0 716 group->ecfp_reduce(H3, H3, group);
michael@0 717
michael@0 718 /* S = py * qz^3 */
michael@0 719 PREFIX(multiply) (S, p->y, q->z3);
michael@0 720 group->ecfp_reduce(S, S, group);
michael@0 721
michael@0 722 /* rz = pz * qz * H */
michael@0 723 PREFIX(multiply) (t0, q->z, H);
michael@0 724 group->ecfp_reduce(t0, t0, group);
michael@0 725 PREFIX(multiply) (t1, t0, p->z);
michael@0 726 group->ecfp_reduce(r->z, t1, group);
michael@0 727
michael@0 728 /* R = (qy * z1^3 - s) */
michael@0 729 PREFIX(multiply) (t0, q->y, p->z3);
michael@0 730 PREFIX(subtractShort) (t0, t0, S);
michael@0 731 group->ecfp_reduce(R, t0, group);
michael@0 732
michael@0 733 /* rx = R^2 - H^3 - 2 * U */
michael@0 734 PREFIX(square) (t0, R);
michael@0 735 PREFIX(subtractShort) (t0, t0, H3);
michael@0 736 PREFIX(subtractShort) (t0, t0, U);
michael@0 737 PREFIX(subtractShort) (t0, t0, U);
michael@0 738 group->ecfp_reduce(r->x, t0, group);
michael@0 739
michael@0 740 /* ry = R(U - rx) - S*H3 */
michael@0 741 PREFIX(subtractShort) (t1, U, r->x);
michael@0 742 PREFIX(multiply) (t0, t1, R);
michael@0 743 PREFIX(multiply) (t1, S, H3);
michael@0 744 PREFIX(subtractLong) (t1, t0, t1);
michael@0 745 group->ecfp_reduce(r->y, t1, group);
michael@0 746
michael@0 747 /* rz2 = rz^2 */
michael@0 748 PREFIX(square) (t0, r->z);
michael@0 749 group->ecfp_reduce(r->z2, t0, group);
michael@0 750
michael@0 751 /* rz3 = rz^3 */
michael@0 752 PREFIX(multiply) (t0, r->z, r->z2);
michael@0 753 group->ecfp_reduce(r->z3, t0, group);
michael@0 754
michael@0 755 CLEANUP:
michael@0 756 return;
michael@0 757 }
michael@0 758
michael@0 759 /* Expects out to be an array of size 16 of Chudnovsky Jacobian points.
michael@0 760 * Fills in Chudnovsky Jacobian form (x, y, z, z^2, z^3), for -15P, -13P,
michael@0 761 * -11P, -9P, -7P, -5P, -3P, -P, P, 3P, 5P, 7P, 9P, 11P, 13P, 15P */
michael@0 762 void PREFIX(precompute_chud) (ecfp_chud_pt * out, const ecfp_aff_pt * p,
michael@0 763 const EC_group_fp * group) {
michael@0 764
michael@0 765 ecfp_chud_pt p2;
michael@0 766
michael@0 767 /* Set out[8] = P */
michael@0 768 PREFIX(copy) (out[8].x, p->x);
michael@0 769 PREFIX(copy) (out[8].y, p->y);
michael@0 770 PREFIX(one) (out[8].z);
michael@0 771 PREFIX(one) (out[8].z2);
michael@0 772 PREFIX(one) (out[8].z3);
michael@0 773
michael@0 774 /* Set p2 = 2P */
michael@0 775 PREFIX(pt_dbl_aff2chud) (p, &p2, group);
michael@0 776
michael@0 777 /* Set 3P, 5P, ..., 15P */
michael@0 778 PREFIX(pt_add_chud) (&out[8], &p2, &out[9], group);
michael@0 779 PREFIX(pt_add_chud) (&out[9], &p2, &out[10], group);
michael@0 780 PREFIX(pt_add_chud) (&out[10], &p2, &out[11], group);
michael@0 781 PREFIX(pt_add_chud) (&out[11], &p2, &out[12], group);
michael@0 782 PREFIX(pt_add_chud) (&out[12], &p2, &out[13], group);
michael@0 783 PREFIX(pt_add_chud) (&out[13], &p2, &out[14], group);
michael@0 784 PREFIX(pt_add_chud) (&out[14], &p2, &out[15], group);
michael@0 785
michael@0 786 /* Set -15P, -13P, ..., -P */
michael@0 787 PREFIX(pt_neg_chud) (&out[8], &out[7]);
michael@0 788 PREFIX(pt_neg_chud) (&out[9], &out[6]);
michael@0 789 PREFIX(pt_neg_chud) (&out[10], &out[5]);
michael@0 790 PREFIX(pt_neg_chud) (&out[11], &out[4]);
michael@0 791 PREFIX(pt_neg_chud) (&out[12], &out[3]);
michael@0 792 PREFIX(pt_neg_chud) (&out[13], &out[2]);
michael@0 793 PREFIX(pt_neg_chud) (&out[14], &out[1]);
michael@0 794 PREFIX(pt_neg_chud) (&out[15], &out[0]);
michael@0 795 }
michael@0 796
michael@0 797 /* Expects out to be an array of size 16 of Jacobian points. Fills in
michael@0 798 * Jacobian form (x, y, z), for O, P, 2P, ... 15P */
michael@0 799 void PREFIX(precompute_jac) (ecfp_jac_pt * precomp, const ecfp_aff_pt * p,
michael@0 800 const EC_group_fp * group) {
michael@0 801 int i;
michael@0 802
michael@0 803 /* fill precomputation table */
michael@0 804 /* set precomp[0] */
michael@0 805 PREFIX(set_pt_inf_jac) (&precomp[0]);
michael@0 806 /* set precomp[1] */
michael@0 807 PREFIX(copy) (precomp[1].x, p->x);
michael@0 808 PREFIX(copy) (precomp[1].y, p->y);
michael@0 809 if (PREFIX(pt_is_inf_aff) (p) == MP_YES) {
michael@0 810 PREFIX(zero) (precomp[1].z);
michael@0 811 } else {
michael@0 812 PREFIX(one) (precomp[1].z);
michael@0 813 }
michael@0 814 /* set precomp[2] */
michael@0 815 group->pt_dbl_jac(&precomp[1], &precomp[2], group);
michael@0 816
michael@0 817 /* set rest of precomp */
michael@0 818 for (i = 3; i < 16; i++) {
michael@0 819 group->pt_add_jac_aff(&precomp[i - 1], p, &precomp[i], group);
michael@0 820 }
michael@0 821 }

mercurial