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