mpi: Add an API for EC math.
[libgcrypt.git] / mpi / ec.c
1 /* ec.c -  Elliptic Curve functions
2  * Copyright (C) 2007 Free Software Foundation, Inc.
3  * Copyright (C) 2013 g10 Code GmbH
4  *
5  * This file is part of Libgcrypt.
6  *
7  * Libgcrypt is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU Lesser General Public License as
9  * published by the Free Software Foundation; either version 2.1 of
10  * the License, or (at your option) any later version.
11  *
12  * Libgcrypt is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with this program; if not, see <http://www.gnu.org/licenses/>.
19  */
20
21 #include <config.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <errno.h>
25
26 #include "mpi-internal.h"
27 #include "longlong.h"
28 #include "g10lib.h"
29 #include "context.h"
30
31
32 #define point_init(a)  _gcry_mpi_point_init ((a))
33 #define point_free(a)  _gcry_mpi_point_free_parts ((a))
34
35
36 /* Object to represent a point in projective coordinates. */
37 /* Currently defined in mpi.h */
38
39 /* This context is used with all our EC functions. */
40 struct mpi_ec_ctx_s
41 {
42   /* Domain parameters.  */
43   gcry_mpi_t p;   /* Prime specifying the field GF(p).  */
44   gcry_mpi_t a;   /* First coefficient of the Weierstrass equation.  */
45
46   int a_is_pminus3;  /* True if A = P - 3. */
47
48   /* Some often used constants.  */
49   gcry_mpi_t one;
50   gcry_mpi_t two;
51   gcry_mpi_t three;
52   gcry_mpi_t four;
53   gcry_mpi_t eight;
54   gcry_mpi_t two_inv_p;
55
56   /* Scratch variables.  */
57   gcry_mpi_t scratch[11];
58
59   /* Helper for fast reduction.  */
60 /*   int nist_nbits; /\* If this is a NIST curve, the number of bits.  *\/ */
61 /*   gcry_mpi_t s[10]; */
62 /*   gcry_mpi_t c; */
63
64 };
65
66
67 /* Create a new point option.  NBITS gives the size in bits of one
68    coordinate; it is only used to pre-allocate some resources and
69    might also be passed as 0 to use a default value.  */
70 mpi_point_t
71 gcry_mpi_point_new (unsigned int nbits)
72 {
73   mpi_point_t p;
74
75   (void)nbits;  /* Currently not used.  */
76
77   p = gcry_xmalloc (sizeof *p);
78   _gcry_mpi_point_init (p);
79   return p;
80 }
81
82
83 /* Release the point object P.  P may be NULL. */
84 void
85 gcry_mpi_point_release (mpi_point_t p)
86 {
87   if (p)
88     {
89       _gcry_mpi_point_free_parts (p);
90       gcry_free (p);
91     }
92 }
93
94
95 /* Initialize the fields of a point object.  gcry_mpi_point_free_parts
96    may be used to release the fields.  */
97 void
98 _gcry_mpi_point_init (mpi_point_t p)
99 {
100   p->x = mpi_new (0);
101   p->y = mpi_new (0);
102   p->z = mpi_new (0);
103 }
104
105
106 /* Release the parts of a point object. */
107 void
108 _gcry_mpi_point_free_parts (mpi_point_t p)
109 {
110   mpi_free (p->x); p->x = NULL;
111   mpi_free (p->y); p->y = NULL;
112   mpi_free (p->z); p->z = NULL;
113 }
114
115
116 /* Set the value from S into D.  */
117 static void
118 point_set (mpi_point_t d, mpi_point_t s)
119 {
120   mpi_set (d->x, s->x);
121   mpi_set (d->y, s->y);
122   mpi_set (d->z, s->z);
123 }
124
125 /* Set the projective coordinates from POINT into X, Y, and Z.  If a
126    coordinate is not required, X, Y, or Z may be passed as NULL.  */
127 void
128 gcry_mpi_point_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
129                     mpi_point_t point)
130 {
131   if (x)
132     mpi_set (x, point->x);
133   if (y)
134     mpi_set (y, point->y);
135   if (z)
136     mpi_set (z, point->z);
137 }
138
139
140 /* Set the projective coordinates from POINT into X, Y, and Z and
141    release POINT.  If a coordinate is not required, X, Y, or Z may be
142    passed as NULL.  */
143 void
144 gcry_mpi_point_snatch_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
145                            mpi_point_t point)
146 {
147   mpi_snatch (x, point->x);
148   mpi_snatch (y, point->y);
149   mpi_snatch (z, point->z);
150   gcry_free (point);
151 }
152
153
154 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
155    coordinate is given as NULL, the value 0 is stored into point.  If
156    POINT is given as NULL a new point object is allocated.  Returns
157    POINT or the newly allocated point object. */
158 mpi_point_t
159 gcry_mpi_point_set (mpi_point_t point,
160                     gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
161 {
162   if (!point)
163     point = gcry_mpi_point_new (0);
164
165   if (x)
166     mpi_set (point->x, x);
167   else
168     mpi_clear (point->x);
169   if (y)
170     mpi_set (point->y, y);
171   else
172     mpi_clear (point->y);
173   if (z)
174     mpi_set (point->z, z);
175   else
176     mpi_clear (point->z);
177
178   return point;
179 }
180
181
182 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
183    coordinate is given as NULL, the value 0 is stored into point.  If
184    POINT is given as NULL a new point object is allocated.  The
185    coordinates X, Y, and Z are released.  Returns POINT or the newly
186    allocated point object. */
187 mpi_point_t
188 gcry_mpi_point_snatch_set (mpi_point_t point,
189                            gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
190 {
191   if (!point)
192     point = gcry_mpi_point_new (0);
193
194   if (x)
195     mpi_snatch (point->x, x);
196   else
197     mpi_clear (point->x);
198   if (y)
199     mpi_snatch (point->y, y);
200   else
201     mpi_clear (point->y);
202   if (z)
203     mpi_snatch (point->z, z);
204   else
205     mpi_clear (point->z);
206
207   return point;
208 }
209
210
211 static void
212 ec_addm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
213 {
214   mpi_addm (w, u, v, ctx->p);
215 }
216
217 static void
218 ec_subm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
219 {
220   mpi_subm (w, u, v, ctx->p);
221 }
222
223 static void
224 ec_mulm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
225 {
226 #if 0
227   /* NOTE: This code works only for limb sizes of 32 bit.  */
228   mpi_limb_t *wp, *sp;
229
230   if (ctx->nist_nbits == 192)
231     {
232       mpi_mul (w, u, v);
233       mpi_resize (w, 12);
234       wp = w->d;
235
236       sp = ctx->s[0]->d;
237       sp[0*2+0] = wp[0*2+0];
238       sp[0*2+1] = wp[0*2+1];
239       sp[1*2+0] = wp[1*2+0];
240       sp[1*2+1] = wp[1*2+1];
241       sp[2*2+0] = wp[2*2+0];
242       sp[2*2+1] = wp[2*2+1];
243
244       sp = ctx->s[1]->d;
245       sp[0*2+0] = wp[3*2+0];
246       sp[0*2+1] = wp[3*2+1];
247       sp[1*2+0] = wp[3*2+0];
248       sp[1*2+1] = wp[3*2+1];
249       sp[2*2+0] = 0;
250       sp[2*2+1] = 0;
251
252       sp = ctx->s[2]->d;
253       sp[0*2+0] = 0;
254       sp[0*2+1] = 0;
255       sp[1*2+0] = wp[4*2+0];
256       sp[1*2+1] = wp[4*2+1];
257       sp[2*2+0] = wp[4*2+0];
258       sp[2*2+1] = wp[4*2+1];
259
260       sp = ctx->s[3]->d;
261       sp[0*2+0] = wp[5*2+0];
262       sp[0*2+1] = wp[5*2+1];
263       sp[1*2+0] = wp[5*2+0];
264       sp[1*2+1] = wp[5*2+1];
265       sp[2*2+0] = wp[5*2+0];
266       sp[2*2+1] = wp[5*2+1];
267
268       ctx->s[0]->nlimbs = 6;
269       ctx->s[1]->nlimbs = 6;
270       ctx->s[2]->nlimbs = 6;
271       ctx->s[3]->nlimbs = 6;
272
273       mpi_add (ctx->c, ctx->s[0], ctx->s[1]);
274       mpi_add (ctx->c, ctx->c, ctx->s[2]);
275       mpi_add (ctx->c, ctx->c, ctx->s[3]);
276
277       while ( mpi_cmp (ctx->c, ctx->p ) >= 0 )
278         mpi_sub ( ctx->c, ctx->c, ctx->p );
279       mpi_set (w, ctx->c);
280     }
281   else if (ctx->nist_nbits == 384)
282     {
283       int i;
284       mpi_mul (w, u, v);
285       mpi_resize (w, 24);
286       wp = w->d;
287
288 #define NEXT(a) do { ctx->s[(a)]->nlimbs = 12; \
289                      sp = ctx->s[(a)]->d; \
290                      i = 0; } while (0)
291 #define X(a) do { sp[i++] = wp[(a)];} while (0)
292 #define X0(a) do { sp[i++] = 0; } while (0)
293       NEXT(0);
294       X(0);X(1);X(2);X(3);X(4);X(5);X(6);X(7);X(8);X(9);X(10);X(11);
295       NEXT(1);
296       X0();X0();X0();X0();X(21);X(22);X(23);X0();X0();X0();X0();X0();
297       NEXT(2);
298       X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);X(20);X(21);X(22);X(23);
299       NEXT(3);
300       X(21);X(22);X(23);X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);X(20);
301       NEXT(4);
302       X0();X(23);X0();X(20);X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);
303       NEXT(5);
304       X0();X0();X0();X0();X(20);X(21);X(22);X(23);X0();X0();X0();X0();
305       NEXT(6);
306       X(20);X0();X0();X(21);X(22);X(23);X0();X0();X0();X0();X0();X0();
307       NEXT(7);
308       X(23);X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);X(20);X(21);X(22);
309       NEXT(8);
310       X0();X(20);X(21);X(22);X(23);X0();X0();X0();X0();X0();X0();X0();
311       NEXT(9);
312       X0();X0();X0();X(23);X(23);X0();X0();X0();X0();X0();X0();X0();
313 #undef X0
314 #undef X
315 #undef NEXT
316       mpi_add (ctx->c, ctx->s[0], ctx->s[1]);
317       mpi_add (ctx->c, ctx->c, ctx->s[1]);
318       mpi_add (ctx->c, ctx->c, ctx->s[2]);
319       mpi_add (ctx->c, ctx->c, ctx->s[3]);
320       mpi_add (ctx->c, ctx->c, ctx->s[4]);
321       mpi_add (ctx->c, ctx->c, ctx->s[5]);
322       mpi_add (ctx->c, ctx->c, ctx->s[6]);
323       mpi_sub (ctx->c, ctx->c, ctx->s[7]);
324       mpi_sub (ctx->c, ctx->c, ctx->s[8]);
325       mpi_sub (ctx->c, ctx->c, ctx->s[9]);
326
327       while ( mpi_cmp (ctx->c, ctx->p ) >= 0 )
328         mpi_sub ( ctx->c, ctx->c, ctx->p );
329       while ( ctx->c->sign )
330         mpi_add ( ctx->c, ctx->c, ctx->p );
331       mpi_set (w, ctx->c);
332     }
333   else
334 #endif /*0*/
335     mpi_mulm (w, u, v, ctx->p);
336 }
337
338 static void
339 ec_powm (gcry_mpi_t w, const gcry_mpi_t b, const gcry_mpi_t e,
340          mpi_ec_t ctx)
341 {
342   mpi_powm (w, b, e, ctx->p);
343 }
344
345 static void
346 ec_invm (gcry_mpi_t x, gcry_mpi_t a, mpi_ec_t ctx)
347 {
348   mpi_invm (x, a, ctx->p);
349 }
350
351
352
353 /* This function initialized a context for elliptic curve based on the
354    field GF(p).  P is the prime specifying this field, A is the first
355    coefficient.  CTX is expected to be zeroized.  */
356 static void
357 ec_p_init (mpi_ec_t ctx, gcry_mpi_t p, gcry_mpi_t a)
358 {
359   int i;
360   gcry_mpi_t tmp;
361
362   mpi_normalize (p);
363   mpi_normalize (a);
364
365   /* Fixme: Do we want to check some constraints? e.g.
366      a < p
367   */
368
369   ctx->p = mpi_copy (p);
370   ctx->a = mpi_copy (a);
371
372   tmp = mpi_alloc_like (ctx->p);
373   mpi_sub_ui (tmp, ctx->p, 3);
374   ctx->a_is_pminus3 = !mpi_cmp (ctx->a, tmp);
375   mpi_free (tmp);
376
377
378   /* Allocate constants.  */
379   ctx->one   = mpi_alloc_set_ui (1);
380   ctx->two   = mpi_alloc_set_ui (2);
381   ctx->three = mpi_alloc_set_ui (3);
382   ctx->four  = mpi_alloc_set_ui (4);
383   ctx->eight = mpi_alloc_set_ui (8);
384   ctx->two_inv_p = mpi_alloc (0);
385   ec_invm (ctx->two_inv_p, ctx->two, ctx);
386
387   /* Allocate scratch variables.  */
388   for (i=0; i< DIM(ctx->scratch); i++)
389     ctx->scratch[i] = mpi_alloc_like (ctx->p);
390
391   /* Prepare for fast reduction.  */
392   /* FIXME: need a test for NIST values.  However it does not gain us
393      any real advantage, for 384 bits it is actually slower than using
394      mpi_mulm.  */
395 /*   ctx->nist_nbits = mpi_get_nbits (ctx->p); */
396 /*   if (ctx->nist_nbits == 192) */
397 /*     { */
398 /*       for (i=0; i < 4; i++) */
399 /*         ctx->s[i] = mpi_new (192); */
400 /*       ctx->c    = mpi_new (192*2); */
401 /*     } */
402 /*   else if (ctx->nist_nbits == 384) */
403 /*     { */
404 /*       for (i=0; i < 10; i++) */
405 /*         ctx->s[i] = mpi_new (384); */
406 /*       ctx->c    = mpi_new (384*2); */
407 /*     } */
408 }
409
410
411 static void
412 ec_deinit (void *opaque)
413 {
414   mpi_ec_t ctx = opaque;
415   int i;
416
417   mpi_free (ctx->p);
418   mpi_free (ctx->a);
419
420   mpi_free (ctx->one);
421   mpi_free (ctx->two);
422   mpi_free (ctx->three);
423   mpi_free (ctx->four);
424   mpi_free (ctx->eight);
425
426   mpi_free (ctx->two_inv_p);
427
428   for (i=0; i< DIM(ctx->scratch); i++)
429     mpi_free (ctx->scratch[i]);
430
431 /*   if (ctx->nist_nbits == 192) */
432 /*     { */
433 /*       for (i=0; i < 4; i++) */
434 /*         mpi_free (ctx->s[i]); */
435 /*       mpi_free (ctx->c); */
436 /*     } */
437 /*   else if (ctx->nist_nbits == 384) */
438 /*     { */
439 /*       for (i=0; i < 10; i++) */
440 /*         mpi_free (ctx->s[i]); */
441 /*       mpi_free (ctx->c); */
442 /*     } */
443 }
444
445
446 /* This function returns a new context for elliptic curve based on the
447    field GF(p).  P is the prime specifying this field, A is the first
448    coefficient.  This function is only used within Libgcrypt and not
449    part of the public API.
450
451    This context needs to be released using _gcry_mpi_ec_free.  */
452 mpi_ec_t
453 _gcry_mpi_ec_p_internal_new (gcry_mpi_t p, gcry_mpi_t a)
454 {
455   mpi_ec_t ctx;
456
457   ctx = gcry_xcalloc (1, sizeof *ctx);
458   ec_p_init (ctx, p, a);
459
460   return ctx;
461 }
462
463
464 void
465 _gcry_mpi_ec_free (mpi_ec_t ctx)
466 {
467   if (ctx)
468     {
469       ec_deinit (ctx);
470       gcry_free (ctx);
471     }
472 }
473
474
475 /* This function returns a new context for elliptic curve operations
476    based on the field GF(p).  P is the prime specifying this field, A
477    is the first coefficient.  This function is part of the public API.
478    On error this function returns NULL and sets ERRNO.
479    The context needs to be released using gcry_ctx_release.  */
480 gcry_ctx_t
481 gcry_mpi_ec_p_new (gcry_mpi_t p, gcry_mpi_t a)
482 {
483   gcry_ctx_t ctx;
484   mpi_ec_t ec;
485
486   if (!p || !a || !mpi_cmp_ui (a, 0))
487     {
488       gpg_err_set_errno (EINVAL);
489       return NULL;
490     }
491
492   ctx = _gcry_ctx_alloc (CONTEXT_TYPE_EC, sizeof *ec, ec_deinit);
493   if (!ctx)
494     return NULL;
495   ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
496   ec_p_init (ec, p, a);
497
498   return ctx;
499 }
500
501
502 /* Compute the affine coordinates from the projective coordinates in
503    POINT.  Set them into X and Y.  If one coordinate is not required,
504    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
505    on success or !0 if POINT is at infinity.  */
506 int
507 _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
508                          mpi_ec_t ctx)
509 {
510   gcry_mpi_t z1, z2, z3;
511
512   if (!mpi_cmp_ui (point->z, 0))
513     return -1;
514
515   z1 = mpi_new (0);
516   z2 = mpi_new (0);
517   ec_invm (z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
518   ec_mulm (z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
519
520   if (x)
521     ec_mulm (x, point->x, z2, ctx);
522
523   if (y)
524     {
525       z3 = mpi_new (0);
526       ec_mulm (z3, z2, z1, ctx);      /* z3 = z^(-3) mod p  */
527       ec_mulm (y, point->y, z3, ctx);
528       mpi_free (z3);
529     }
530
531   mpi_free (z2);
532   mpi_free (z1);
533   return 0;
534 }
535
536
537 \f
538 /*  RESULT = 2 * POINT  */
539 void
540 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
541 {
542 #define x3 (result->x)
543 #define y3 (result->y)
544 #define z3 (result->z)
545 #define t1 (ctx->scratch[0])
546 #define t2 (ctx->scratch[1])
547 #define t3 (ctx->scratch[2])
548 #define l1 (ctx->scratch[3])
549 #define l2 (ctx->scratch[4])
550 #define l3 (ctx->scratch[5])
551
552   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
553     {
554       /* P_y == 0 || P_z == 0 => [1:1:0] */
555       mpi_set_ui (x3, 1);
556       mpi_set_ui (y3, 1);
557       mpi_set_ui (z3, 0);
558     }
559   else
560     {
561       if (ctx->a_is_pminus3)  /* Use the faster case.  */
562         {
563           /* L1 = 3(X - Z^2)(X + Z^2) */
564           /*                          T1: used for Z^2. */
565           /*                          T2: used for the right term.  */
566           ec_powm (t1, point->z, ctx->two, ctx);
567           ec_subm (l1, point->x, t1, ctx);
568           ec_mulm (l1, l1, ctx->three, ctx);
569           ec_addm (t2, point->x, t1, ctx);
570           ec_mulm (l1, l1, t2, ctx);
571         }
572       else /* Standard case. */
573         {
574           /* L1 = 3X^2 + aZ^4 */
575           /*                          T1: used for aZ^4. */
576           ec_powm (l1, point->x, ctx->two, ctx);
577           ec_mulm (l1, l1, ctx->three, ctx);
578           ec_powm (t1, point->z, ctx->four, ctx);
579           ec_mulm (t1, t1, ctx->a, ctx);
580           ec_addm (l1, l1, t1, ctx);
581         }
582       /* Z3 = 2YZ */
583       ec_mulm (z3, point->y, point->z, ctx);
584       ec_mulm (z3, z3, ctx->two, ctx);
585
586       /* L2 = 4XY^2 */
587       /*                              T2: used for Y2; required later. */
588       ec_powm (t2, point->y, ctx->two, ctx);
589       ec_mulm (l2, t2, point->x, ctx);
590       ec_mulm (l2, l2, ctx->four, ctx);
591
592       /* X3 = L1^2 - 2L2 */
593       /*                              T1: used for L2^2. */
594       ec_powm (x3, l1, ctx->two, ctx);
595       ec_mulm (t1, l2, ctx->two, ctx);
596       ec_subm (x3, x3, t1, ctx);
597
598       /* L3 = 8Y^4 */
599       /*                              T2: taken from above. */
600       ec_powm (t2, t2, ctx->two, ctx);
601       ec_mulm (l3, t2, ctx->eight, ctx);
602
603       /* Y3 = L1(L2 - X3) - L3 */
604       ec_subm (y3, l2, x3, ctx);
605       ec_mulm (y3, y3, l1, ctx);
606       ec_subm (y3, y3, l3, ctx);
607     }
608
609 #undef x3
610 #undef y3
611 #undef z3
612 #undef t1
613 #undef t2
614 #undef t3
615 #undef l1
616 #undef l2
617 #undef l3
618 }
619
620
621
622 /* RESULT = P1 + P2 */
623 void
624 _gcry_mpi_ec_add_points (mpi_point_t result,
625                          mpi_point_t p1, mpi_point_t p2,
626                          mpi_ec_t ctx)
627 {
628 #define x1 (p1->x    )
629 #define y1 (p1->y    )
630 #define z1 (p1->z    )
631 #define x2 (p2->x    )
632 #define y2 (p2->y    )
633 #define z2 (p2->z    )
634 #define x3 (result->x)
635 #define y3 (result->y)
636 #define z3 (result->z)
637 #define l1 (ctx->scratch[0])
638 #define l2 (ctx->scratch[1])
639 #define l3 (ctx->scratch[2])
640 #define l4 (ctx->scratch[3])
641 #define l5 (ctx->scratch[4])
642 #define l6 (ctx->scratch[5])
643 #define l7 (ctx->scratch[6])
644 #define l8 (ctx->scratch[7])
645 #define l9 (ctx->scratch[8])
646 #define t1 (ctx->scratch[9])
647 #define t2 (ctx->scratch[10])
648
649   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
650     {
651       /* Same point; need to call the duplicate function.  */
652       _gcry_mpi_ec_dup_point (result, p1, ctx);
653     }
654   else if (!mpi_cmp_ui (z1, 0))
655     {
656       /* P1 is at infinity.  */
657       mpi_set (x3, p2->x);
658       mpi_set (y3, p2->y);
659       mpi_set (z3, p2->z);
660     }
661   else if (!mpi_cmp_ui (z2, 0))
662     {
663       /* P2 is at infinity.  */
664       mpi_set (x3, p1->x);
665       mpi_set (y3, p1->y);
666       mpi_set (z3, p1->z);
667     }
668   else
669     {
670       int z1_is_one = !mpi_cmp_ui (z1, 1);
671       int z2_is_one = !mpi_cmp_ui (z2, 1);
672
673       /* l1 = x1 z2^2  */
674       /* l2 = x2 z1^2  */
675       if (z2_is_one)
676         mpi_set (l1, x1);
677       else
678         {
679           ec_powm (l1, z2, ctx->two, ctx);
680           ec_mulm (l1, l1, x1, ctx);
681         }
682       if (z1_is_one)
683         mpi_set (l2, x1);
684       else
685         {
686           ec_powm (l2, z1, ctx->two, ctx);
687           ec_mulm (l2, l2, x2, ctx);
688         }
689       /* l3 = l1 - l2 */
690       ec_subm (l3, l1, l2, ctx);
691       /* l4 = y1 z2^3  */
692       ec_powm (l4, z2, ctx->three, ctx);
693       ec_mulm (l4, l4, y1, ctx);
694       /* l5 = y2 z1^3  */
695       ec_powm (l5, z1, ctx->three, ctx);
696       ec_mulm (l5, l5, y2, ctx);
697       /* l6 = l4 - l5  */
698       ec_subm (l6, l4, l5, ctx);
699
700       if (!mpi_cmp_ui (l3, 0))
701         {
702           if (!mpi_cmp_ui (l6, 0))
703             {
704               /* P1 and P2 are the same - use duplicate function.  */
705               _gcry_mpi_ec_dup_point (result, p1, ctx);
706             }
707           else
708             {
709               /* P1 is the inverse of P2.  */
710               mpi_set_ui (x3, 1);
711               mpi_set_ui (y3, 1);
712               mpi_set_ui (z3, 0);
713             }
714         }
715       else
716         {
717           /* l7 = l1 + l2  */
718           ec_addm (l7, l1, l2, ctx);
719           /* l8 = l4 + l5  */
720           ec_addm (l8, l4, l5, ctx);
721           /* z3 = z1 z2 l3  */
722           ec_mulm (z3, z1, z2, ctx);
723           ec_mulm (z3, z3, l3, ctx);
724           /* x3 = l6^2 - l7 l3^2  */
725           ec_powm (t1, l6, ctx->two, ctx);
726           ec_powm (t2, l3, ctx->two, ctx);
727           ec_mulm (t2, t2, l7, ctx);
728           ec_subm (x3, t1, t2, ctx);
729           /* l9 = l7 l3^2 - 2 x3  */
730           ec_mulm (t1, x3, ctx->two, ctx);
731           ec_subm (l9, t2, t1, ctx);
732           /* y3 = (l9 l6 - l8 l3^3)/2  */
733           ec_mulm (l9, l9, l6, ctx);
734           ec_powm (t1, l3, ctx->three, ctx); /* fixme: Use saved value*/
735           ec_mulm (t1, t1, l8, ctx);
736           ec_subm (y3, l9, t1, ctx);
737           ec_mulm (y3, y3, ctx->two_inv_p, ctx);
738         }
739     }
740
741 #undef x1
742 #undef y1
743 #undef z1
744 #undef x2
745 #undef y2
746 #undef z2
747 #undef x3
748 #undef y3
749 #undef z3
750 #undef l1
751 #undef l2
752 #undef l3
753 #undef l4
754 #undef l5
755 #undef l6
756 #undef l7
757 #undef l8
758 #undef l9
759 #undef t1
760 #undef t2
761 }
762
763
764
765 /* Scalar point multiplication - the main function for ECC.  If takes
766    an integer SCALAR and a POINT as well as the usual context CTX.
767    RESULT will be set to the resulting point. */
768 void
769 _gcry_mpi_ec_mul_point (mpi_point_t result,
770                         gcry_mpi_t scalar, mpi_point_t point,
771                         mpi_ec_t ctx)
772 {
773 #if 0
774   /* Simple left to right binary method.  GECC Algorithm 3.27 */
775   unsigned int nbits;
776   int i;
777
778   nbits = mpi_get_nbits (scalar);
779   mpi_set_ui (result->x, 1);
780   mpi_set_ui (result->y, 1);
781   mpi_set_ui (result->z, 0);
782
783   for (i=nbits-1; i >= 0; i--)
784     {
785       _gcry_mpi_ec_dup_point (result, result, ctx);
786       if (mpi_test_bit (scalar, i) == 1)
787         _gcry_mpi_ec_add_points (result, result, point, ctx);
788     }
789
790 #else
791   gcry_mpi_t x1, y1, z1, k, h, yy;
792   unsigned int i, loops;
793   mpi_point_struct p1, p2, p1inv;
794
795   x1 = mpi_alloc_like (ctx->p);
796   y1 = mpi_alloc_like (ctx->p);
797   h  = mpi_alloc_like (ctx->p);
798   k  = mpi_copy (scalar);
799   yy = mpi_copy (point->y);
800
801   if ( mpi_is_neg (k) )
802     {
803       k->sign = 0;
804       ec_invm (yy, yy, ctx);
805     }
806
807   if (!mpi_cmp_ui (point->z, 1))
808     {
809       mpi_set (x1, point->x);
810       mpi_set (y1, yy);
811     }
812   else
813     {
814       gcry_mpi_t z2, z3;
815
816       z2 = mpi_alloc_like (ctx->p);
817       z3 = mpi_alloc_like (ctx->p);
818       ec_mulm (z2, point->z, point->z, ctx);
819       ec_mulm (z3, point->z, z2, ctx);
820       ec_invm (z2, z2, ctx);
821       ec_mulm (x1, point->x, z2, ctx);
822       ec_invm (z3, z3, ctx);
823       ec_mulm (y1, yy, z3, ctx);
824       mpi_free (z2);
825       mpi_free (z3);
826     }
827   z1 = mpi_copy (ctx->one);
828
829   mpi_mul (h, k, ctx->three); /* h = 3k */
830   loops = mpi_get_nbits (h);
831
832   mpi_set (result->x, point->x);
833   mpi_set (result->y, yy); mpi_free (yy); yy = NULL;
834   mpi_set (result->z, point->z);
835
836   p1.x = x1; x1 = NULL;
837   p1.y = y1; y1 = NULL;
838   p1.z = z1; z1 = NULL;
839   point_init (&p2);
840   point_init (&p1inv);
841
842   for (i=loops-2; i > 0; i--)
843     {
844       _gcry_mpi_ec_dup_point (result, result, ctx);
845       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
846         {
847           point_set (&p2, result);
848           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
849         }
850       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
851         {
852           point_set (&p2, result);
853           /* Invert point: y = p - y mod p  */
854           point_set (&p1inv, &p1);
855           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
856           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
857         }
858     }
859
860   point_free (&p1);
861   point_free (&p2);
862   point_free (&p1inv);
863   mpi_free (h);
864   mpi_free (k);
865 #endif
866 }