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. */
39struct 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.  */
69void
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. */
79void
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.  */
88static void
89point_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
98static void
99ec_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
104static void
105ec_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
110static void
111ec_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
225static void
226ec_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
232static void
233ec_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.  */
245mpi_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
304void
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.  */
346int
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
380/*  RESULT = 2 * POINT  */
381void
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 */
465void
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. */
610void
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