mpi: Add mpi_snatch and change an internal typedef.
[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
66 /* Initialize the fields of a point object.  gcry_mpi_point_free_parts
67    may be used to release the fields.  */
68 void
69 _gcry_mpi_point_init (mpi_point_t p)
70 {
71   p->x = mpi_new (0);
72   p->y = mpi_new (0);
73   p->z = mpi_new (0);
74 }
75
76
77 /* Release the parts of a point object. */
78 void
79 _gcry_mpi_point_free_parts (mpi_point_t p)
80 {
81   mpi_free (p->x); p->x = NULL;
82   mpi_free (p->y); p->y = NULL;
83   mpi_free (p->z); p->z = NULL;
84 }
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
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 \f
379 /*  RESULT = 2 * POINT  */
380 void
381 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
382 {
383 #define x3 (result->x)
384 #define y3 (result->y)
385 #define z3 (result->z)
386 #define t1 (ctx->scratch[0])
387 #define t2 (ctx->scratch[1])
388 #define t3 (ctx->scratch[2])
389 #define l1 (ctx->scratch[3])
390 #define l2 (ctx->scratch[4])
391 #define l3 (ctx->scratch[5])
392
393   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
394     {
395       /* P_y == 0 || P_z == 0 => [1:1:0] */
396       mpi_set_ui (x3, 1);
397       mpi_set_ui (y3, 1);
398       mpi_set_ui (z3, 0);
399     }
400   else
401     {
402       if (ctx->a_is_pminus3)  /* Use the faster case.  */
403         {
404           /* L1 = 3(X - Z^2)(X + Z^2) */
405           /*                          T1: used for Z^2. */
406           /*                          T2: used for the right term.  */
407           ec_powm (t1, point->z, ctx->two, ctx);
408           ec_subm (l1, point->x, t1, ctx);
409           ec_mulm (l1, l1, ctx->three, ctx);
410           ec_addm (t2, point->x, t1, ctx);
411           ec_mulm (l1, l1, t2, ctx);
412         }
413       else /* Standard case. */
414         {
415           /* L1 = 3X^2 + aZ^4 */
416           /*                          T1: used for aZ^4. */
417           ec_powm (l1, point->x, ctx->two, ctx);
418           ec_mulm (l1, l1, ctx->three, ctx);
419           ec_powm (t1, point->z, ctx->four, ctx);
420           ec_mulm (t1, t1, ctx->a, ctx);
421           ec_addm (l1, l1, t1, ctx);
422         }
423       /* Z3 = 2YZ */
424       ec_mulm (z3, point->y, point->z, ctx);
425       ec_mulm (z3, z3, ctx->two, ctx);
426
427       /* L2 = 4XY^2 */
428       /*                              T2: used for Y2; required later. */
429       ec_powm (t2, point->y, ctx->two, ctx);
430       ec_mulm (l2, t2, point->x, ctx);
431       ec_mulm (l2, l2, ctx->four, ctx);
432
433       /* X3 = L1^2 - 2L2 */
434       /*                              T1: used for L2^2. */
435       ec_powm (x3, l1, ctx->two, ctx);
436       ec_mulm (t1, l2, ctx->two, ctx);
437       ec_subm (x3, x3, t1, ctx);
438
439       /* L3 = 8Y^4 */
440       /*                              T2: taken from above. */
441       ec_powm (t2, t2, ctx->two, ctx);
442       ec_mulm (l3, t2, ctx->eight, ctx);
443
444       /* Y3 = L1(L2 - X3) - L3 */
445       ec_subm (y3, l2, x3, ctx);
446       ec_mulm (y3, y3, l1, ctx);
447       ec_subm (y3, y3, l3, ctx);
448     }
449
450 #undef x3
451 #undef y3
452 #undef z3
453 #undef t1
454 #undef t2
455 #undef t3
456 #undef l1
457 #undef l2
458 #undef l3
459 }
460
461
462
463 /* RESULT = P1 + P2 */
464 void
465 _gcry_mpi_ec_add_points (mpi_point_t result,
466                          mpi_point_t p1, mpi_point_t p2,
467                          mpi_ec_t ctx)
468 {
469 #define x1 (p1->x    )
470 #define y1 (p1->y    )
471 #define z1 (p1->z    )
472 #define x2 (p2->x    )
473 #define y2 (p2->y    )
474 #define z2 (p2->z    )
475 #define x3 (result->x)
476 #define y3 (result->y)
477 #define z3 (result->z)
478 #define l1 (ctx->scratch[0])
479 #define l2 (ctx->scratch[1])
480 #define l3 (ctx->scratch[2])
481 #define l4 (ctx->scratch[3])
482 #define l5 (ctx->scratch[4])
483 #define l6 (ctx->scratch[5])
484 #define l7 (ctx->scratch[6])
485 #define l8 (ctx->scratch[7])
486 #define l9 (ctx->scratch[8])
487 #define t1 (ctx->scratch[9])
488 #define t2 (ctx->scratch[10])
489
490   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
491     {
492       /* Same point; need to call the duplicate function.  */
493       _gcry_mpi_ec_dup_point (result, p1, ctx);
494     }
495   else if (!mpi_cmp_ui (z1, 0))
496     {
497       /* P1 is at infinity.  */
498       mpi_set (x3, p2->x);
499       mpi_set (y3, p2->y);
500       mpi_set (z3, p2->z);
501     }
502   else if (!mpi_cmp_ui (z2, 0))
503     {
504       /* P2 is at infinity.  */
505       mpi_set (x3, p1->x);
506       mpi_set (y3, p1->y);
507       mpi_set (z3, p1->z);
508     }
509   else
510     {
511       int z1_is_one = !mpi_cmp_ui (z1, 1);
512       int z2_is_one = !mpi_cmp_ui (z2, 1);
513
514       /* l1 = x1 z2^2  */
515       /* l2 = x2 z1^2  */
516       if (z2_is_one)
517         mpi_set (l1, x1);
518       else
519         {
520           ec_powm (l1, z2, ctx->two, ctx);
521           ec_mulm (l1, l1, x1, ctx);
522         }
523       if (z1_is_one)
524         mpi_set (l2, x1);
525       else
526         {
527           ec_powm (l2, z1, ctx->two, ctx);
528           ec_mulm (l2, l2, x2, ctx);
529         }
530       /* l3 = l1 - l2 */
531       ec_subm (l3, l1, l2, ctx);
532       /* l4 = y1 z2^3  */
533       ec_powm (l4, z2, ctx->three, ctx);
534       ec_mulm (l4, l4, y1, ctx);
535       /* l5 = y2 z1^3  */
536       ec_powm (l5, z1, ctx->three, ctx);
537       ec_mulm (l5, l5, y2, ctx);
538       /* l6 = l4 - l5  */
539       ec_subm (l6, l4, l5, ctx);
540
541       if (!mpi_cmp_ui (l3, 0))
542         {
543           if (!mpi_cmp_ui (l6, 0))
544             {
545               /* P1 and P2 are the same - use duplicate function.  */
546               _gcry_mpi_ec_dup_point (result, p1, ctx);
547             }
548           else
549             {
550               /* P1 is the inverse of P2.  */
551               mpi_set_ui (x3, 1);
552               mpi_set_ui (y3, 1);
553               mpi_set_ui (z3, 0);
554             }
555         }
556       else
557         {
558           /* l7 = l1 + l2  */
559           ec_addm (l7, l1, l2, ctx);
560           /* l8 = l4 + l5  */
561           ec_addm (l8, l4, l5, ctx);
562           /* z3 = z1 z2 l3  */
563           ec_mulm (z3, z1, z2, ctx);
564           ec_mulm (z3, z3, l3, ctx);
565           /* x3 = l6^2 - l7 l3^2  */
566           ec_powm (t1, l6, ctx->two, ctx);
567           ec_powm (t2, l3, ctx->two, ctx);
568           ec_mulm (t2, t2, l7, ctx);
569           ec_subm (x3, t1, t2, ctx);
570           /* l9 = l7 l3^2 - 2 x3  */
571           ec_mulm (t1, x3, ctx->two, ctx);
572           ec_subm (l9, t2, t1, ctx);
573           /* y3 = (l9 l6 - l8 l3^3)/2  */
574           ec_mulm (l9, l9, l6, ctx);
575           ec_powm (t1, l3, ctx->three, ctx); /* fixme: Use saved value*/
576           ec_mulm (t1, t1, l8, ctx);
577           ec_subm (y3, l9, t1, ctx);
578           ec_mulm (y3, y3, ctx->two_inv_p, ctx);
579         }
580     }
581
582 #undef x1
583 #undef y1
584 #undef z1
585 #undef x2
586 #undef y2
587 #undef z2
588 #undef x3
589 #undef y3
590 #undef z3
591 #undef l1
592 #undef l2
593 #undef l3
594 #undef l4
595 #undef l5
596 #undef l6
597 #undef l7
598 #undef l8
599 #undef l9
600 #undef t1
601 #undef t2
602 }
603
604
605
606 /* Scalar point multiplication - the main function for ECC.  If takes
607    an integer SCALAR and a POINT as well as the usual context CTX.
608    RESULT will be set to the resulting point. */
609 void
610 _gcry_mpi_ec_mul_point (mpi_point_t result,
611                         gcry_mpi_t scalar, mpi_point_t point,
612                         mpi_ec_t ctx)
613 {
614 #if 0
615   /* Simple left to right binary method.  GECC Algorithm 3.27 */
616   unsigned int nbits;
617   int i;
618
619   nbits = mpi_get_nbits (scalar);
620   mpi_set_ui (result->x, 1);
621   mpi_set_ui (result->y, 1);
622   mpi_set_ui (result->z, 0);
623
624   for (i=nbits-1; i >= 0; i--)
625     {
626       _gcry_mpi_ec_dup_point (result, result, ctx);
627       if (mpi_test_bit (scalar, i) == 1)
628         _gcry_mpi_ec_add_points (result, result, point, ctx);
629     }
630
631 #else
632   gcry_mpi_t x1, y1, z1, k, h, yy;
633   unsigned int i, loops;
634   mpi_point_struct p1, p2, p1inv;
635
636   x1 = mpi_alloc_like (ctx->p);
637   y1 = mpi_alloc_like (ctx->p);
638   h  = mpi_alloc_like (ctx->p);
639   k  = mpi_copy (scalar);
640   yy = mpi_copy (point->y);
641
642   if ( mpi_is_neg (k) )
643     {
644       k->sign = 0;
645       ec_invm (yy, yy, ctx);
646     }
647
648   if (!mpi_cmp_ui (point->z, 1))
649     {
650       mpi_set (x1, point->x);
651       mpi_set (y1, yy);
652     }
653   else
654     {
655       gcry_mpi_t z2, z3;
656
657       z2 = mpi_alloc_like (ctx->p);
658       z3 = mpi_alloc_like (ctx->p);
659       ec_mulm (z2, point->z, point->z, ctx);
660       ec_mulm (z3, point->z, z2, ctx);
661       ec_invm (z2, z2, ctx);
662       ec_mulm (x1, point->x, z2, ctx);
663       ec_invm (z3, z3, ctx);
664       ec_mulm (y1, yy, z3, ctx);
665       mpi_free (z2);
666       mpi_free (z3);
667     }
668   z1 = mpi_copy (ctx->one);
669
670   mpi_mul (h, k, ctx->three); /* h = 3k */
671   loops = mpi_get_nbits (h);
672
673   mpi_set (result->x, point->x);
674   mpi_set (result->y, yy); mpi_free (yy); yy = NULL;
675   mpi_set (result->z, point->z);
676
677   p1.x = x1; x1 = NULL;
678   p1.y = y1; y1 = NULL;
679   p1.z = z1; z1 = NULL;
680   point_init (&p2);
681   point_init (&p1inv);
682
683   for (i=loops-2; i > 0; i--)
684     {
685       _gcry_mpi_ec_dup_point (result, result, ctx);
686       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
687         {
688           point_set (&p2, result);
689           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
690         }
691       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
692         {
693           point_set (&p2, result);
694           /* Invert point: y = p - y mod p  */
695           point_set (&p1inv, &p1);
696           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
697           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
698         }
699     }
700
701   point_free (&p1);
702   point_free (&p2);
703   point_free (&p1inv);
704   mpi_free (h);
705   mpi_free (k);
706 #endif
707 }