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