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