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