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