ecc: Fix memory leaks on error.
[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 #include <errno.h>
25
26 #include "mpi-internal.h"
27 #include "longlong.h"
28 #include "g10lib.h"
29 #include "context.h"
30 #include "ec-context.h"
31 #include "ec-internal.h"
32
33
34 #define point_init(a)  _gcry_mpi_point_init ((a))
35 #define point_free(a)  _gcry_mpi_point_free_parts ((a))
36
37
38 /* Print a point using the log fucntions.  If CTX is not NULL affine
39    coordinates will be printed.  */
40 void
41 _gcry_mpi_point_log (const char *name, mpi_point_t point, mpi_ec_t ctx)
42 {
43   gcry_mpi_t x, y;
44   char buf[100];
45
46   if (!point)
47     {
48       snprintf (buf, sizeof buf - 1, "%s.*", name);
49       log_mpidump (buf, NULL);
50       return;
51     }
52   snprintf (buf, sizeof buf - 1, "%s.X", name);
53
54   if (ctx)
55     {
56       x = mpi_new (0);
57       y = mpi_new (0);
58     }
59   if (!ctx || _gcry_mpi_ec_get_affine (x, y, point, ctx))
60     {
61       log_mpidump (buf, point->x);
62       buf[strlen(buf)-1] = 'Y';
63       log_mpidump (buf, point->y);
64       buf[strlen(buf)-1] = 'Z';
65       log_mpidump (buf, point->z);
66     }
67   else
68     {
69       buf[strlen(buf)-1] = 'x';
70       log_mpidump (buf, x);
71       buf[strlen(buf)-1] = 'y';
72       log_mpidump (buf, y);
73
74     }
75   if (ctx)
76     {
77       _gcry_mpi_release (x);
78       _gcry_mpi_release (y);
79     }
80 }
81
82
83 /* Create a new point option.  NBITS gives the size in bits of one
84    coordinate; it is only used to pre-allocate some resources and
85    might also be passed as 0 to use a default value.  */
86 mpi_point_t
87 _gcry_mpi_point_new (unsigned int nbits)
88 {
89   mpi_point_t p;
90
91   (void)nbits;  /* Currently not used.  */
92
93   p = xmalloc (sizeof *p);
94   _gcry_mpi_point_init (p);
95   return p;
96 }
97
98
99 /* Release the point object P.  P may be NULL. */
100 void
101 _gcry_mpi_point_release (mpi_point_t p)
102 {
103   if (p)
104     {
105       _gcry_mpi_point_free_parts (p);
106       xfree (p);
107     }
108 }
109
110
111 /* Initialize the fields of a point object.  gcry_mpi_point_free_parts
112    may be used to release the fields.  */
113 void
114 _gcry_mpi_point_init (mpi_point_t p)
115 {
116   p->x = mpi_new (0);
117   p->y = mpi_new (0);
118   p->z = mpi_new (0);
119 }
120
121
122 /* Release the parts of a point object. */
123 void
124 _gcry_mpi_point_free_parts (mpi_point_t p)
125 {
126   mpi_free (p->x); p->x = NULL;
127   mpi_free (p->y); p->y = NULL;
128   mpi_free (p->z); p->z = NULL;
129 }
130
131
132 /* Set the value from S into D.  */
133 static void
134 point_set (mpi_point_t d, mpi_point_t s)
135 {
136   mpi_set (d->x, s->x);
137   mpi_set (d->y, s->y);
138   mpi_set (d->z, s->z);
139 }
140
141
142 static void
143 point_resize (mpi_point_t p, mpi_ec_t ctx)
144 {
145   /*
146    * For now, we allocate enough limbs for our EC computation of ec_*.
147    * Once we will improve ec_* to be constant size (and constant
148    * time), NLIMBS can be ctx->p->nlimbs.
149    */
150   size_t nlimbs = 2*ctx->p->nlimbs+1;
151
152   mpi_resize (p->x, nlimbs);
153   if (ctx->model != MPI_EC_MONTGOMERY)
154     mpi_resize (p->y, nlimbs);
155   mpi_resize (p->z, nlimbs);
156 }
157
158
159 static void
160 point_swap_cond (mpi_point_t d, mpi_point_t s, unsigned long swap,
161                  mpi_ec_t ctx)
162 {
163   mpi_swap_cond (d->x, s->x, swap);
164   if (ctx->model != MPI_EC_MONTGOMERY)
165     mpi_swap_cond (d->y, s->y, swap);
166   mpi_swap_cond (d->z, s->z, swap);
167 }
168
169
170 /* Set the projective coordinates from POINT into X, Y, and Z.  If a
171    coordinate is not required, X, Y, or Z may be passed as NULL.  */
172 void
173 _gcry_mpi_point_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
174                      mpi_point_t point)
175 {
176   if (x)
177     mpi_set (x, point->x);
178   if (y)
179     mpi_set (y, point->y);
180   if (z)
181     mpi_set (z, point->z);
182 }
183
184
185 /* Set the projective coordinates from POINT into X, Y, and Z and
186    release POINT.  If a coordinate is not required, X, Y, or Z may be
187    passed as NULL.  */
188 void
189 _gcry_mpi_point_snatch_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
190                             mpi_point_t point)
191 {
192   mpi_snatch (x, point->x);
193   mpi_snatch (y, point->y);
194   mpi_snatch (z, point->z);
195   xfree (point);
196 }
197
198
199 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
200    coordinate is given as NULL, the value 0 is stored into point.  If
201    POINT is given as NULL a new point object is allocated.  Returns
202    POINT or the newly allocated point object. */
203 mpi_point_t
204 _gcry_mpi_point_set (mpi_point_t point,
205                      gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
206 {
207   if (!point)
208     point = mpi_point_new (0);
209
210   if (x)
211     mpi_set (point->x, x);
212   else
213     mpi_clear (point->x);
214   if (y)
215     mpi_set (point->y, y);
216   else
217     mpi_clear (point->y);
218   if (z)
219     mpi_set (point->z, z);
220   else
221     mpi_clear (point->z);
222
223   return point;
224 }
225
226
227 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
228    coordinate is given as NULL, the value 0 is stored into point.  If
229    POINT is given as NULL a new point object is allocated.  The
230    coordinates X, Y, and Z are released.  Returns POINT or the newly
231    allocated point object. */
232 mpi_point_t
233 _gcry_mpi_point_snatch_set (mpi_point_t point,
234                             gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
235 {
236   if (!point)
237     point = mpi_point_new (0);
238
239   if (x)
240     mpi_snatch (point->x, x);
241   else
242     mpi_clear (point->x);
243   if (y)
244     mpi_snatch (point->y, y);
245   else
246     mpi_clear (point->y);
247   if (z)
248     mpi_snatch (point->z, z);
249   else
250     mpi_clear (point->z);
251
252   return point;
253 }
254
255
256 /* W = W mod P.  */
257 static void
258 ec_mod (gcry_mpi_t w, mpi_ec_t ec)
259 {
260   if (0 && ec->dialect == ECC_DIALECT_ED25519)
261     _gcry_mpi_ec_ed25519_mod (w);
262   else if (ec->t.p_barrett)
263     _gcry_mpi_mod_barrett (w, w, ec->t.p_barrett);
264   else
265     _gcry_mpi_mod (w, w, ec->p);
266 }
267
268 static void
269 ec_addm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
270 {
271   mpi_add (w, u, v);
272   ec_mod (w, ctx);
273 }
274
275 static void
276 ec_subm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ec)
277 {
278   (void)ec;
279   mpi_sub (w, u, v);
280   /*ec_mod (w, ec);*/
281 }
282
283 static void
284 ec_mulm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
285 {
286   mpi_mul (w, u, v);
287   ec_mod (w, ctx);
288 }
289
290 /* W = 2 * U mod P.  */
291 static void
292 ec_mul2 (gcry_mpi_t w, gcry_mpi_t u, mpi_ec_t ctx)
293 {
294   mpi_lshift (w, u, 1);
295   ec_mod (w, ctx);
296 }
297
298 static void
299 ec_powm (gcry_mpi_t w, const gcry_mpi_t b, const gcry_mpi_t e,
300          mpi_ec_t ctx)
301 {
302   mpi_powm (w, b, e, ctx->p);
303   /* _gcry_mpi_abs (w); */
304 }
305
306
307 /* Shortcut for
308      ec_powm (B, B, mpi_const (MPI_C_TWO), ctx);
309    for easier optimization.  */
310 static void
311 ec_pow2 (gcry_mpi_t w, const gcry_mpi_t b, mpi_ec_t ctx)
312 {
313   /* Using mpi_mul is slightly faster (at least on amd64).  */
314   /* mpi_powm (w, b, mpi_const (MPI_C_TWO), ctx->p); */
315   ec_mulm (w, b, b, ctx);
316 }
317
318
319 /* Shortcut for
320      ec_powm (B, B, mpi_const (MPI_C_THREE), ctx);
321    for easier optimization.  */
322 static void
323 ec_pow3 (gcry_mpi_t w, const gcry_mpi_t b, mpi_ec_t ctx)
324 {
325   mpi_powm (w, b, mpi_const (MPI_C_THREE), ctx->p);
326 }
327
328
329 static void
330 ec_invm (gcry_mpi_t x, gcry_mpi_t a, mpi_ec_t ctx)
331 {
332   if (!mpi_invm (x, a, ctx->p))
333     {
334       log_error ("ec_invm: inverse does not exist:\n");
335       log_mpidump ("  a", a);
336       log_mpidump ("  p", ctx->p);
337     }
338 }
339
340
341 /* Force recomputation of all helper variables.  */
342 void
343 _gcry_mpi_ec_get_reset (mpi_ec_t ec)
344 {
345   ec->t.valid.a_is_pminus3 = 0;
346   ec->t.valid.two_inv_p = 0;
347 }
348
349
350 /* Accessor for helper variable.  */
351 static int
352 ec_get_a_is_pminus3 (mpi_ec_t ec)
353 {
354   gcry_mpi_t tmp;
355
356   if (!ec->t.valid.a_is_pminus3)
357     {
358       ec->t.valid.a_is_pminus3 = 1;
359       tmp = mpi_alloc_like (ec->p);
360       mpi_sub_ui (tmp, ec->p, 3);
361       ec->t.a_is_pminus3 = !mpi_cmp (ec->a, tmp);
362       mpi_free (tmp);
363     }
364
365   return ec->t.a_is_pminus3;
366 }
367
368
369 /* Accessor for helper variable.  */
370 static gcry_mpi_t
371 ec_get_two_inv_p (mpi_ec_t ec)
372 {
373   if (!ec->t.valid.two_inv_p)
374     {
375       ec->t.valid.two_inv_p = 1;
376       if (!ec->t.two_inv_p)
377         ec->t.two_inv_p = mpi_alloc (0);
378       ec_invm (ec->t.two_inv_p, mpi_const (MPI_C_TWO), ec);
379     }
380   return ec->t.two_inv_p;
381 }
382
383
384
385 /* This function initialized a context for elliptic curve based on the
386    field GF(p).  P is the prime specifying this field, A is the first
387    coefficient.  CTX is expected to be zeroized.  */
388 static void
389 ec_p_init (mpi_ec_t ctx, enum gcry_mpi_ec_models model,
390            enum ecc_dialects dialect,
391            int flags,
392            gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
393 {
394   int i;
395   static int use_barrett;
396
397   if (!use_barrett)
398     {
399       if (getenv ("GCRYPT_BARRETT"))
400         use_barrett = 1;
401       else
402         use_barrett = -1;
403     }
404
405   /* Fixme: Do we want to check some constraints? e.g.  a < p  */
406
407   ctx->model = model;
408   ctx->dialect = dialect;
409   ctx->flags = flags;
410   if (dialect == ECC_DIALECT_ED25519)
411     ctx->nbits = 256;
412   else
413     ctx->nbits = mpi_get_nbits (p);
414   ctx->p = mpi_copy (p);
415   ctx->a = mpi_copy (a);
416   ctx->b = mpi_copy (b);
417
418   ctx->t.p_barrett = use_barrett > 0? _gcry_mpi_barrett_init (ctx->p, 0):NULL;
419
420   _gcry_mpi_ec_get_reset (ctx);
421
422   /* Allocate scratch variables.  */
423   for (i=0; i< DIM(ctx->t.scratch); i++)
424     ctx->t.scratch[i] = mpi_alloc_like (ctx->p);
425
426   /* Prepare for fast reduction.  */
427   /* FIXME: need a test for NIST values.  However it does not gain us
428      any real advantage, for 384 bits it is actually slower than using
429      mpi_mulm.  */
430 /*   ctx->nist_nbits = mpi_get_nbits (ctx->p); */
431 /*   if (ctx->nist_nbits == 192) */
432 /*     { */
433 /*       for (i=0; i < 4; i++) */
434 /*         ctx->s[i] = mpi_new (192); */
435 /*       ctx->c    = mpi_new (192*2); */
436 /*     } */
437 /*   else if (ctx->nist_nbits == 384) */
438 /*     { */
439 /*       for (i=0; i < 10; i++) */
440 /*         ctx->s[i] = mpi_new (384); */
441 /*       ctx->c    = mpi_new (384*2); */
442 /*     } */
443 }
444
445
446 static void
447 ec_deinit (void *opaque)
448 {
449   mpi_ec_t ctx = opaque;
450   int i;
451
452   _gcry_mpi_barrett_free (ctx->t.p_barrett);
453
454   /* Domain parameter.  */
455   mpi_free (ctx->p);
456   mpi_free (ctx->a);
457   mpi_free (ctx->b);
458   _gcry_mpi_point_release (ctx->G);
459   mpi_free (ctx->n);
460   mpi_free (ctx->h);
461
462   /* The key.  */
463   _gcry_mpi_point_release (ctx->Q);
464   mpi_free (ctx->d);
465
466   /* Private data of ec.c.  */
467   mpi_free (ctx->t.two_inv_p);
468
469   for (i=0; i< DIM(ctx->t.scratch); i++)
470     mpi_free (ctx->t.scratch[i]);
471
472 /*   if (ctx->nist_nbits == 192) */
473 /*     { */
474 /*       for (i=0; i < 4; i++) */
475 /*         mpi_free (ctx->s[i]); */
476 /*       mpi_free (ctx->c); */
477 /*     } */
478 /*   else if (ctx->nist_nbits == 384) */
479 /*     { */
480 /*       for (i=0; i < 10; i++) */
481 /*         mpi_free (ctx->s[i]); */
482 /*       mpi_free (ctx->c); */
483 /*     } */
484 }
485
486
487 /* This function returns a new context for elliptic curve based on the
488    field GF(p).  P is the prime specifying this field, A is the first
489    coefficient, B is the second coefficient, and MODEL is the model
490    for the curve.  This function is only used within Libgcrypt and not
491    part of the public API.
492
493    This context needs to be released using _gcry_mpi_ec_free.  */
494 mpi_ec_t
495 _gcry_mpi_ec_p_internal_new (enum gcry_mpi_ec_models model,
496                              enum ecc_dialects dialect,
497                              int flags,
498                              gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
499 {
500   mpi_ec_t ctx;
501
502   ctx = xcalloc (1, sizeof *ctx);
503   ec_p_init (ctx, model, dialect, flags, p, a, b);
504
505   return ctx;
506 }
507
508
509 /* This is a variant of _gcry_mpi_ec_p_internal_new which returns an
510    public context and does some error checking on the supplied
511    arguments.  On success the new context is stored at R_CTX and 0 is
512    returned; on error NULL is stored at R_CTX and an error code is
513    returned.
514
515    The context needs to be released using gcry_ctx_release.  */
516 gpg_err_code_t
517 _gcry_mpi_ec_p_new (gcry_ctx_t *r_ctx,
518                     enum gcry_mpi_ec_models model,
519                     enum ecc_dialects dialect,
520                     int flags,
521                     gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
522 {
523   gcry_ctx_t ctx;
524   mpi_ec_t ec;
525
526   *r_ctx = NULL;
527   if (!p || !a)
528     return GPG_ERR_EINVAL;
529
530   ctx = _gcry_ctx_alloc (CONTEXT_TYPE_EC, sizeof *ec, ec_deinit);
531   if (!ctx)
532     return gpg_err_code_from_syserror ();
533   ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
534   ec_p_init (ec, model, dialect, flags, p, a, b);
535
536   *r_ctx = ctx;
537   return 0;
538 }
539
540
541 void
542 _gcry_mpi_ec_free (mpi_ec_t ctx)
543 {
544   if (ctx)
545     {
546       ec_deinit (ctx);
547       xfree (ctx);
548     }
549 }
550
551
552 gcry_mpi_t
553 _gcry_mpi_ec_get_mpi (const char *name, gcry_ctx_t ctx, int copy)
554 {
555   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
556
557   return _gcry_ecc_get_mpi (name, ec, copy);
558 }
559
560
561 gcry_mpi_point_t
562 _gcry_mpi_ec_get_point (const char *name, gcry_ctx_t ctx, int copy)
563 {
564   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
565
566   (void)copy;  /* Not used.  */
567
568   return _gcry_ecc_get_point (name, ec);
569 }
570
571
572 gpg_err_code_t
573 _gcry_mpi_ec_set_mpi (const char *name, gcry_mpi_t newvalue,
574                       gcry_ctx_t ctx)
575 {
576   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
577
578   return _gcry_ecc_set_mpi (name, newvalue, ec);
579 }
580
581
582 gpg_err_code_t
583 _gcry_mpi_ec_set_point (const char *name, gcry_mpi_point_t newvalue,
584                         gcry_ctx_t ctx)
585 {
586   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
587
588   return _gcry_ecc_set_point (name, newvalue, ec);
589 }
590
591
592 /* Given an encoded point in the MPI VALUE and a context EC, decode
593  * the point according to the context and store it in RESULT.  On
594  * error an error code is return but RESULT might have been changed.
595  * If no context is given the function tries to decode VALUE by
596  * assuming a 0x04 prefixed uncompressed encoding.  */
597 gpg_err_code_t
598 _gcry_mpi_ec_decode_point (mpi_point_t result, gcry_mpi_t value, mpi_ec_t ec)
599 {
600   gcry_err_code_t rc;
601
602   if (ec && ec->dialect == ECC_DIALECT_ED25519)
603     rc = _gcry_ecc_eddsa_decodepoint (value, ec, result, NULL, NULL);
604   else if (ec && ec->model == MPI_EC_MONTGOMERY)
605     rc = _gcry_ecc_mont_decodepoint (value, ec, result);
606   else
607     rc = _gcry_ecc_os2ec (result, value);
608
609   return rc;
610 }
611
612
613 /* Compute the affine coordinates from the projective coordinates in
614    POINT.  Set them into X and Y.  If one coordinate is not required,
615    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
616    on success or !0 if POINT is at infinity.  */
617 int
618 _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
619                          mpi_ec_t ctx)
620 {
621   if (!mpi_cmp_ui (point->z, 0))
622     return -1;
623
624   switch (ctx->model)
625     {
626     case MPI_EC_WEIERSTRASS: /* Using Jacobian coordinates.  */
627       {
628         gcry_mpi_t z1, z2, z3;
629
630         z1 = mpi_new (0);
631         z2 = mpi_new (0);
632         ec_invm (z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
633         ec_mulm (z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
634
635         if (x)
636           ec_mulm (x, point->x, z2, ctx);
637
638         if (y)
639           {
640             z3 = mpi_new (0);
641             ec_mulm (z3, z2, z1, ctx);      /* z3 = z^(-3) mod p  */
642             ec_mulm (y, point->y, z3, ctx);
643             mpi_free (z3);
644           }
645
646         mpi_free (z2);
647         mpi_free (z1);
648       }
649       return 0;
650
651     case MPI_EC_MONTGOMERY:
652       {
653         if (x)
654           mpi_set (x, point->x);
655
656         if (y)
657           {
658             log_fatal ("%s: Getting Y-coordinate on %s is not supported\n",
659                        "_gcry_mpi_ec_get_affine", "Montgomery");
660             return -1;
661           }
662       }
663       return 0;
664
665     case MPI_EC_EDWARDS:
666       {
667         gcry_mpi_t z;
668
669         z = mpi_new (0);
670         ec_invm (z, point->z, ctx);
671
672         if (x)
673           ec_mulm (x, point->x, z, ctx);
674         if (y)
675           ec_mulm (y, point->y, z, ctx);
676
677         _gcry_mpi_release (z);
678       }
679       return 0;
680
681     default:
682       return -1;
683     }
684 }
685
686
687 \f
688 /*  RESULT = 2 * POINT  (Weierstrass version). */
689 static void
690 dup_point_weierstrass (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
691 {
692 #define x3 (result->x)
693 #define y3 (result->y)
694 #define z3 (result->z)
695 #define t1 (ctx->t.scratch[0])
696 #define t2 (ctx->t.scratch[1])
697 #define t3 (ctx->t.scratch[2])
698 #define l1 (ctx->t.scratch[3])
699 #define l2 (ctx->t.scratch[4])
700 #define l3 (ctx->t.scratch[5])
701
702   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
703     {
704       /* P_y == 0 || P_z == 0 => [1:1:0] */
705       mpi_set_ui (x3, 1);
706       mpi_set_ui (y3, 1);
707       mpi_set_ui (z3, 0);
708     }
709   else
710     {
711       if (ec_get_a_is_pminus3 (ctx))  /* Use the faster case.  */
712         {
713           /* L1 = 3(X - Z^2)(X + Z^2) */
714           /*                          T1: used for Z^2. */
715           /*                          T2: used for the right term.  */
716           ec_pow2 (t1, point->z, ctx);
717           ec_subm (l1, point->x, t1, ctx);
718           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
719           ec_addm (t2, point->x, t1, ctx);
720           ec_mulm (l1, l1, t2, ctx);
721         }
722       else /* Standard case. */
723         {
724           /* L1 = 3X^2 + aZ^4 */
725           /*                          T1: used for aZ^4. */
726           ec_pow2 (l1, point->x, ctx);
727           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
728           ec_powm (t1, point->z, mpi_const (MPI_C_FOUR), ctx);
729           ec_mulm (t1, t1, ctx->a, ctx);
730           ec_addm (l1, l1, t1, ctx);
731         }
732       /* Z3 = 2YZ */
733       ec_mulm (z3, point->y, point->z, ctx);
734       ec_mul2 (z3, z3, ctx);
735
736       /* L2 = 4XY^2 */
737       /*                              T2: used for Y2; required later. */
738       ec_pow2 (t2, point->y, ctx);
739       ec_mulm (l2, t2, point->x, ctx);
740       ec_mulm (l2, l2, mpi_const (MPI_C_FOUR), ctx);
741
742       /* X3 = L1^2 - 2L2 */
743       /*                              T1: used for L2^2. */
744       ec_pow2 (x3, l1, ctx);
745       ec_mul2 (t1, l2, ctx);
746       ec_subm (x3, x3, t1, ctx);
747
748       /* L3 = 8Y^4 */
749       /*                              T2: taken from above. */
750       ec_pow2 (t2, t2, ctx);
751       ec_mulm (l3, t2, mpi_const (MPI_C_EIGHT), ctx);
752
753       /* Y3 = L1(L2 - X3) - L3 */
754       ec_subm (y3, l2, x3, ctx);
755       ec_mulm (y3, y3, l1, ctx);
756       ec_subm (y3, y3, l3, ctx);
757     }
758
759 #undef x3
760 #undef y3
761 #undef z3
762 #undef t1
763 #undef t2
764 #undef t3
765 #undef l1
766 #undef l2
767 #undef l3
768 }
769
770
771 /*  RESULT = 2 * POINT  (Montgomery version). */
772 static void
773 dup_point_montgomery (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
774 {
775   (void)result;
776   (void)point;
777   (void)ctx;
778   log_fatal ("%s: %s not yet supported\n",
779              "_gcry_mpi_ec_dup_point", "Montgomery");
780 }
781
782
783 /*  RESULT = 2 * POINT  (Twisted Edwards version). */
784 static void
785 dup_point_edwards (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
786 {
787 #define X1 (point->x)
788 #define Y1 (point->y)
789 #define Z1 (point->z)
790 #define X3 (result->x)
791 #define Y3 (result->y)
792 #define Z3 (result->z)
793 #define B (ctx->t.scratch[0])
794 #define C (ctx->t.scratch[1])
795 #define D (ctx->t.scratch[2])
796 #define E (ctx->t.scratch[3])
797 #define F (ctx->t.scratch[4])
798 #define H (ctx->t.scratch[5])
799 #define J (ctx->t.scratch[6])
800
801   /* Compute: (X_3 : Y_3 : Z_3) = 2( X_1 : Y_1 : Z_1 ) */
802
803   /* B = (X_1 + Y_1)^2  */
804   ec_addm (B, X1, Y1, ctx);
805   ec_pow2 (B, B, ctx);
806
807   /* C = X_1^2 */
808   /* D = Y_1^2 */
809   ec_pow2 (C, X1, ctx);
810   ec_pow2 (D, Y1, ctx);
811
812   /* E = aC */
813   if (ctx->dialect == ECC_DIALECT_ED25519)
814     {
815       mpi_set (E, C);
816       _gcry_mpi_neg (E, E);
817     }
818   else
819     ec_mulm (E, ctx->a, C, ctx);
820
821   /* F = E + D */
822   ec_addm (F, E, D, ctx);
823
824   /* H = Z_1^2 */
825   ec_pow2 (H, Z1, ctx);
826
827   /* J = F - 2H */
828   ec_mul2 (J, H, ctx);
829   ec_subm (J, F, J, ctx);
830
831   /* X_3 = (B - C - D) · J */
832   ec_subm (X3, B, C, ctx);
833   ec_subm (X3, X3, D, ctx);
834   ec_mulm (X3, X3, J, ctx);
835
836   /* Y_3 = F · (E - D) */
837   ec_subm (Y3, E, D, ctx);
838   ec_mulm (Y3, Y3, F, ctx);
839
840   /* Z_3 = F · J */
841   ec_mulm (Z3, F, J, ctx);
842
843 #undef X1
844 #undef Y1
845 #undef Z1
846 #undef X3
847 #undef Y3
848 #undef Z3
849 #undef B
850 #undef C
851 #undef D
852 #undef E
853 #undef F
854 #undef H
855 #undef J
856 }
857
858
859 /*  RESULT = 2 * POINT  */
860 void
861 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
862 {
863   switch (ctx->model)
864     {
865     case MPI_EC_WEIERSTRASS:
866       dup_point_weierstrass (result, point, ctx);
867       break;
868     case MPI_EC_MONTGOMERY:
869       dup_point_montgomery (result, point, ctx);
870       break;
871     case MPI_EC_EDWARDS:
872       dup_point_edwards (result, point, ctx);
873       break;
874     }
875 }
876
877
878 /* RESULT = P1 + P2  (Weierstrass version).*/
879 static void
880 add_points_weierstrass (mpi_point_t result,
881                         mpi_point_t p1, mpi_point_t p2,
882                         mpi_ec_t ctx)
883 {
884 #define x1 (p1->x    )
885 #define y1 (p1->y    )
886 #define z1 (p1->z    )
887 #define x2 (p2->x    )
888 #define y2 (p2->y    )
889 #define z2 (p2->z    )
890 #define x3 (result->x)
891 #define y3 (result->y)
892 #define z3 (result->z)
893 #define l1 (ctx->t.scratch[0])
894 #define l2 (ctx->t.scratch[1])
895 #define l3 (ctx->t.scratch[2])
896 #define l4 (ctx->t.scratch[3])
897 #define l5 (ctx->t.scratch[4])
898 #define l6 (ctx->t.scratch[5])
899 #define l7 (ctx->t.scratch[6])
900 #define l8 (ctx->t.scratch[7])
901 #define l9 (ctx->t.scratch[8])
902 #define t1 (ctx->t.scratch[9])
903 #define t2 (ctx->t.scratch[10])
904
905   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
906     {
907       /* Same point; need to call the duplicate function.  */
908       _gcry_mpi_ec_dup_point (result, p1, ctx);
909     }
910   else if (!mpi_cmp_ui (z1, 0))
911     {
912       /* P1 is at infinity.  */
913       mpi_set (x3, p2->x);
914       mpi_set (y3, p2->y);
915       mpi_set (z3, p2->z);
916     }
917   else if (!mpi_cmp_ui (z2, 0))
918     {
919       /* P2 is at infinity.  */
920       mpi_set (x3, p1->x);
921       mpi_set (y3, p1->y);
922       mpi_set (z3, p1->z);
923     }
924   else
925     {
926       int z1_is_one = !mpi_cmp_ui (z1, 1);
927       int z2_is_one = !mpi_cmp_ui (z2, 1);
928
929       /* l1 = x1 z2^2  */
930       /* l2 = x2 z1^2  */
931       if (z2_is_one)
932         mpi_set (l1, x1);
933       else
934         {
935           ec_pow2 (l1, z2, ctx);
936           ec_mulm (l1, l1, x1, ctx);
937         }
938       if (z1_is_one)
939         mpi_set (l2, x2);
940       else
941         {
942           ec_pow2 (l2, z1, ctx);
943           ec_mulm (l2, l2, x2, ctx);
944         }
945       /* l3 = l1 - l2 */
946       ec_subm (l3, l1, l2, ctx);
947       /* l4 = y1 z2^3  */
948       ec_powm (l4, z2, mpi_const (MPI_C_THREE), ctx);
949       ec_mulm (l4, l4, y1, ctx);
950       /* l5 = y2 z1^3  */
951       ec_powm (l5, z1, mpi_const (MPI_C_THREE), ctx);
952       ec_mulm (l5, l5, y2, ctx);
953       /* l6 = l4 - l5  */
954       ec_subm (l6, l4, l5, ctx);
955
956       if (!mpi_cmp_ui (l3, 0))
957         {
958           if (!mpi_cmp_ui (l6, 0))
959             {
960               /* P1 and P2 are the same - use duplicate function.  */
961               _gcry_mpi_ec_dup_point (result, p1, ctx);
962             }
963           else
964             {
965               /* P1 is the inverse of P2.  */
966               mpi_set_ui (x3, 1);
967               mpi_set_ui (y3, 1);
968               mpi_set_ui (z3, 0);
969             }
970         }
971       else
972         {
973           /* l7 = l1 + l2  */
974           ec_addm (l7, l1, l2, ctx);
975           /* l8 = l4 + l5  */
976           ec_addm (l8, l4, l5, ctx);
977           /* z3 = z1 z2 l3  */
978           ec_mulm (z3, z1, z2, ctx);
979           ec_mulm (z3, z3, l3, ctx);
980           /* x3 = l6^2 - l7 l3^2  */
981           ec_pow2 (t1, l6, ctx);
982           ec_pow2 (t2, l3, ctx);
983           ec_mulm (t2, t2, l7, ctx);
984           ec_subm (x3, t1, t2, ctx);
985           /* l9 = l7 l3^2 - 2 x3  */
986           ec_mul2 (t1, x3, ctx);
987           ec_subm (l9, t2, t1, ctx);
988           /* y3 = (l9 l6 - l8 l3^3)/2  */
989           ec_mulm (l9, l9, l6, ctx);
990           ec_powm (t1, l3, mpi_const (MPI_C_THREE), ctx); /* fixme: Use saved value*/
991           ec_mulm (t1, t1, l8, ctx);
992           ec_subm (y3, l9, t1, ctx);
993           ec_mulm (y3, y3, ec_get_two_inv_p (ctx), ctx);
994         }
995     }
996
997 #undef x1
998 #undef y1
999 #undef z1
1000 #undef x2
1001 #undef y2
1002 #undef z2
1003 #undef x3
1004 #undef y3
1005 #undef z3
1006 #undef l1
1007 #undef l2
1008 #undef l3
1009 #undef l4
1010 #undef l5
1011 #undef l6
1012 #undef l7
1013 #undef l8
1014 #undef l9
1015 #undef t1
1016 #undef t2
1017 }
1018
1019
1020 /* RESULT = P1 + P2  (Montgomery version).*/
1021 static void
1022 add_points_montgomery (mpi_point_t result,
1023                        mpi_point_t p1, mpi_point_t p2,
1024                        mpi_ec_t ctx)
1025 {
1026   (void)result;
1027   (void)p1;
1028   (void)p2;
1029   (void)ctx;
1030   log_fatal ("%s: %s not yet supported\n",
1031              "_gcry_mpi_ec_add_points", "Montgomery");
1032 }
1033
1034
1035 /* RESULT = P1 + P2  (Twisted Edwards version).*/
1036 static void
1037 add_points_edwards (mpi_point_t result,
1038                     mpi_point_t p1, mpi_point_t p2,
1039                     mpi_ec_t ctx)
1040 {
1041 #define X1 (p1->x)
1042 #define Y1 (p1->y)
1043 #define Z1 (p1->z)
1044 #define X2 (p2->x)
1045 #define Y2 (p2->y)
1046 #define Z2 (p2->z)
1047 #define X3 (result->x)
1048 #define Y3 (result->y)
1049 #define Z3 (result->z)
1050 #define A (ctx->t.scratch[0])
1051 #define B (ctx->t.scratch[1])
1052 #define C (ctx->t.scratch[2])
1053 #define D (ctx->t.scratch[3])
1054 #define E (ctx->t.scratch[4])
1055 #define F (ctx->t.scratch[5])
1056 #define G (ctx->t.scratch[6])
1057 #define tmp (ctx->t.scratch[7])
1058
1059   /* Compute: (X_3 : Y_3 : Z_3) = (X_1 : Y_1 : Z_1) + (X_2 : Y_2 : Z_3)  */
1060
1061   /* A = Z1 · Z2 */
1062   ec_mulm (A, Z1, Z2, ctx);
1063
1064   /* B = A^2 */
1065   ec_pow2 (B, A, ctx);
1066
1067   /* C = X1 · X2 */
1068   ec_mulm (C, X1, X2, ctx);
1069
1070   /* D = Y1 · Y2 */
1071   ec_mulm (D, Y1, Y2, ctx);
1072
1073   /* E = d · C · D */
1074   ec_mulm (E, ctx->b, C, ctx);
1075   ec_mulm (E, E, D, ctx);
1076
1077   /* F = B - E */
1078   ec_subm (F, B, E, ctx);
1079
1080   /* G = B + E */
1081   ec_addm (G, B, E, ctx);
1082
1083   /* X_3 = A · F · ((X_1 + Y_1) · (X_2 + Y_2) - C - D) */
1084   ec_addm (tmp, X1, Y1, ctx);
1085   ec_addm (X3, X2, Y2, ctx);
1086   ec_mulm (X3, X3, tmp, ctx);
1087   ec_subm (X3, X3, C, ctx);
1088   ec_subm (X3, X3, D, ctx);
1089   ec_mulm (X3, X3, F, ctx);
1090   ec_mulm (X3, X3, A, ctx);
1091
1092   /* Y_3 = A · G · (D - aC) */
1093   if (ctx->dialect == ECC_DIALECT_ED25519)
1094     {
1095       /* Using ec_addm (Y3, D, C, ctx) is possible but a litte bit
1096          slower because a subm does currently skip the mod step.  */
1097       mpi_set (Y3, C);
1098       _gcry_mpi_neg (Y3, Y3);
1099       ec_subm (Y3, D, Y3, ctx);
1100     }
1101   else
1102     {
1103       ec_mulm (Y3, ctx->a, C, ctx);
1104       ec_subm (Y3, D, Y3, ctx);
1105     }
1106   ec_mulm (Y3, Y3, G, ctx);
1107   ec_mulm (Y3, Y3, A, ctx);
1108
1109   /* Z_3 = F · G */
1110   ec_mulm (Z3, F, G, ctx);
1111
1112
1113 #undef X1
1114 #undef Y1
1115 #undef Z1
1116 #undef X2
1117 #undef Y2
1118 #undef Z2
1119 #undef X3
1120 #undef Y3
1121 #undef Z3
1122 #undef A
1123 #undef B
1124 #undef C
1125 #undef D
1126 #undef E
1127 #undef F
1128 #undef G
1129 #undef tmp
1130 }
1131
1132
1133 /* Compute a step of Montgomery Ladder (only use X and Z in the point).
1134    Inputs:  P1, P2, and x-coordinate of DIF = P1 - P1.
1135    Outputs: PRD = 2 * P1 and  SUM = P1 + P2. */
1136 static void
1137 montgomery_ladder (mpi_point_t prd, mpi_point_t sum,
1138                    mpi_point_t p1, mpi_point_t p2, gcry_mpi_t dif_x,
1139                    mpi_ec_t ctx)
1140 {
1141   ec_addm (sum->x, p2->x, p2->z, ctx);
1142   ec_subm (p2->z, p2->x, p2->z, ctx);
1143   ec_addm (prd->x, p1->x, p1->z, ctx);
1144   ec_subm (p1->z, p1->x, p1->z, ctx);
1145   ec_mulm (p2->x, p1->z, sum->x, ctx);
1146   ec_mulm (p2->z, prd->x, p2->z, ctx);
1147   ec_pow2 (p1->x, prd->x, ctx);
1148   ec_pow2 (p1->z, p1->z, ctx);
1149   ec_addm (sum->x, p2->x, p2->z, ctx);
1150   ec_subm (p2->z, p2->x, p2->z, ctx);
1151   ec_mulm (prd->x, p1->x, p1->z, ctx);
1152   ec_subm (p1->z, p1->x, p1->z, ctx);
1153   ec_pow2 (sum->x, sum->x, ctx);
1154   ec_pow2 (sum->z, p2->z, ctx);
1155   ec_mulm (prd->z, p1->z, ctx->a, ctx); /* CTX->A: (a-2)/4 */
1156   ec_mulm (sum->z, sum->z, dif_x, ctx);
1157   ec_addm (prd->z, p1->x, prd->z, ctx);
1158   ec_mulm (prd->z, prd->z, p1->z, ctx);
1159 }
1160
1161
1162 /* RESULT = P1 + P2 */
1163 void
1164 _gcry_mpi_ec_add_points (mpi_point_t result,
1165                          mpi_point_t p1, mpi_point_t p2,
1166                          mpi_ec_t ctx)
1167 {
1168   switch (ctx->model)
1169     {
1170     case MPI_EC_WEIERSTRASS:
1171       add_points_weierstrass (result, p1, p2, ctx);
1172       break;
1173     case MPI_EC_MONTGOMERY:
1174       add_points_montgomery (result, p1, p2, ctx);
1175       break;
1176     case MPI_EC_EDWARDS:
1177       add_points_edwards (result, p1, p2, ctx);
1178       break;
1179     }
1180 }
1181
1182
1183 /* RESULT = P1 - P2  (Weierstrass version).*/
1184 static void
1185 sub_points_weierstrass (mpi_point_t result,
1186                         mpi_point_t p1, mpi_point_t p2,
1187                         mpi_ec_t ctx)
1188 {
1189   (void)result;
1190   (void)p1;
1191   (void)p2;
1192   (void)ctx;
1193   log_fatal ("%s: %s not yet supported\n",
1194              "_gcry_mpi_ec_sub_points", "Weierstrass");
1195 }
1196
1197
1198 /* RESULT = P1 - P2  (Montgomery version).*/
1199 static void
1200 sub_points_montgomery (mpi_point_t result,
1201                        mpi_point_t p1, mpi_point_t p2,
1202                        mpi_ec_t ctx)
1203 {
1204   (void)result;
1205   (void)p1;
1206   (void)p2;
1207   (void)ctx;
1208   log_fatal ("%s: %s not yet supported\n",
1209              "_gcry_mpi_ec_sub_points", "Montgomery");
1210 }
1211
1212
1213 /* RESULT = P1 - P2  (Twisted Edwards version).*/
1214 static void
1215 sub_points_edwards (mpi_point_t result,
1216                     mpi_point_t p1, mpi_point_t p2,
1217                     mpi_ec_t ctx)
1218 {
1219   mpi_point_t p2i = _gcry_mpi_point_new (0);
1220   point_set (p2i, p2);
1221   _gcry_mpi_neg (p2i->x, p2i->x);
1222   add_points_edwards (result, p1, p2i, ctx);
1223   _gcry_mpi_point_release (p2i);
1224 }
1225
1226
1227 /* RESULT = P1 - P2 */
1228 void
1229 _gcry_mpi_ec_sub_points (mpi_point_t result,
1230                          mpi_point_t p1, mpi_point_t p2,
1231                          mpi_ec_t ctx)
1232 {
1233   switch (ctx->model)
1234     {
1235     case MPI_EC_WEIERSTRASS:
1236       sub_points_weierstrass (result, p1, p2, ctx);
1237       break;
1238     case MPI_EC_MONTGOMERY:
1239       sub_points_montgomery (result, p1, p2, ctx);
1240       break;
1241     case MPI_EC_EDWARDS:
1242       sub_points_edwards (result, p1, p2, ctx);
1243       break;
1244     }
1245 }
1246
1247
1248 /* Scalar point multiplication - the main function for ECC.  If takes
1249    an integer SCALAR and a POINT as well as the usual context CTX.
1250    RESULT will be set to the resulting point. */
1251 void
1252 _gcry_mpi_ec_mul_point (mpi_point_t result,
1253                         gcry_mpi_t scalar, mpi_point_t point,
1254                         mpi_ec_t ctx)
1255 {
1256   gcry_mpi_t x1, y1, z1, k, h, yy;
1257   unsigned int i, loops;
1258   mpi_point_struct p1, p2, p1inv;
1259
1260   if (ctx->model == MPI_EC_EDWARDS
1261       || (ctx->model == MPI_EC_WEIERSTRASS
1262           && mpi_is_secure (scalar)))
1263     {
1264       /* Simple left to right binary method.  GECC Algorithm 3.27 */
1265       unsigned int nbits;
1266       int j;
1267
1268       nbits = mpi_get_nbits (scalar);
1269       if (ctx->model == MPI_EC_WEIERSTRASS)
1270         {
1271           mpi_set_ui (result->x, 1);
1272           mpi_set_ui (result->y, 1);
1273           mpi_set_ui (result->z, 0);
1274         }
1275       else
1276         {
1277           mpi_set_ui (result->x, 0);
1278           mpi_set_ui (result->y, 1);
1279           mpi_set_ui (result->z, 1);
1280         }
1281
1282       if (mpi_is_secure (scalar))
1283         {
1284           /* If SCALAR is in secure memory we assume that it is the
1285              secret key we use constant time operation.  */
1286           mpi_point_struct tmppnt;
1287
1288           point_init (&tmppnt);
1289           point_resize (result, ctx);
1290           point_resize (&tmppnt, ctx);
1291           for (j=nbits-1; j >= 0; j--)
1292             {
1293               _gcry_mpi_ec_dup_point (result, result, ctx);
1294               _gcry_mpi_ec_add_points (&tmppnt, result, point, ctx);
1295               point_swap_cond (result, &tmppnt, mpi_test_bit (scalar, j), ctx);
1296             }
1297           point_free (&tmppnt);
1298         }
1299       else
1300         {
1301           for (j=nbits-1; j >= 0; j--)
1302             {
1303               _gcry_mpi_ec_dup_point (result, result, ctx);
1304               if (mpi_test_bit (scalar, j))
1305                 _gcry_mpi_ec_add_points (result, result, point, ctx);
1306             }
1307         }
1308       return;
1309     }
1310   else if (ctx->model == MPI_EC_MONTGOMERY)
1311     {
1312       unsigned int nbits;
1313       int j;
1314       mpi_point_struct p1_, p2_;
1315       mpi_point_t q1, q2, prd, sum;
1316       unsigned long sw;
1317
1318       /* Compute scalar point multiplication with Montgomery Ladder.
1319          Note that we don't use Y-coordinate in the points at all.
1320          RESULT->Y will be filled by zero.  */
1321
1322       nbits = mpi_get_nbits (scalar);
1323       point_init (&p1);
1324       point_init (&p2);
1325       point_init (&p1_);
1326       point_init (&p2_);
1327       mpi_set_ui (p1.x, 1);
1328       mpi_free (p2.x);
1329       p2.x  = mpi_copy (point->x);
1330       mpi_set_ui (p2.z, 1);
1331
1332       point_resize (&p1, ctx);
1333       point_resize (&p2, ctx);
1334       point_resize (&p1_, ctx);
1335       point_resize (&p2_, ctx);
1336
1337       q1 = &p1;
1338       q2 = &p2;
1339       prd = &p1_;
1340       sum = &p2_;
1341
1342       for (j=nbits-1; j >= 0; j--)
1343         {
1344           mpi_point_t t;
1345
1346           sw = mpi_test_bit (scalar, j);
1347           point_swap_cond (q1, q2, sw, ctx);
1348           montgomery_ladder (prd, sum, q1, q2, point->x, ctx);
1349           point_swap_cond (prd, sum, sw, ctx);
1350           t = q1;  q1 = prd;  prd = t;
1351           t = q2;  q2 = sum;  sum = t;
1352         }
1353
1354       mpi_clear (result->y);
1355       sw = (nbits & 1);
1356       point_swap_cond (&p1, &p1_, sw, ctx);
1357
1358       if (p1.z->nlimbs == 0)
1359         {
1360           mpi_set_ui (result->x, 1);
1361           mpi_set_ui (result->z, 0);
1362         }
1363       else
1364         {
1365           z1 = mpi_new (0);
1366           ec_invm (z1, p1.z, ctx);
1367           ec_mulm (result->x, p1.x, z1, ctx);
1368           mpi_set_ui (result->z, 1);
1369           mpi_free (z1);
1370         }
1371
1372       point_free (&p1);
1373       point_free (&p2);
1374       point_free (&p1_);
1375       point_free (&p2_);
1376       return;
1377     }
1378
1379   x1 = mpi_alloc_like (ctx->p);
1380   y1 = mpi_alloc_like (ctx->p);
1381   h  = mpi_alloc_like (ctx->p);
1382   k  = mpi_copy (scalar);
1383   yy = mpi_copy (point->y);
1384
1385   if ( mpi_has_sign (k) )
1386     {
1387       k->sign = 0;
1388       ec_invm (yy, yy, ctx);
1389     }
1390
1391   if (!mpi_cmp_ui (point->z, 1))
1392     {
1393       mpi_set (x1, point->x);
1394       mpi_set (y1, yy);
1395     }
1396   else
1397     {
1398       gcry_mpi_t z2, z3;
1399
1400       z2 = mpi_alloc_like (ctx->p);
1401       z3 = mpi_alloc_like (ctx->p);
1402       ec_mulm (z2, point->z, point->z, ctx);
1403       ec_mulm (z3, point->z, z2, ctx);
1404       ec_invm (z2, z2, ctx);
1405       ec_mulm (x1, point->x, z2, ctx);
1406       ec_invm (z3, z3, ctx);
1407       ec_mulm (y1, yy, z3, ctx);
1408       mpi_free (z2);
1409       mpi_free (z3);
1410     }
1411   z1 = mpi_copy (mpi_const (MPI_C_ONE));
1412
1413   mpi_mul (h, k, mpi_const (MPI_C_THREE)); /* h = 3k */
1414   loops = mpi_get_nbits (h);
1415   if (loops < 2)
1416     {
1417       /* If SCALAR is zero, the above mpi_mul sets H to zero and thus
1418          LOOPs will be zero.  To avoid an underflow of I in the main
1419          loop we set LOOP to 2 and the result to (0,0,0).  */
1420       loops = 2;
1421       mpi_clear (result->x);
1422       mpi_clear (result->y);
1423       mpi_clear (result->z);
1424     }
1425   else
1426     {
1427       mpi_set (result->x, point->x);
1428       mpi_set (result->y, yy);
1429       mpi_set (result->z, point->z);
1430     }
1431   mpi_free (yy); yy = NULL;
1432
1433   p1.x = x1; x1 = NULL;
1434   p1.y = y1; y1 = NULL;
1435   p1.z = z1; z1 = NULL;
1436   point_init (&p2);
1437   point_init (&p1inv);
1438
1439   /* Invert point: y = p - y mod p  */
1440   point_set (&p1inv, &p1);
1441   ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
1442
1443   for (i=loops-2; i > 0; i--)
1444     {
1445       _gcry_mpi_ec_dup_point (result, result, ctx);
1446       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
1447         {
1448           point_set (&p2, result);
1449           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
1450         }
1451       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
1452         {
1453           point_set (&p2, result);
1454           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
1455         }
1456     }
1457
1458   point_free (&p1);
1459   point_free (&p2);
1460   point_free (&p1inv);
1461   mpi_free (h);
1462   mpi_free (k);
1463 }
1464
1465
1466 /* Return true if POINT is on the curve described by CTX.  */
1467 int
1468 _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
1469 {
1470   int res = 0;
1471   gcry_mpi_t x, y, w;
1472
1473   x = mpi_new (0);
1474   y = mpi_new (0);
1475   w = mpi_new (0);
1476
1477   switch (ctx->model)
1478     {
1479     case MPI_EC_WEIERSTRASS:
1480       {
1481         gcry_mpi_t xxx;
1482
1483         if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1484           goto leave;
1485
1486         xxx = mpi_new (0);
1487
1488         /* y^2 == x^3 + a·x + b */
1489         ec_pow2 (y, y, ctx);
1490
1491         ec_pow3 (xxx, x, ctx);
1492         ec_mulm (w, ctx->a, x, ctx);
1493         ec_addm (w, w, ctx->b, ctx);
1494         ec_addm (w, w, xxx, ctx);
1495
1496         if (!mpi_cmp (y, w))
1497           res = 1;
1498
1499         _gcry_mpi_release (xxx);
1500       }
1501       break;
1502     case MPI_EC_MONTGOMERY:
1503       {
1504 #define xx y
1505         /* With Montgomery curve, only X-coordinate is valid.  */
1506         if (_gcry_mpi_ec_get_affine (x, NULL, point, ctx))
1507           goto leave;
1508
1509         /* The equation is: b * y^2 == x^3 + a · x^2 + x */
1510         /* We check if right hand is quadratic residue or not by
1511            Euler's criterion.  */
1512         /* CTX->A has (a-2)/4 and CTX->B has b^-1 */
1513         ec_mulm (w, ctx->a, mpi_const (MPI_C_FOUR), ctx);
1514         ec_addm (w, w, mpi_const (MPI_C_TWO), ctx);
1515         ec_mulm (w, w, x, ctx);
1516         ec_pow2 (xx, x, ctx);
1517         ec_addm (w, w, xx, ctx);
1518         ec_addm (w, w, mpi_const (MPI_C_ONE), ctx);
1519         ec_mulm (w, w, x, ctx);
1520         ec_mulm (w, w, ctx->b, ctx);
1521 #undef xx
1522         /* Compute Euler's criterion: w^(p-1)/2 */
1523 #define p_minus1 y
1524         ec_subm (p_minus1, ctx->p, mpi_const (MPI_C_ONE), ctx);
1525         mpi_rshift (p_minus1, p_minus1, 1);
1526         ec_powm (w, w, p_minus1, ctx);
1527
1528         res = !mpi_cmp_ui (w, 1);
1529 #undef p_minus1
1530       }
1531       break;
1532     case MPI_EC_EDWARDS:
1533       {
1534         if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1535           goto leave;
1536
1537         /* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
1538         ec_pow2 (x, x, ctx);
1539         ec_pow2 (y, y, ctx);
1540         if (ctx->dialect == ECC_DIALECT_ED25519)
1541           {
1542             mpi_set (w, x);
1543             _gcry_mpi_neg (w, w);
1544           }
1545         else
1546           ec_mulm (w, ctx->a, x, ctx);
1547         ec_addm (w, w, y, ctx);
1548         ec_subm (w, w, mpi_const (MPI_C_ONE), ctx);
1549         ec_mulm (x, x, y, ctx);
1550         ec_mulm (x, x, ctx->b, ctx);
1551         ec_subm (w, w, x, ctx);
1552         if (!mpi_cmp_ui (w, 0))
1553           res = 1;
1554       }
1555       break;
1556     }
1557
1558  leave:
1559   _gcry_mpi_release (w);
1560   _gcry_mpi_release (x);
1561   _gcry_mpi_release (y);
1562
1563   return res;
1564 }