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