crypto/ec/ec_point.c

   1 /*
   2  *
   3  *      ec_point.c
   4  *
   5  *      Elliptic Curve Arithmetic Functions
   6  *
   7  *      Copyright (C) Lenka Fibikova 2000
   8  *
   9  *
  10  */
  11
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <assert.h>
  15 #include <memory.h>
  16
  17 #include <openssl/bn.h>
  18
  19 #include "bn_modfs.h"
  20 #include "bn_mont2.h"
  21 #include "ec.h"
  22
  23 EC_POINT *ECP_new()
  24 {
  25         EC_POINT *ret;
  26
  27         ret=(EC_POINT *)malloc(sizeof(EC_POINT));
  28         if (ret == NULL) return NULL;
  29         ret->X = BN_new();
  30         ret->Y = BN_new();
  31         ret->Z = BN_new();
  32         ret->is_in_mont = 0;
  33
  34         if (ret->X == NULL || ret->Y == NULL || ret->Z == NULL)
  35         {
  36                 if (ret->X != NULL) BN_free(ret->X);
  37                 if (ret->Y != NULL) BN_free(ret->Y);
  38                 if (ret->Z != NULL) BN_free(ret->Z);
  39                 free(ret);
  40                 return(NULL);
  41         }
  42         return(ret);
  43 }
  44
  45 void ECP_clear_free(EC_POINT *P)
  46 {
  47         if (P == NULL) return;
  48
  49         P->is_in_mont = 0;
  50         if (P->X != NULL) BN_clear_free(P->X);
  51         if (P->Y != NULL) BN_clear_free(P->Y);
  52         if (P->Z != NULL) BN_clear_free(P->Z);
  53         free(P);
  54 }
  55
  56 void ECP_clear_free_precompute(ECP_PRECOMPUTE *prec)
  57 {
  58         int i;
  59         int max;
  60
  61         if (prec == NULL) return;
  62         if (prec->Pi != NULL)
  63         {
  64                 max = 1;
  65                 max <<= (prec->r - 1);
  66
  67                 for (i = 0; i < max; i++)
  68                 {
  69                         if (prec->Pi[i] != NULL) ECP_clear_free(prec->Pi[i]);
  70                 }
  71         }
  72         free(prec);
  73 }
  74
  75 int ECP_is_on_ec(EC_POINT *P, EC *E, BN_CTX *ctx)
  76 {
  77         BIGNUM *n0, *n1, *n2, *p;
  78         int Pnorm;
  79
  80         assert(P != NULL);
  81         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
  82
  83         assert(E != NULL);
  84         assert(E->A != NULL && E->B != NULL && E->p != NULL);
  85
  86         assert(ctx != NULL);
  87
  88         assert(!P->is_in_mont);
  89
  90         if (ECP_is_infty(P)) return 1;
  91
  92         n0 = ctx->bn[ctx->tos];
  93         n1 = ctx->bn[ctx->tos + 1];
  94         n2 = ctx->bn[ctx->tos + 2];
  95         ctx->tos += 3;
  96
  97
  98         p = E->p;
  99
 100         Pnorm = (ECP_is_norm(P));
 101
 102         if (!Pnorm)
 103         {
 104                 if (!BN_mod_mul(n0, P->Z, P->Z, p, ctx)) goto err;
 105                 if (!BN_mod_mul(n1, n0, n0, p, ctx)) goto err;
 106                 if (!BN_mod_mul(n2, n0, n1, p, ctx)) goto err;
 107         }
 108
 109         if (!BN_mod_mul(n0, P->X, P->X, p, ctx)) goto err;
 110         if (!BN_mod_mul(n0, n0, P->X, p, ctx)) goto err;
 111
 112         if (Pnorm)
 113         {
 114                 if (!BN_mod_mul(n1, P->X, E->A, p, ctx)) goto err;
 115         }
 116         else
 117         {
 118                 if (!BN_mod_mul(n1, n1, P->X, p, ctx)) goto err;
 119                 if (!BN_mod_mul(n1, n1, E->A, p, ctx)) goto err;
 120         }
 121         if (!BN_mod_add(n0, n0, n1, p, ctx)) goto err;
 122
 123         if (Pnorm)
 124         {
 125                 if (!BN_mod_add(n0, n0, E->B, p, ctx)) goto err;
 126         }
 127         else
 128         {
 129                 if (!BN_mod_mul(n2, n2, E->B,  p, ctx)) goto err;
 130                 if (!BN_mod_add(n0, n0, n2, p, ctx)) goto err;
 131         }
 132
 133         if (!BN_mod_mul(n1, P->Y, P->Y, p, ctx)) goto err;
 134
 135         if (BN_cmp(n0, n1))
 136         {
 137                 ctx->tos -= 3;
 138                 return 0;
 139         }
 140
 141         ctx->tos -= 3;
 142         return 1;
 143
 144 err:
 145         ctx->tos -= 3;
 146         return -1;
 147 }
 148
 149
 150 EC_POINT *ECP_generate(BIGNUM *x, BIGNUM *z,EC *E, BN_CTX *ctx)
 151 /* x == NULL || z = 0  -> point of infinity     */
 152 /* z == NULL || z = 1  -> normalized            */
 153 {
 154         BIGNUM *n0, *n1;
 155         EC_POINT *ret;
 156         int Pnorm, Pinfty, X0, A0;
 157
 158         assert(E != NULL);
 159         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
 160
 161         assert(ctx != NULL);
 162
 163         Pinfty = (x == NULL);
 164         Pnorm = (z == NULL);
 165         if (!Pnorm)
 166         {
 167                 Pnorm = BN_is_one(z);
 168                 Pinfty = (Pinfty || BN_is_zero(z));
 169         }
 170
 171         if (Pinfty)
 172         {
 173                 if ((ret = ECP_new()) == NULL) return NULL;
 174                 if (!BN_zero(ret->Z))
 175                 {
 176                         ECP_clear_free(ret);
 177                         return NULL;
 178                 }
 179                 return ret;
 180         }
 181
 182         X0 = BN_is_zero(x);
 183         A0 = BN_is_zero(E->A);
 184
 185         if ((ret = ECP_new()) == NULL) return NULL;
 186
 187         ret->is_in_mont = 0;
 188
 189         n0 = ctx->bn[ctx->tos];
 190         n1 = ctx->bn[ctx->tos + 1];
 191         if (!BN_zero(n0)) return NULL;
 192         if (!BN_zero(n1)) return NULL;
 193
 194         ctx->tos += 2;
 195
 196         if (!X0)
 197         {
 198                 if (!BN_mod_sqr(n0, x, E->p, ctx)) goto err;
 199                 if (!BN_mod_mul(n0, n0, x, E->p, ctx)) goto err;        /* x^3 */
 200         }
 201
 202         if (!X0 && !A0)
 203         {
 204                 if (!BN_mod_mul(n1, E->A, x, E->p, ctx)) goto err;      /* Ax */
 205                 if (!BN_mod_add(n0, n0, n1, E->p, ctx)) goto err;       /* x^3 + Ax */
 206         }
 207
 208         if (!BN_is_zero(E->B))
 209                 if (!BN_mod_add(n0, n0, E->B, E->p, ctx)) goto err;     /* x^3 + Ax +B */
 210
 211         if (!BN_mod_sqrt(ret->Y, n0, E->p, ctx)) goto err;
 212         if (BN_copy(ret->X, x) == NULL) goto err;
 213
 214         if (Pnorm)
 215         {
 216                 if (!BN_one(ret->Z)) goto err;
 217         }
 218         else
 219         {
 220                 if (BN_copy(ret->Z, z) == NULL) goto err;
 221                 if (!BN_mod_sqr(n0, z, E->p, ctx)) goto err;
 222                 if (!BN_mod_mul(ret->X, ret->X, n0, E->p, ctx)) goto err;
 223                 if (!BN_mod_mul(n0, n0, z, E->p, ctx)) goto err;
 224                 if (!BN_mod_mul(ret->Y, ret->Y, n0, E->p, ctx)) goto err;
 225         }
 226
 227 #ifdef TEST
 228         if (!ECP_is_on_ec(ret, E, ctx)) goto err;
 229 #endif
 230
 231         ctx->tos -= 2;
 232         return ret;
 233
 234 err:
 235         if (ret != NULL) ECP_clear_free(ret);
 236         ctx->tos -= 2;
 237         return NULL;
 238 }
 239
 240 int ECP_ecp2bin(EC_POINT *P, unsigned char *to, int form)
 241 /* form =       1 ... compressed
 242                         2 ... uncompressed
 243                         3 ... hybrid */
 244 {
 245         int bytes, bx, by;
 246
 247         assert (P != NULL);
 248         assert (P->X != NULL && P->Y != NULL && P->Z != NULL);
 249         assert (!P->is_in_mont);
 250         assert (ECP_is_norm(P) || ECP_is_infty(P));
 251         assert (to != NULL);
 252         assert (form > 0 && form < 4);
 253
 254         if (BN_is_zero(P->Z))
 255         {
 256                 to[0] = 0;
 257                 return 1;
 258         }
 259
 260         bx = BN_num_bytes(P->X);
 261         if (form == 1 ) bytes = bx + 1;
 262         else
 263         {
 264                 by = BN_num_bytes(P->Y);
 265                 bytes = (bx > by ? bx : by);
 266                 bytes = bytes * 2 + 1;
 267         }
 268         memset(to, 0, bytes);
 269
 270         switch (form)
 271         {
 272         case 1: to[0] = 2;      break;
 273         case 2: to[0] = 4;      break;
 274         case 3: to[0] = 6;      break;
 275         }
 276         if (form != 2) to[0] += BN_is_bit_set(P->Y, 0);
 277
 278
 279         if ((BN_bn2bin(P->X, to + 1)) != bx) return 0;
 280         if (form != 1)
 281         {
 282                 if ((BN_bn2bin(P->Y, to + bx + 1)) != by) return 0;
 283         }
 284
 285         return bytes;
 286 }
 287
 288 int ECP_bin2ecp(unsigned char *from, int len, EC_POINT *P, EC *E, BN_CTX *ctx)
 289 {
 290         int y;
 291         BIGNUM *x;
 292         EC_POINT *pp;
 293
 294         assert (E != NULL);
 295         assert (E->A != NULL && E->B != NULL && E->p != NULL);
 296         assert (!E->is_in_mont);
 297
 298         assert (ctx != NULL);
 299         assert (from != NULL);
 300         assert (P != NULL);
 301         assert (P->X != NULL && P->Y != NULL && P->Z != NULL);
 302
 303         if (len == 1 && from[0] != 0) return 0;
 304
 305         if (len == 0 || len == 1)
 306         {
 307                 if (!BN_zero(P->Z)) return 0;
 308                 return 1;
 309         }
 310
 311         switch (from[0])
 312         {
 313         case 2:
 314         case 3:
 315                 y = from[0] - 2;
 316                 if ((x = BN_new()) == NULL) return 0;
 317                 if (BN_bin2bn(from + 1, len - 1, x) == NULL) return 0;
 318
 319                 pp = ECP_generate(x, NULL, E, ctx);
 320                 BN_clear_free(x);
 321                 if (pp == NULL) return 0;
 322
 323                 ECP_copy(P, pp);
 324                 ECP_clear_free(pp);
 325
 326                 if (BN_is_bit_set(P->Y, 0) != y)
 327                         if (!BN_sub(P->Y, E->p, P->Y)) return 0;
 328                 break;
 329
 330         case 4:
 331         case 6:
 332         case 7:
 333                 y = (len - 1)/2;
 334                 if (BN_bin2bn(from + 1, y, P->X) == NULL) return 0;
 335                 if (BN_bin2bn(from + y + 1, y, P->Y) == NULL) return 0;
 336                 if (!BN_set_word(P->Z, 1)) return 0;
 337                 break;
 338
 339         default:
 340                 assert(0);
 341
 342         }
 343
 344         if (!ECP_is_on_ec(P, E, ctx)) return 0;
 345         return 1;
 346 }
 347
 348 int ECP_normalize(EC_POINT *P, EC *E, BN_CTX *ctx)
 349 {
 350         BIGNUM *z, *zm;
 351
 352         assert (P != NULL);
 353         assert (P->X != NULL && P->Y != NULL && P->Z != NULL);
 354
 355         assert (E != NULL);
 356         assert (E->A != NULL && E->B != NULL && E->p != NULL);
 357
 358         assert (ctx != NULL);
 359
 360         if (ECP_is_norm(P)) return 1;
 361         if (ECP_is_infty(P)) return 0;
 362
 363         if ((zm = BN_mod_inverse(P->Z, P->Z, E->p, ctx)) == NULL) return 0;
 364
 365         assert(!P->is_in_mont);
 366
 367
 368         z = ctx->bn[ctx->tos];
 369         ctx->tos++;
 370
 371         if (!BN_mod_mul(z, zm, zm, E->p, ctx)) goto err;
 372         if (!BN_mod_mul(P->X, P->X, z, E->p, ctx)) goto err;
 373
 374         if (!BN_mod_mul(z, z, zm, E->p, ctx)) goto err;
 375         if (!BN_mod_mul(P->Y, P->Y, z, E->p, ctx)) goto err;
 376
 377         if (!BN_one(P->Z)) goto err;
 378
 379         if (zm != NULL) BN_clear_free(zm);
 380
 381         ctx->tos--;
 382         return 1;
 383
 384 err:
 385         if (zm != NULL) BN_clear_free(zm);
 386         ctx->tos--;
 387         return 0;
 388 }
 389
 390 int ECP_copy(EC_POINT *R, EC_POINT *P)
 391 {
 392         assert(P != NULL);
 393         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 394
 395         assert(R != NULL);
 396         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
 397
 398         if (BN_copy(R->X, P->X) == NULL) return 0;
 399         if (BN_copy(R->Y, P->Y) == NULL) return 0;
 400         if (BN_copy(R->Z, P->Z) == NULL) return 0;
 401         R->is_in_mont = P->is_in_mont;
 402
 403         return 1;
 404 }
 405
 406 EC_POINT *ECP_dup(EC_POINT *P)
 407 {
 408         EC_POINT *ret;
 409
 410         ret = ECP_new();
 411         if (ret == NULL) return NULL;
 412
 413         if (!ECP_copy(ret, P))
 414         {
 415                 ECP_clear_free(ret);
 416                 return(NULL);
 417         }
 418
 419         return(ret);
 420 }
 421
 422
 423 EC_POINT *ECP_minus(EC_POINT *P, BIGNUM *p) /* mont || non-mont */
 424 {
 425         EC_POINT *ret;
 426
 427         assert(P != NULL);
 428         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 429
 430         assert(p != NULL);
 431
 432         assert(BN_cmp(P->Y, p) < 0);
 433
 434         ret = ECP_dup(P);
 435         if (ret == NULL) return NULL;
 436
 437         if (BN_is_zero(ret->Y)) return ret;
 438
 439         if (!BN_sub(ret->Y, p, ret->Y))
 440         {
 441                 ECP_clear_free(ret);
 442                 return NULL;
 443         }
 444
 445         return ret;
 446 }
 447
 448
 449 #ifdef SIMPLE
 450 int ECP_cmp(EC_POINT *P, EC_POINT *Q, BIGNUM *p, BN_CTX *ctx)
 451 /* return values:
 452         -2 ... error
 453          0 ... P = Q
 454         -1 ... P = -Q
 455          1 ... else
 456 */
 457 {
 458         BIGNUM *n0, *n1, *n2, *n3, *n4;
 459         int Pnorm, Qnorm;
 460
 461         assert(P != NULL);
 462         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 463
 464         assert(Q != NULL);
 465         assert(Q->X != NULL && Q->Y != NULL && Q->Z != NULL);
 466
 467         assert(p != NULL);
 468         assert(ctx != NULL);
 469
 470         assert(!P->is_in_mont);
 471         assert(!Q->is_in_mont);
 472
 473         if (ECP_is_infty(P) && ECP_is_infty(Q)) return 0;
 474         if (ECP_is_infty(P) || ECP_is_infty(Q)) return 1;
 475
 476
 477         Pnorm = (ECP_is_norm(P));
 478         Qnorm = (ECP_is_norm(Q));
 479
 480         n0 = ctx->bn[ctx->tos];
 481         n1 = ctx->bn[ctx->tos + 1];
 482         n2 = ctx->bn[ctx->tos + 2];
 483         n3 = ctx->bn[ctx->tos + 3];
 484         n4 = ctx->bn[ctx->tos + 4];
 485         ctx->tos += 5;
 486
 487         if (Qnorm)
 488         {
 489                 if (BN_copy(n1, P->X) == NULL) goto err;                        /* L1 = x_p */
 490                 if (BN_copy(n2, P->Y) == NULL) goto err;                        /* L2 = y_p */
 491         }
 492         else
 493         {
 494                 if (!BN_sqr(n0, Q->Z, ctx)) goto err;
 495                 if (!BN_mod_mul(n1, P->X, n0, p, ctx)) goto err;        /* L1 = x_p * z_q^2 */
 496
 497                 if (!BN_mod_mul(n0, n0, Q->Z, p, ctx)) goto err;
 498                 if (!BN_mod_mul(n2, P->Y, n0, p, ctx)) goto err;        /* L2 = y_p * z_q^3 */
 499         }
 500
 501         if (Pnorm)
 502         {
 503                 if (BN_copy(n3, Q->X) == NULL) goto err;                        /* L3 = x_q */
 504                 if (BN_copy(n4, Q->Y) == NULL) goto err;                        /* L4 = y_q */
 505         }
 506         else
 507         {
 508                 if (!BN_sqr(n0, P->Z, ctx)) goto err;
 509                 if (!BN_mod_mul(n3, Q->X, n0, p, ctx)) goto err;        /* L3 = x_q * z_p^2 */
 510
 511                 if (!BN_mod_mul(n0, n0, P->Z, p, ctx)) goto err;
 512                 if (!BN_mod_mul(n4, Q->Y, n0, p, ctx)) goto err;        /* L4 = y_q * z_p^3 */
 513         }
 514
 515         if (!BN_mod_sub(n0, n1, n3, p, ctx)) goto err;                  /* L5 = L1 - L3 */
 516
 517         if (!BN_is_zero(n0))
 518         {
 519                 ctx->tos -= 5;
 520                 return 1;
 521         }
 522
 523         if (!BN_mod_sub(n0, n2, n4, p, ctx)) goto err;                  /* L6 = L2 - L4 */
 524
 525         if (!BN_is_zero(n0))
 526         {
 527                 ctx->tos -= 5;
 528                 return -1;
 529         }
 530
 531         ctx->tos -= 5;
 532         return 0;
 533
 534 err:
 535         ctx->tos -= 5;
 536         return -2;
 537 }
 538
 539 int ECP_double(EC_POINT *R, EC_POINT *P, EC *E, BN_CTX *ctx)
 540 /* R <- 2P (on E) */
 541 {
 542         BIGNUM *n0, *n1, *n2, *n3, *p;
 543         int Pnorm, A0;
 544
 545         assert(P != NULL);
 546         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 547
 548         assert(R != NULL);
 549         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
 550
 551         assert(E != NULL);
 552         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
 553
 554         assert(ctx != NULL);
 555
 556         assert(!P->is_in_mont);
 557
 558         if (ECP_is_infty(P))
 559         {
 560                 if (!BN_zero(R->Z)) return 0;
 561                 return 1;
 562         }
 563
 564         Pnorm = (ECP_is_norm(P));
 565         A0 = (BN_is_zero(E->A));
 566
 567         n0 = ctx->bn[ctx->tos];
 568         n1 = ctx->bn[ctx->tos + 1];
 569         n2 = ctx->bn[ctx->tos + 2];
 570         n3 = ctx->bn[ctx->tos + 3];
 571         ctx->tos += 4;
 572
 573         p = E->p;
 574
 575         /* L1 */
 576         if (Pnorm || A0)
 577         {
 578                 if (!BN_mod_sqr(n1, P->X, p, ctx)) goto err;
 579                 if (!BN_mul_word(n1, 3)) goto err;
 580                 if (!A0)                                   /* if A = 0: L1 = 3 * x^2 + a * z^4 = 3 * x ^2    */
 581                         if (!BN_mod_add(n1, n1, E->A, p, ctx)) goto err; /* L1 = 3 * x^2 + a * z^4 = 3 * x^2 + a */
 582         }
 583         else
 584         {
 585                 if (!BN_mod_sqr(n0, P->Z, p, ctx)) goto err;
 586                 if (!BN_mod_mul(n0, n0, n0, p, ctx)) goto err;
 587                 if (!BN_mod_mul(n0, n0, E->A, p, ctx)) goto err;
 588                 if (!BN_mod_sqr(n1, P->X, p, ctx)) goto err;
 589                 if (!BN_mul_word(n1, 3)) goto err;
 590                 if (!BN_mod_add(n1, n1, n0, p, ctx)) goto err;          /* L1 = 3 * x^2 + a * z^4 */
 591         }
 592
 593         /* Z */
 594         if (Pnorm)
 595         {
 596                 if (BN_copy(n0, P->Y) == NULL) goto err;
 597         }
 598         else
 599         {
 600                 if (!BN_mod_mul(n0, P->Y, P->Z, p, ctx)) goto err;
 601         }
 602         if (!BN_lshift1(n0, n0)) goto err;
 603         if (!BN_smod(R->Z, n0, p, ctx)) goto err;                               /* Z = 2 * y * z */
 604
 605         /* L2 */
 606         if (!BN_mod_sqr(n3, P->Y, p, ctx)) goto err;
 607         if (!BN_mod_mul(n2, P->X, n3, p, ctx)) goto err;
 608         if (!BN_lshift(n2, n2, 2)) goto err;
 609         if (!BN_smod(n2, n2, p, ctx)) goto err;                                 /* L2 = 4 * x * y^2 */
 610
 611         /* X */
 612         if (!BN_lshift1(n0, n2)) goto err;
 613         if (!BN_mod_sqr(R->X, n1, p, ctx)) goto err;
 614         if (!BN_mod_sub(R->X, R->X, n0, p, ctx)) goto err;              /* X = L1^2 - 2 * L2 */
 615
 616         /* L3 */
 617         if (!BN_mod_sqr(n0, n3, p, ctx)) goto err;
 618         if (!BN_lshift(n3, n0, 3)) goto err;
 619         if (!BN_smod(n3, n3, p, ctx)) goto err;                                 /* L3 = 8 * y^4 */
 620
 621         /* Y */
 622         if (!BN_mod_sub(n0, n2, R->X, p, ctx)) goto err;
 623         if (!BN_mod_mul(n0, n1, n0, p, ctx)) goto err;
 624         if (!BN_mod_sub(R->Y, n0, n3, p, ctx)) goto err;                /* Y = L1 * (L2 - X) - L3 */
 625
 626
 627 #ifdef TEST
 628         if (!ECP_is_on_ec(R, E, ctx)) return 0;
 629 #endif
 630
 631         ctx->tos -= 4;
 632         return 1;
 633
 634 err:
 635         ctx->tos -= 4;
 636         return 0;
 637 }
 638
 639 int ECP_add(EC_POINT *R, EC_POINT *P, EC_POINT *Q, EC *E, BN_CTX *ctx)
 640 /* R <- P + Q (on E) */
 641 {
 642         BIGNUM *n0, *n1, *n2, *n3, *n4, *n5, *n6, *p;
 643         int Pnorm, Qnorm;
 644
 645         assert(P != NULL);
 646         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 647
 648         assert(Q != NULL);
 649         assert(Q->X != NULL && Q->Y != NULL && Q->Z != NULL);
 650
 651         assert(R != NULL);
 652         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
 653
 654         assert(E != NULL);
 655         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
 656         assert(!BN_is_zero(E->h));;
 657
 658         assert(ctx != NULL);
 659
 660         assert(!P->is_in_mont);
 661         assert(!Q->is_in_mont);
 662
 663         if (P == Q) return ECP_double(R, P, E, ctx);
 664
 665         if (ECP_is_infty(P)) return ECP_copy(R, Q);
 666         if (ECP_is_infty(Q)) return ECP_copy(R, P);
 667
 668         Pnorm = (ECP_is_norm(P));
 669         Qnorm = (ECP_is_norm(Q));
 670
 671         n0 = ctx->bn[ctx->tos];
 672         n1 = ctx->bn[ctx->tos + 1];
 673         n2 = ctx->bn[ctx->tos + 2];
 674         n3 = ctx->bn[ctx->tos + 3];
 675         n4 = ctx->bn[ctx->tos + 4];
 676         n5 = ctx->bn[ctx->tos + 5];
 677         n6 = ctx->bn[ctx->tos + 6];
 678         ctx->tos += 7;
 679         p = E->p;
 680
 681         /* L1; L2 */
 682         if (Qnorm)
 683         {
 684                 if (BN_copy(n1, P->X) == NULL) goto err;         /* L1 = x_p */
 685                 if (BN_copy(n2, P->Y) == NULL) goto err;         /* L2 = y_p */
 686         }
 687         else
 688         {
 689                 if (!BN_sqr(n0, Q->Z, ctx)) goto err;
 690                 if (!BN_mod_mul(n1, P->X, n0, p, ctx)) goto err; /* L1 = x_p * z_q^2 */
 691
 692                 if (!BN_mod_mul(n0, n0, Q->Z, p, ctx)) goto err;
 693                 if (!BN_mod_mul(n2, P->Y, n0, p, ctx)) goto err; /* L2 = y_p * z_q^3 */
 694         }
 695
 696         /* L3; L4 */
 697         if (Pnorm)
 698         {
 699                 if (BN_copy(n3, Q->X) == NULL) goto err;         /* L3 = x_q */
 700                 if (BN_copy(n4, Q->Y) == NULL) goto err;         /* L4 = y_q */
 701         }
 702         else
 703         {
 704                 if (!BN_sqr(n0, P->Z, ctx)) goto err;
 705                 if (!BN_mod_mul(n3, Q->X, n0, p, ctx)) goto err; /* L3 = x_q * z_p^2 */
 706
 707                 if (!BN_mod_mul(n0, n0, P->Z, p, ctx)) goto err;
 708                 if (!BN_mod_mul(n4, Q->Y, n0, p, ctx)) goto err; /* L4 = y_q * z_p^3 */
 709         }
 710
 711         /* L5; L6 */
 712         if (!BN_mod_sub(n5, n1, n3, p, ctx)) goto err;          /* L5 = L1 - L3 */
 713         if (!BN_mod_sub(n6, n2, n4, p, ctx)) goto err;          /* L6 = L2 - L4 */
 714
 715         /* pata */
 716         if (BN_is_zero(n5))
 717         {
 718                 if (BN_is_zero(n6))     /* P = Q => P + Q = 2P */
 719                 {
 720                         ctx->tos -= 7;
 721                         return ECP_double(R, P, E, ctx);
 722                 }
 723                 else                             /* P = -Q => P + Q = \infty */
 724                 {
 725                         ctx->tos -= 7;
 726                         if (!BN_zero(R->Z)) return 0;
 727                         return 1;
 728                 }
 729         }
 730
 731         /* L7; L8 */
 732         if (!BN_mod_add(n1, n1, n3, p, ctx)) goto err;          /* L7 = L1 + L3 */
 733         if (!BN_mod_add(n2, n2, n4, p, ctx)) goto err;          /* L8 = L2 + L4 */
 734
 735         /* Z */
 736         if (Pnorm)
 737         {
 738                 if (BN_copy(n0, Q->Z) == NULL) goto err;
 739         }
 740         else
 741         {
 742                 if (!BN_mod_mul(n0, P->Z, Q->Z, p, ctx)) goto err;
 743         }
 744         if (!BN_mod_mul(R->Z, n0, n5, p, ctx)) goto err;        /* Z = z_p * z_q * L_5 */
 745
 746         /* X */
 747         if (!BN_mod_sqr(n0, n6, p, ctx)) goto err;
 748         if (!BN_mod_sqr(n4, n5, p, ctx)) goto err;
 749         if (!BN_mod_mul(n3, n1, n4, p, ctx)) goto err;
 750         if (!BN_mod_sub(R->X, n0, n3, p, ctx)) goto err;        /* X = L6^2 - L5^2 * L7 */
 751
 752         /* L9 */
 753         if (!BN_lshift1(n0, R->X)) goto err;
 754         if (!BN_mod_sub(n0, n3, n0, p, ctx)) goto err;          /* L9 = L5^2 * L7 - 2X */
 755
 756         /* Y */
 757         if (!BN_mod_mul(n0, n0, n6, p, ctx)) goto err;
 758         if (!BN_mod_mul(n5, n4, n5, p, ctx)) goto err;
 759         if (!BN_mod_mul(n1, n2, n5, p, ctx)) goto err;
 760         if (!BN_mod_sub(n0, n0, n1, p, ctx)) goto err;
 761         if (!BN_mod_mul(R->Y, n0, E->h, p, ctx)) goto err;      /* Y = (L6 * L9 - L8 * L5^3) / 2 */
 762
 763
 764
 765 #ifdef TEST
 766         if (!ECP_is_on_ec(R, E, ctx)) return 0;
 767 #endif
 768
 769         ctx->tos -= 7;
 770         return 1;
 771
 772 err:
 773         ctx->tos -= 7;
 774         return 0;
 775 }
 776
 777
 778 ECP_PRECOMPUTE *ECP_precompute(int r, EC_POINT *P, EC *E, BN_CTX *ctx)
 779 {
 780         ECP_PRECOMPUTE *ret;
 781         EC_POINT *P2;
 782         int i, max;
 783
 784         assert(r > 2);
 785         assert(!P->is_in_mont);
 786         assert(!E->is_in_mont);
 787
 788         ret=(ECP_PRECOMPUTE *)malloc(sizeof(ECP_PRECOMPUTE));
 789         if (ret == NULL) return NULL;
 790
 791         max = 1;
 792         max <<= (r - 1);
 793
 794         ret->r = 0;
 795
 796         ret->Pi=(EC_POINT **)malloc(sizeof(EC_POINT *) * max);
 797         if (ret->Pi == NULL) goto err;
 798
 799
 800         /* P2 = [2]P */
 801         if ((P2 = ECP_new()) == NULL) goto err;
 802         if (!ECP_double(P2, P, E, ctx)) goto err;
 803
 804         /* P_0 = P */
 805         if((ret->Pi[0] = ECP_dup(P)) == NULL) goto err;
 806
 807
 808         /* P_i = P_(i-1) + P2 */
 809         for (i = 1; i < max; i++)
 810         {
 811                 if ((ret->Pi[i] = ECP_new()) == NULL) goto err;
 812
 813                 if (!ECP_add(ret->Pi[i], P2, ret->Pi[i - 1], E, ctx)) goto err;
 814         }
 815
 816         ret->r = r;
 817         ECP_clear_free(P2);
 818
 819         return ret;
 820
 821 err:
 822         ECP_clear_free(P2);
 823         ECP_clear_free_precompute(ret);
 824         return NULL;
 825 }
 826
 827 int ECP_multiply(EC_POINT *R, BIGNUM *k, ECP_PRECOMPUTE *prec, EC *E, BN_CTX *ctx)
 828 /* R = [k]P */
 829 {
 830         int j;
 831         int t, nextw, h, r;
 832
 833         assert(R != NULL);
 834         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
 835
 836         assert(E != NULL);
 837         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
 838
 839         assert(k != NULL);
 840         assert(!k->neg);
 841
 842         assert(ctx != NULL);
 843         assert(prec != NULL);
 844
 845         assert(!E->is_in_mont);
 846
 847         if (BN_is_zero(k))
 848         {
 849                 if (!BN_zero(R->Z)) return 0;
 850                 R->is_in_mont = 0;
 851                 return 1;
 852         }
 853
 854
 855         j = BN_num_bits(k);
 856         j--;
 857
 858         r = prec->r;
 859
 860         if (!BN_zero(R->Z)) return 0;
 861         R->is_in_mont = 0;
 862
 863         while(j >= 0)
 864         {
 865                 if (!BN_is_bit_set(k, j))
 866                 {
 867                         if (!ECP_double(R, R, E, ctx)) return 0;
 868                         j--;
 869                 }
 870                 else
 871                 {
 872                         nextw = j - r;
 873                         if (nextw < -1) nextw = -1;
 874                         t = nextw + 1;
 875                         while(!BN_is_bit_set(k, t))
 876                         {
 877                                 t++;
 878                         }
 879
 880                         if (!ECP_double(R, R, E, ctx)) return 0;
 881
 882                         j--;
 883                         if (j < t) h = 0;
 884                         else
 885                         {
 886                                 h = 1;
 887                                 for(; j > t; j--)
 888                                 {
 889                                         h <<= 1;
 890                                         if (BN_is_bit_set(k, j)) h++;
 891                                         if (!ECP_double(R, R, E, ctx)) return 0;
 892                                 }
 893                                 if (!ECP_double(R, R, E, ctx)) return 0;
 894                                 j--;
 895                         }
 896
 897                         if (!ECP_add(R, R, prec->Pi[h], E, ctx)) return 0;
 898
 899                         for (; j > nextw; j--)
 900                         {
 901                                 if (!ECP_double(R, R, E, ctx)) return 0;
 902                         }
 903
 904                 }
 905         }
 906
 907         return 1;
 908 }
 909
 910 #endif /* SIMPLE */
 911
 912 #ifdef MONTGOMERY
 913
 914 int ECP_to_montgomery(EC_POINT *P, BN_MONTGOMERY *mont, BN_CTX *ctx)
 915 {
 916
 917         assert(P != NULL);
 918         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 919
 920         assert(mont != NULL);
 921         assert(mont->p != NULL);
 922
 923         assert(ctx != NULL);
 924
 925         if (P->is_in_mont) return 1;
 926
 927         if (!BN_lshift(P->X, P->X, mont->R_num_bits)) return 0;
 928         if (!BN_mod(P->X, P->X, mont->p, ctx)) return 0;
 929
 930         if (!BN_lshift(P->Y, P->Y, mont->R_num_bits)) return 0;
 931         if (!BN_mod(P->Y, P->Y, mont->p, ctx)) return 0;
 932
 933         if (!BN_lshift(P->Z, P->Z, mont->R_num_bits)) return 0;
 934         if (!BN_mod(P->Z, P->Z, mont->p, ctx)) return 0;
 935
 936         P->is_in_mont = 1;
 937         return 1;
 938 }
 939
 940
 941 int ECP_from_montgomery(EC_POINT *P, BN_MONTGOMERY *mont, BN_CTX *ctx)
 942 {
 943
 944         assert(P != NULL);
 945         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 946
 947         assert(mont != NULL);
 948         assert(mont->p != NULL);
 949
 950         assert(ctx != NULL);
 951
 952         if (!P->is_in_mont) return 1;
 953
 954         if (!BN_mont_red(P->X, mont, ctx)) return 0;
 955         if (!BN_mont_red(P->Y, mont, ctx)) return 0;
 956         if (!BN_mont_red(P->Z, mont, ctx)) return 0;
 957
 958         P->is_in_mont = 0;
 959         return 1;
 960 }
 961
 962 int ECP_mont_cmp(EC_POINT *P, EC_POINT *Q, BN_MONTGOMERY *mont, BN_CTX *ctx)
 963 /* return values:
 964         -2 ... error
 965          0 ... P = Q
 966         -1 ... P = -Q
 967          1 ... else
 968 */
 969 {
 970         BIGNUM *n0, *n1, *n2, *n3, *n4, *n5, *p;
 971
 972         assert(P != NULL);
 973         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
 974
 975         assert(Q != NULL);
 976         assert(Q->X != NULL && Q->Y != NULL && Q->Z != NULL);
 977
 978         assert(mont != NULL);
 979         assert(mont->p != NULL);
 980
 981         assert(ctx != NULL);
 982
 983         if (!P->is_in_mont)
 984                 if (!ECP_to_montgomery(P, mont, ctx)) return 0;
 985
 986         if (!Q->is_in_mont)
 987                 if (!ECP_to_montgomery(Q, mont, ctx)) return 0;
 988
 989
 990         if (ECP_is_infty(P) && ECP_is_infty(Q)) return 0;
 991         if (ECP_is_infty(P) || ECP_is_infty(Q)) return 1;
 992
 993
 994         n0 = ctx->bn[ctx->tos];
 995         n1 = ctx->bn[ctx->tos + 1];
 996         n2 = ctx->bn[ctx->tos + 2];
 997         n3 = ctx->bn[ctx->tos + 3];
 998         n4 = ctx->bn[ctx->tos + 4];
 999         n5 = ctx->bn[ctx->tos + 5];
1000         ctx->tos += 6;
1001
1002         p = mont->p;
1003
1004
1005         if (!BN_mont_mod_mul(n5, Q->Z, Q->Z, mont, ctx)) goto err;
1006         if (!BN_mont_mod_mul(n1, P->X, n5, mont, ctx)) goto err;        /* L1 = x_p * z_q^2 */
1007
1008         if (!BN_mont_mod_mul(n0, n5, Q->Z, mont, ctx)) goto err;
1009         if (!BN_mont_mod_mul(n2, P->Y, n0, mont, ctx)) goto err;        /* L2 = y_p * z_q^3 */
1010
1011         if (!BN_mont_mod_mul(n5, P->Z, P->Z, mont, ctx)) goto err;
1012         if (!BN_mont_mod_mul(n3, Q->X, n5, mont, ctx)) goto err;        /* L3 = x_q * z_p^2 */
1013
1014         if (!BN_mont_mod_mul(n0, n5, P->Z, mont, ctx)) goto err;
1015         if (!BN_mont_mod_mul(n4, Q->Y, n0, mont, ctx)) goto err;        /* L4 = y_q * z_p^3 */
1016
1017
1018         if (!BN_mod_sub_quick(n0, n1, n3, p)) goto err;                 /* L5 = L1 - L3 */
1019
1020         if (!BN_is_zero(n0))
1021         {
1022                 ctx->tos -= 6;
1023                 return 1;
1024         }
1025
1026         if (!BN_mod_sub_quick(n0, n2, n4, p)) goto err;                 /* L6 = L2 - L4 */
1027
1028         if (!BN_is_zero(n0))
1029         {
1030                 ctx->tos -= 6;
1031                 return -1;
1032         }
1033
1034         ctx->tos -= 6;
1035         return 0;
1036
1037 err:
1038         ctx->tos -= 6;
1039         return -2;
1040 }
1041
1042
1043 int ECP_mont_double(EC_POINT *R, EC_POINT *P, EC *E, BN_MONTGOMERY *mont, BN_CTX *ctx)
1044 /* R <- 2P (on E) */
1045 {
1046         BIGNUM *n0, *n1, *n2, *n3, *p;
1047
1048         assert(P != NULL);
1049         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
1050
1051         assert(R != NULL);
1052         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
1053
1054         assert(E != NULL);
1055         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
1056
1057         assert(ctx != NULL);
1058
1059         if (!P->is_in_mont)
1060                 if (!ECP_to_montgomery(P, mont, ctx)) return 0;
1061
1062         if (!E->is_in_mont)
1063                 if (!EC_to_montgomery(E, mont, ctx)) return 0;
1064
1065         R->is_in_mont = 1;
1066
1067         if (ECP_is_infty(P))
1068         {
1069                 if (!BN_zero(R->Z)) return 0;
1070                 return 1;
1071         }
1072
1073
1074         n0 = ctx->bn[ctx->tos];
1075         n1 = ctx->bn[ctx->tos + 1];
1076         n2 = ctx->bn[ctx->tos + 2];
1077         n3 = ctx->bn[ctx->tos + 3];
1078
1079         ctx->tos += 4;
1080
1081         p = E->p;
1082
1083         /* L1 */
1084         if (!BN_mont_mod_mul(n0, P->Z, P->Z, mont, ctx)) goto err;
1085         if (!BN_mont_mod_mul(n2, n0, n0, mont, ctx)) goto err;
1086         if (!BN_mont_mod_mul(n0, n2, E->A, mont, ctx)) goto err;
1087         if (!BN_mont_mod_mul(n1, P->X, P->X, mont, ctx)) goto err;
1088         if (!BN_mod_lshift1_quick(n2, n1, p)) goto err;
1089         if (!BN_mod_add_quick(n1, n1, n2, p)) goto err;
1090         if (!BN_mod_add_quick(n1, n1, n0, p)) goto err;         /* L1 = 3 * x^2 + a * z^4 */
1091
1092         /* Z */
1093         if (!BN_mont_mod_mul(n0, P->Y, P->Z, mont, ctx)) goto err;
1094         if (!BN_mod_lshift1_quick(R->Z, n0, p)) goto err;               /* Z = 2 * y * z */
1095
1096         /* L2 */
1097         if (!BN_mont_mod_mul(n3, P->Y, P->Y, mont, ctx)) goto err;
1098         if (!BN_mont_mod_mul(n2, P->X, n3, mont, ctx)) goto err;
1099         if (!BN_mod_lshift_quick(n2, n2, 2, p)) goto err;               /* L2 = 4 * x * y^2 */
1100
1101         /* X */
1102         if (!BN_mod_lshift1_quick(n0, n2, p)) goto err;
1103         if (!BN_mont_mod_mul(R->X, n1, n1, mont, ctx)) goto err;
1104         if (!BN_mod_sub_quick(R->X, R->X, n0, p)) goto err;     /* X = L1^2 - 2 * L2 */
1105
1106         /* L3 */
1107         if (!BN_mont_mod_mul(n0, n3, n3, mont, ctx)) goto err;
1108         if (!BN_mod_lshift_quick(n3, n0, 3, p)) goto err;               /* L3 = 8 * y^4 */
1109
1110
1111         /* Y */
1112         if (!BN_mod_sub_quick(n2, n2, R->X, p)) goto err;
1113         if (!BN_mont_mod_mul(n0, n1, n2, mont, ctx)) goto err;
1114         if (!BN_mod_sub_quick(R->Y, n0, n3, p)) goto err;               /* Y = L1 * (L2 - X) - L3 */
1115
1116         ctx->tos -= 4;
1117         return 1;
1118
1119 err:
1120         ctx->tos -= 4;
1121         return 0;
1122 }
1123
1124
1125 int ECP_mont_add(EC_POINT *R, EC_POINT *P, EC_POINT *Q, EC *E, BN_MONTGOMERY *mont, BN_CTX *ctx)
1126 /* R <- P + Q (on E) */
1127 {
1128         BIGNUM *n0, *n1, *n2, *n3, *n4, *n5, *n6, *p;
1129
1130         assert(P != NULL);
1131         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
1132
1133         assert(Q != NULL);
1134         assert(Q->X != NULL && Q->Y != NULL && Q->Z != NULL);
1135
1136         assert(R != NULL);
1137         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
1138
1139         assert(E != NULL);
1140         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
1141         assert(!BN_is_zero(E->h));;
1142
1143         assert(ctx != NULL);
1144
1145         if (!Q->is_in_mont)
1146                 if (!ECP_to_montgomery(Q, mont, ctx)) return 0;
1147
1148         if (!P->is_in_mont)
1149                 if (!ECP_to_montgomery(P, mont, ctx)) return 0;
1150
1151         if (!E->is_in_mont)
1152                 if (!EC_to_montgomery(E, mont, ctx)) return 0;
1153
1154         if (P == Q) return ECP_mont_double(R, P, E, mont, ctx);
1155
1156         if (ECP_is_infty(P)) return ECP_copy(R, Q);
1157         if (ECP_is_infty(Q)) return ECP_copy(R, P);
1158
1159
1160         n0 = ctx->bn[ctx->tos];
1161         n1 = ctx->bn[ctx->tos + 1];
1162         n2 = ctx->bn[ctx->tos + 2];
1163         n3 = ctx->bn[ctx->tos + 3];
1164         n4 = ctx->bn[ctx->tos + 4];
1165         n5 = ctx->bn[ctx->tos + 5];
1166         n6 = ctx->bn[ctx->tos + 6];
1167         ctx->tos += 7;
1168
1169
1170         p = E->p;
1171
1172         R->is_in_mont = 1;
1173
1174         /* L1; L2 */
1175         if (!BN_mont_mod_mul(n6, Q->Z, Q->Z, mont, ctx)) goto err;
1176         if (!BN_mont_mod_mul(n1, P->X, n6, mont, ctx)) goto err;        /* L1 = x_p * z_q^2 */
1177
1178         if (!BN_mont_mod_mul(n0, n6, Q->Z, mont, ctx)) goto err;
1179         if (!BN_mont_mod_mul(n2, P->Y, n0, mont, ctx)) goto err;        /* L2 = y_p * z_q^3 */
1180
1181
1182         /* L3; L4 */
1183         if (!BN_mont_mod_mul(n6, P->Z, P->Z, mont, ctx)) goto err;
1184         if (!BN_mont_mod_mul(n3, Q->X, n6, mont, ctx)) goto err;        /* L3 = x_q * z_p^2 */
1185
1186         if (!BN_mont_mod_mul(n0, n6, P->Z, mont, ctx)) goto err;
1187         if (!BN_mont_mod_mul(n4, Q->Y, n0, mont, ctx)) goto err;        /* L4 = y_q * z_p^3 */
1188
1189
1190         /* L5; L6 */
1191         if (!BN_mod_sub_quick(n5, n1, n3, p)) goto err;                 /* L5 = L1 - L3 */
1192         if (!BN_mod_sub_quick(n6, n2, n4, p)) goto err;                 /*L6 = L2 - L4 */
1193
1194
1195         /* pata */
1196         if (BN_is_zero(n5))
1197         {
1198                 if (BN_is_zero(n6))  /* P = Q => P + Q = 2P */
1199                 {
1200                         ctx->tos -= 7;
1201                         return ECP_mont_double(R, P, E, mont, ctx);
1202                 }
1203                 else                             /* P = -Q => P + Q = \infty */
1204                 {
1205                         ctx->tos -= 7;
1206                         if (!BN_zero(R->Z)) return 0;
1207                         return 1;
1208                 }
1209         }
1210
1211         /* L7; L8 */
1212         if (!BN_mod_add_quick(n1, n1, n3, p)) goto err;                 /* L7 = L1 + L3 */
1213         if (!BN_mod_add_quick(n2, n2, n4, p)) goto err;                 /* L8 = L2 + L4 */
1214
1215
1216         /* Z */
1217         if (!BN_mont_mod_mul(n0, P->Z, Q->Z, mont, ctx)) goto err;
1218         if (!BN_mont_mod_mul(R->Z, n0, n5, mont, ctx)) goto err;        /* Z = z_p * z_q * L_5 */
1219
1220
1221         /* X */
1222         if (!BN_mont_mod_mul(n0, n6, n6, mont, ctx)) goto err;
1223         if (!BN_mont_mod_mul(n4, n5, n5, mont, ctx)) goto err;
1224         if (!BN_mont_mod_mul(n3, n1, n4, mont, ctx)) goto err;
1225         if (!BN_mod_sub_quick(R->X, n0, n3, p)) goto err;                       /* X = L6^2 - L5^2 * L7 */
1226
1227
1228         /* L9 */
1229         if (!BN_mod_lshift1_quick(n0, R->X, p)) goto err;
1230         if (!BN_mod_sub_quick(n3, n3, n0, p)) goto err;                 /* L9 = L5^2 * L7 - 2X */
1231
1232
1233         /* Y */
1234         if (!BN_mont_mod_mul(n0, n3, n6, mont, ctx)) goto err;
1235         if (!BN_mont_mod_mul(n6, n4, n5, mont, ctx)) goto err;
1236         if (!BN_mont_mod_mul(n1, n2, n6, mont, ctx)) goto err;
1237         if (!BN_mod_sub_quick(n0, n0, n1, p)) goto err;
1238         if (!BN_mont_mod_mul(R->Y, n0, E->h, mont, ctx)) goto err;      /* Y = (L6 * L9 - L8 * L5^3) / 2 */
1239
1240
1241         ctx->tos -= 7;
1242         return 1;
1243
1244 err:
1245         ctx->tos -= 7;
1246         return 0;
1247 }
1248
1249
1250 ECP_PRECOMPUTE *ECP_mont_precompute(int r, EC_POINT *P, EC *E, BN_MONTGOMERY *mont, BN_CTX *ctx)
1251 {
1252         ECP_PRECOMPUTE *ret;
1253         EC_POINT *P2;
1254         int i, max;
1255
1256         assert(r > 2);
1257         assert(r < sizeof(unsigned int) * 8 - 1);
1258
1259         assert(mont != NULL);
1260         assert(mont->p != NULL);
1261
1262         if (!P->is_in_mont)
1263                 if (!ECP_to_montgomery(P, mont, ctx)) return 0;
1264
1265         if (!E->is_in_mont)
1266                 if (!EC_to_montgomery(E, mont, ctx)) return 0;
1267
1268         ret=(ECP_PRECOMPUTE *)malloc(sizeof(ECP_PRECOMPUTE));
1269         if (ret == NULL) return NULL;
1270
1271         max = 1;
1272         max <<= (r - 1);
1273
1274         ret->r = 0;
1275
1276         ret->Pi=(EC_POINT **)malloc(sizeof(EC_POINT *) * max);
1277         if (ret->Pi == NULL) goto err;
1278
1279
1280         /* P2 = [2]P */
1281         if ((P2 = ECP_new()) == NULL) goto err;
1282         if (!ECP_mont_double(P2, P, E, mont, ctx)) goto err;
1283
1284         /* P_0 = P */
1285         if((ret->Pi[0] = ECP_dup(P)) == NULL) goto err;
1286
1287
1288         /* P_i = P_(i-1) + P2 */
1289         for (i = 1; i < max; i++)
1290         {
1291                 if ((ret->Pi[i] = ECP_new()) == NULL) goto err;
1292                 if (!ECP_mont_add(ret->Pi[i], P2, ret->Pi[i - 1], E, mont, ctx)) goto err;
1293         }
1294
1295         ret->r = r;
1296         ECP_clear_free(P2);
1297
1298         return ret;
1299
1300 err:
1301         ECP_clear_free(P2);
1302         ECP_clear_free_precompute(ret);
1303         return NULL;
1304 }
1305
1306 int ECP_mont_multiply(EC_POINT *R, BIGNUM *k, ECP_PRECOMPUTE *prec, EC *E, BN_MONTGOMERY *mont, BN_CTX *ctx)
1307 /* R = [k]P   P = prec->Pi[0]*/
1308 {
1309         int j;
1310         int t, nextw, h, r;
1311
1312         assert(R != NULL);
1313         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
1314
1315         assert(E != NULL);
1316         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
1317
1318         assert(k != NULL);
1319         assert(!k->neg);
1320
1321         assert(ctx != NULL);
1322         assert(prec != NULL);
1323
1324         assert(mont != NULL);
1325         assert(mont->p != NULL);
1326
1327         if (!E->is_in_mont)
1328                 if (!EC_to_montgomery(E, mont, ctx)) return 0;
1329
1330
1331         if (BN_is_zero(k))
1332         {
1333                 if (!BN_zero(R->Z)) return 0;
1334                 R->is_in_mont = 1;
1335                 return 1;
1336         }
1337
1338         j = BN_num_bits(k);
1339         j--;
1340
1341         r = prec->r;
1342
1343         if (!BN_zero(R->Z)) return 0;
1344         R->is_in_mont = 1;
1345
1346         while(j >= 0)
1347         {
1348                 if (!BN_is_bit_set(k, j))
1349                 {
1350                         if (!ECP_mont_double(R, R, E, mont, ctx)) return 0;
1351                         j--;
1352                 }
1353                 else
1354                 {
1355                         nextw = j - r;
1356                         if (nextw < -1) nextw = -1;
1357                         t = nextw + 1;
1358                         while(!BN_is_bit_set(k, t))
1359                         {
1360                                 t++;
1361                         }
1362
1363                         if (!ECP_mont_double(R, R, E, mont, ctx)) return 0;
1364
1365                         j--;
1366                         if (j < t) h = 0;
1367                         else
1368                         {
1369                                 h = 1;
1370                                 for(; j > t; j--)
1371                                 {
1372                                         h <<= 1;
1373                                         if (BN_is_bit_set(k, j)) h++;
1374                                         if (!ECP_mont_double(R, R, E, mont, ctx)) return 0;
1375                                 }
1376                                 if (!ECP_mont_double(R, R, E, mont, ctx)) return 0;
1377                                 j--;
1378                         }
1379
1380                         if (!ECP_mont_add(R, R, prec->Pi[h], E, mont, ctx)) return 0;
1381
1382                         for (; j > nextw; j--)
1383                         {
1384                                 if (!ECP_mont_double(R, R, E, mont, ctx)) return 0;
1385                         }
1386
1387                 }
1388         }
1389
1390         return 1;
1391 }
1392
1393
1394 int ECP_mont_multiply2(EC_POINT *R, BIGNUM *k, EC_POINT *P, EC *E, BN_MONTGOMERY *mont, BN_CTX *ctx)
1395 /* R = [k]P */
1396 {
1397         int j, hj, kj;
1398         BIGNUM *h;
1399         EC_POINT *mP;
1400
1401         assert(R != NULL);
1402         assert(R->X != NULL && R->Y != NULL && R->Z != NULL);
1403
1404         assert(P != NULL);
1405         assert(P->X != NULL && P->Y != NULL && P->Z != NULL);
1406
1407         assert(E != NULL);
1408         assert(E->A != NULL && E->B != NULL && E->p != NULL && E->h != NULL);
1409
1410         assert(k != NULL);
1411         assert(!k->neg);
1412
1413         assert(ctx != NULL);
1414
1415         assert(mont != NULL);
1416         assert(mont->p != NULL);
1417
1418         if (!E->is_in_mont)
1419                 if (!EC_to_montgomery(E, mont, ctx)) return 0;
1420
1421         if (!P->is_in_mont)
1422                 if (!ECP_to_montgomery(P, mont, ctx)) return 0;
1423
1424
1425         if (BN_is_zero(k))
1426         {
1427                 if (!BN_zero(R->Z)) return 0;
1428                 R->is_in_mont = 1;
1429                 return 1;
1430         }
1431
1432         if ((h = BN_dup(k)) == NULL) return 0;
1433
1434         if (!BN_lshift1(h, h)) goto err;
1435         if (!BN_add(h, h, k)) goto err;
1436
1437         if (!ECP_copy(R, P)) goto err;
1438         if ((mP = ECP_mont_minus(P, mont)) == NULL) goto err;
1439
1440         for(j = BN_num_bits(h) - 2; j > 0; j--)
1441         {
1442                 if (!ECP_mont_double(R, R, E, mont, ctx)) goto err;
1443                 kj = BN_is_bit_set(k, j);
1444                 hj = BN_is_bit_set(h, j);
1445                 if (hj == 1 && kj == 0)
1446                         if (!ECP_mont_add(R, R, P, E, mont, ctx)) goto err;
1447                 if (hj == 0 && kj == 1)
1448                         if (!ECP_mont_add(R, R, mP, E, mont, ctx)) goto err;
1449         }
1450
1451         if (h != NULL) BN_free(h);
1452         if (mP != NULL) ECP_clear_free(mP);
1453         return 1;
1454
1455 err:
1456         if (h != NULL) BN_free(h);
1457         if (mP != NULL) ECP_clear_free(mP);
1458         return 0;
1459 }
1460
1461 #endif /* MONTGOMERY */