ecc: multiplication of Edwards curve to be constant-time.
[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 /* Compute the affine coordinates from the projective coordinates in
593    POINT.  Set them into X and Y.  If one coordinate is not required,
594    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
595    on success or !0 if POINT is at infinity.  */
596 int
597 _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
598                          mpi_ec_t ctx)
599 {
600   if (!mpi_cmp_ui (point->z, 0))
601     return -1;
602
603   switch (ctx->model)
604     {
605     case MPI_EC_WEIERSTRASS: /* Using Jacobian coordinates.  */
606       {
607         gcry_mpi_t z1, z2, z3;
608
609         z1 = mpi_new (0);
610         z2 = mpi_new (0);
611         ec_invm (z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
612         ec_mulm (z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
613
614         if (x)
615           ec_mulm (x, point->x, z2, ctx);
616
617         if (y)
618           {
619             z3 = mpi_new (0);
620             ec_mulm (z3, z2, z1, ctx);      /* z3 = z^(-3) mod p  */
621             ec_mulm (y, point->y, z3, ctx);
622             mpi_free (z3);
623           }
624
625         mpi_free (z2);
626         mpi_free (z1);
627       }
628       return 0;
629
630     case MPI_EC_MONTGOMERY:
631       {
632         if (x)
633           mpi_set (x, point->x);
634
635         if (y)
636           {
637             log_fatal ("%s: Getting Y-coordinate on %s is not supported\n",
638                        "_gcry_mpi_ec_get_affine", "Montgomery");
639             return -1;
640           }
641       }
642       return 0;
643
644     case MPI_EC_EDWARDS:
645       {
646         gcry_mpi_t z;
647
648         z = mpi_new (0);
649         ec_invm (z, point->z, ctx);
650
651         if (x)
652           ec_mulm (x, point->x, z, ctx);
653         if (y)
654           ec_mulm (y, point->y, z, ctx);
655
656         _gcry_mpi_release (z);
657       }
658       return 0;
659
660     default:
661       return -1;
662     }
663 }
664
665
666 \f
667 /*  RESULT = 2 * POINT  (Weierstrass version). */
668 static void
669 dup_point_weierstrass (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
670 {
671 #define x3 (result->x)
672 #define y3 (result->y)
673 #define z3 (result->z)
674 #define t1 (ctx->t.scratch[0])
675 #define t2 (ctx->t.scratch[1])
676 #define t3 (ctx->t.scratch[2])
677 #define l1 (ctx->t.scratch[3])
678 #define l2 (ctx->t.scratch[4])
679 #define l3 (ctx->t.scratch[5])
680
681   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
682     {
683       /* P_y == 0 || P_z == 0 => [1:1:0] */
684       mpi_set_ui (x3, 1);
685       mpi_set_ui (y3, 1);
686       mpi_set_ui (z3, 0);
687     }
688   else
689     {
690       if (ec_get_a_is_pminus3 (ctx))  /* Use the faster case.  */
691         {
692           /* L1 = 3(X - Z^2)(X + Z^2) */
693           /*                          T1: used for Z^2. */
694           /*                          T2: used for the right term.  */
695           ec_pow2 (t1, point->z, ctx);
696           ec_subm (l1, point->x, t1, ctx);
697           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
698           ec_addm (t2, point->x, t1, ctx);
699           ec_mulm (l1, l1, t2, ctx);
700         }
701       else /* Standard case. */
702         {
703           /* L1 = 3X^2 + aZ^4 */
704           /*                          T1: used for aZ^4. */
705           ec_pow2 (l1, point->x, ctx);
706           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
707           ec_powm (t1, point->z, mpi_const (MPI_C_FOUR), ctx);
708           ec_mulm (t1, t1, ctx->a, ctx);
709           ec_addm (l1, l1, t1, ctx);
710         }
711       /* Z3 = 2YZ */
712       ec_mulm (z3, point->y, point->z, ctx);
713       ec_mul2 (z3, z3, ctx);
714
715       /* L2 = 4XY^2 */
716       /*                              T2: used for Y2; required later. */
717       ec_pow2 (t2, point->y, ctx);
718       ec_mulm (l2, t2, point->x, ctx);
719       ec_mulm (l2, l2, mpi_const (MPI_C_FOUR), ctx);
720
721       /* X3 = L1^2 - 2L2 */
722       /*                              T1: used for L2^2. */
723       ec_pow2 (x3, l1, ctx);
724       ec_mul2 (t1, l2, ctx);
725       ec_subm (x3, x3, t1, ctx);
726
727       /* L3 = 8Y^4 */
728       /*                              T2: taken from above. */
729       ec_pow2 (t2, t2, ctx);
730       ec_mulm (l3, t2, mpi_const (MPI_C_EIGHT), ctx);
731
732       /* Y3 = L1(L2 - X3) - L3 */
733       ec_subm (y3, l2, x3, ctx);
734       ec_mulm (y3, y3, l1, ctx);
735       ec_subm (y3, y3, l3, ctx);
736     }
737
738 #undef x3
739 #undef y3
740 #undef z3
741 #undef t1
742 #undef t2
743 #undef t3
744 #undef l1
745 #undef l2
746 #undef l3
747 }
748
749
750 /*  RESULT = 2 * POINT  (Montgomery version). */
751 static void
752 dup_point_montgomery (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
753 {
754   (void)result;
755   (void)point;
756   (void)ctx;
757   log_fatal ("%s: %s not yet supported\n",
758              "_gcry_mpi_ec_dup_point", "Montgomery");
759 }
760
761
762 /*  RESULT = 2 * POINT  (Twisted Edwards version). */
763 static void
764 dup_point_edwards (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
765 {
766 #define X1 (point->x)
767 #define Y1 (point->y)
768 #define Z1 (point->z)
769 #define X3 (result->x)
770 #define Y3 (result->y)
771 #define Z3 (result->z)
772 #define B (ctx->t.scratch[0])
773 #define C (ctx->t.scratch[1])
774 #define D (ctx->t.scratch[2])
775 #define E (ctx->t.scratch[3])
776 #define F (ctx->t.scratch[4])
777 #define H (ctx->t.scratch[5])
778 #define J (ctx->t.scratch[6])
779
780   /* Compute: (X_3 : Y_3 : Z_3) = 2( X_1 : Y_1 : Z_1 ) */
781
782   /* B = (X_1 + Y_1)^2  */
783   ec_addm (B, X1, Y1, ctx);
784   ec_pow2 (B, B, ctx);
785
786   /* C = X_1^2 */
787   /* D = Y_1^2 */
788   ec_pow2 (C, X1, ctx);
789   ec_pow2 (D, Y1, ctx);
790
791   /* E = aC */
792   if (ctx->dialect == ECC_DIALECT_ED25519)
793     {
794       mpi_set (E, C);
795       _gcry_mpi_neg (E, E);
796     }
797   else
798     ec_mulm (E, ctx->a, C, ctx);
799
800   /* F = E + D */
801   ec_addm (F, E, D, ctx);
802
803   /* H = Z_1^2 */
804   ec_pow2 (H, Z1, ctx);
805
806   /* J = F - 2H */
807   ec_mul2 (J, H, ctx);
808   ec_subm (J, F, J, ctx);
809
810   /* X_3 = (B - C - D) · J */
811   ec_subm (X3, B, C, ctx);
812   ec_subm (X3, X3, D, ctx);
813   ec_mulm (X3, X3, J, ctx);
814
815   /* Y_3 = F · (E - D) */
816   ec_subm (Y3, E, D, ctx);
817   ec_mulm (Y3, Y3, F, ctx);
818
819   /* Z_3 = F · J */
820   ec_mulm (Z3, F, J, ctx);
821
822 #undef X1
823 #undef Y1
824 #undef Z1
825 #undef X3
826 #undef Y3
827 #undef Z3
828 #undef B
829 #undef C
830 #undef D
831 #undef E
832 #undef F
833 #undef H
834 #undef J
835 }
836
837
838 /*  RESULT = 2 * POINT  */
839 void
840 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
841 {
842   switch (ctx->model)
843     {
844     case MPI_EC_WEIERSTRASS:
845       dup_point_weierstrass (result, point, ctx);
846       break;
847     case MPI_EC_MONTGOMERY:
848       dup_point_montgomery (result, point, ctx);
849       break;
850     case MPI_EC_EDWARDS:
851       dup_point_edwards (result, point, ctx);
852       break;
853     }
854 }
855
856
857 /* RESULT = P1 + P2  (Weierstrass version).*/
858 static void
859 add_points_weierstrass (mpi_point_t result,
860                         mpi_point_t p1, mpi_point_t p2,
861                         mpi_ec_t ctx)
862 {
863 #define x1 (p1->x    )
864 #define y1 (p1->y    )
865 #define z1 (p1->z    )
866 #define x2 (p2->x    )
867 #define y2 (p2->y    )
868 #define z2 (p2->z    )
869 #define x3 (result->x)
870 #define y3 (result->y)
871 #define z3 (result->z)
872 #define l1 (ctx->t.scratch[0])
873 #define l2 (ctx->t.scratch[1])
874 #define l3 (ctx->t.scratch[2])
875 #define l4 (ctx->t.scratch[3])
876 #define l5 (ctx->t.scratch[4])
877 #define l6 (ctx->t.scratch[5])
878 #define l7 (ctx->t.scratch[6])
879 #define l8 (ctx->t.scratch[7])
880 #define l9 (ctx->t.scratch[8])
881 #define t1 (ctx->t.scratch[9])
882 #define t2 (ctx->t.scratch[10])
883
884   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
885     {
886       /* Same point; need to call the duplicate function.  */
887       _gcry_mpi_ec_dup_point (result, p1, ctx);
888     }
889   else if (!mpi_cmp_ui (z1, 0))
890     {
891       /* P1 is at infinity.  */
892       mpi_set (x3, p2->x);
893       mpi_set (y3, p2->y);
894       mpi_set (z3, p2->z);
895     }
896   else if (!mpi_cmp_ui (z2, 0))
897     {
898       /* P2 is at infinity.  */
899       mpi_set (x3, p1->x);
900       mpi_set (y3, p1->y);
901       mpi_set (z3, p1->z);
902     }
903   else
904     {
905       int z1_is_one = !mpi_cmp_ui (z1, 1);
906       int z2_is_one = !mpi_cmp_ui (z2, 1);
907
908       /* l1 = x1 z2^2  */
909       /* l2 = x2 z1^2  */
910       if (z2_is_one)
911         mpi_set (l1, x1);
912       else
913         {
914           ec_pow2 (l1, z2, ctx);
915           ec_mulm (l1, l1, x1, ctx);
916         }
917       if (z1_is_one)
918         mpi_set (l2, x2);
919       else
920         {
921           ec_pow2 (l2, z1, ctx);
922           ec_mulm (l2, l2, x2, ctx);
923         }
924       /* l3 = l1 - l2 */
925       ec_subm (l3, l1, l2, ctx);
926       /* l4 = y1 z2^3  */
927       ec_powm (l4, z2, mpi_const (MPI_C_THREE), ctx);
928       ec_mulm (l4, l4, y1, ctx);
929       /* l5 = y2 z1^3  */
930       ec_powm (l5, z1, mpi_const (MPI_C_THREE), ctx);
931       ec_mulm (l5, l5, y2, ctx);
932       /* l6 = l4 - l5  */
933       ec_subm (l6, l4, l5, ctx);
934
935       if (!mpi_cmp_ui (l3, 0))
936         {
937           if (!mpi_cmp_ui (l6, 0))
938             {
939               /* P1 and P2 are the same - use duplicate function.  */
940               _gcry_mpi_ec_dup_point (result, p1, ctx);
941             }
942           else
943             {
944               /* P1 is the inverse of P2.  */
945               mpi_set_ui (x3, 1);
946               mpi_set_ui (y3, 1);
947               mpi_set_ui (z3, 0);
948             }
949         }
950       else
951         {
952           /* l7 = l1 + l2  */
953           ec_addm (l7, l1, l2, ctx);
954           /* l8 = l4 + l5  */
955           ec_addm (l8, l4, l5, ctx);
956           /* z3 = z1 z2 l3  */
957           ec_mulm (z3, z1, z2, ctx);
958           ec_mulm (z3, z3, l3, ctx);
959           /* x3 = l6^2 - l7 l3^2  */
960           ec_pow2 (t1, l6, ctx);
961           ec_pow2 (t2, l3, ctx);
962           ec_mulm (t2, t2, l7, ctx);
963           ec_subm (x3, t1, t2, ctx);
964           /* l9 = l7 l3^2 - 2 x3  */
965           ec_mul2 (t1, x3, ctx);
966           ec_subm (l9, t2, t1, ctx);
967           /* y3 = (l9 l6 - l8 l3^3)/2  */
968           ec_mulm (l9, l9, l6, ctx);
969           ec_powm (t1, l3, mpi_const (MPI_C_THREE), ctx); /* fixme: Use saved value*/
970           ec_mulm (t1, t1, l8, ctx);
971           ec_subm (y3, l9, t1, ctx);
972           ec_mulm (y3, y3, ec_get_two_inv_p (ctx), ctx);
973         }
974     }
975
976 #undef x1
977 #undef y1
978 #undef z1
979 #undef x2
980 #undef y2
981 #undef z2
982 #undef x3
983 #undef y3
984 #undef z3
985 #undef l1
986 #undef l2
987 #undef l3
988 #undef l4
989 #undef l5
990 #undef l6
991 #undef l7
992 #undef l8
993 #undef l9
994 #undef t1
995 #undef t2
996 }
997
998
999 /* RESULT = P1 + P2  (Montgomery version).*/
1000 static void
1001 add_points_montgomery (mpi_point_t result,
1002                        mpi_point_t p1, mpi_point_t p2,
1003                        mpi_ec_t ctx)
1004 {
1005   (void)result;
1006   (void)p1;
1007   (void)p2;
1008   (void)ctx;
1009   log_fatal ("%s: %s not yet supported\n",
1010              "_gcry_mpi_ec_add_points", "Montgomery");
1011 }
1012
1013
1014 /* RESULT = P1 + P2  (Twisted Edwards version).*/
1015 static void
1016 add_points_edwards (mpi_point_t result,
1017                     mpi_point_t p1, mpi_point_t p2,
1018                     mpi_ec_t ctx)
1019 {
1020 #define X1 (p1->x)
1021 #define Y1 (p1->y)
1022 #define Z1 (p1->z)
1023 #define X2 (p2->x)
1024 #define Y2 (p2->y)
1025 #define Z2 (p2->z)
1026 #define X3 (result->x)
1027 #define Y3 (result->y)
1028 #define Z3 (result->z)
1029 #define A (ctx->t.scratch[0])
1030 #define B (ctx->t.scratch[1])
1031 #define C (ctx->t.scratch[2])
1032 #define D (ctx->t.scratch[3])
1033 #define E (ctx->t.scratch[4])
1034 #define F (ctx->t.scratch[5])
1035 #define G (ctx->t.scratch[6])
1036 #define tmp (ctx->t.scratch[7])
1037
1038   /* Compute: (X_3 : Y_3 : Z_3) = (X_1 : Y_1 : Z_1) + (X_2 : Y_2 : Z_3)  */
1039
1040   /* A = Z1 · Z2 */
1041   ec_mulm (A, Z1, Z2, ctx);
1042
1043   /* B = A^2 */
1044   ec_pow2 (B, A, ctx);
1045
1046   /* C = X1 · X2 */
1047   ec_mulm (C, X1, X2, ctx);
1048
1049   /* D = Y1 · Y2 */
1050   ec_mulm (D, Y1, Y2, ctx);
1051
1052   /* E = d · C · D */
1053   ec_mulm (E, ctx->b, C, ctx);
1054   ec_mulm (E, E, D, ctx);
1055
1056   /* F = B - E */
1057   ec_subm (F, B, E, ctx);
1058
1059   /* G = B + E */
1060   ec_addm (G, B, E, ctx);
1061
1062   /* X_3 = A · F · ((X_1 + Y_1) · (X_2 + Y_2) - C - D) */
1063   ec_addm (tmp, X1, Y1, ctx);
1064   ec_addm (X3, X2, Y2, ctx);
1065   ec_mulm (X3, X3, tmp, ctx);
1066   ec_subm (X3, X3, C, ctx);
1067   ec_subm (X3, X3, D, ctx);
1068   ec_mulm (X3, X3, F, ctx);
1069   ec_mulm (X3, X3, A, ctx);
1070
1071   /* Y_3 = A · G · (D - aC) */
1072   if (ctx->dialect == ECC_DIALECT_ED25519)
1073     {
1074       /* Using ec_addm (Y3, D, C, ctx) is possible but a litte bit
1075          slower because a subm does currently skip the mod step.  */
1076       mpi_set (Y3, C);
1077       _gcry_mpi_neg (Y3, Y3);
1078       ec_subm (Y3, D, Y3, ctx);
1079     }
1080   else
1081     {
1082       ec_mulm (Y3, ctx->a, C, ctx);
1083       ec_subm (Y3, D, Y3, ctx);
1084     }
1085   ec_mulm (Y3, Y3, G, ctx);
1086   ec_mulm (Y3, Y3, A, ctx);
1087
1088   /* Z_3 = F · G */
1089   ec_mulm (Z3, F, G, ctx);
1090
1091
1092 #undef X1
1093 #undef Y1
1094 #undef Z1
1095 #undef X2
1096 #undef Y2
1097 #undef Z2
1098 #undef X3
1099 #undef Y3
1100 #undef Z3
1101 #undef A
1102 #undef B
1103 #undef C
1104 #undef D
1105 #undef E
1106 #undef F
1107 #undef G
1108 #undef tmp
1109 }
1110
1111
1112 /* Compute a step of Montgomery Ladder (only use X and Z in the point).
1113    Inputs:  P1, P2, and x-coordinate of DIF = P1 - P1.
1114    Outputs: PRD = 2 * P1 and  SUM = P1 + P2. */
1115 static void
1116 montgomery_ladder (mpi_point_t prd, mpi_point_t sum,
1117                    mpi_point_t p1, mpi_point_t p2, gcry_mpi_t dif_x,
1118                    mpi_ec_t ctx)
1119 {
1120   ec_addm (sum->x, p2->x, p2->z, ctx);
1121   ec_subm (p2->z, p2->x, p2->z, ctx);
1122   ec_addm (prd->x, p1->x, p1->z, ctx);
1123   ec_subm (p1->z, p1->x, p1->z, ctx);
1124   ec_mulm (p2->x, p1->z, sum->x, ctx);
1125   ec_mulm (p2->z, prd->x, p2->z, ctx);
1126   ec_pow2 (p1->x, prd->x, ctx);
1127   ec_pow2 (p1->z, p1->z, ctx);
1128   ec_addm (sum->x, p2->x, p2->z, ctx);
1129   ec_subm (p2->z, p2->x, p2->z, ctx);
1130   ec_mulm (prd->x, p1->x, p1->z, ctx);
1131   ec_subm (p1->z, p1->x, p1->z, ctx);
1132   ec_pow2 (sum->x, sum->x, ctx);
1133   ec_pow2 (sum->z, p2->z, ctx);
1134   ec_mulm (prd->z, p1->z, ctx->a, ctx); /* CTX->A: (a-2)/4 */
1135   ec_mulm (sum->z, sum->z, dif_x, ctx);
1136   ec_addm (prd->z, p1->x, prd->z, ctx);
1137   ec_mulm (prd->z, prd->z, p1->z, ctx);
1138 }
1139
1140
1141 /* RESULT = P1 + P2 */
1142 void
1143 _gcry_mpi_ec_add_points (mpi_point_t result,
1144                          mpi_point_t p1, mpi_point_t p2,
1145                          mpi_ec_t ctx)
1146 {
1147   switch (ctx->model)
1148     {
1149     case MPI_EC_WEIERSTRASS:
1150       add_points_weierstrass (result, p1, p2, ctx);
1151       break;
1152     case MPI_EC_MONTGOMERY:
1153       add_points_montgomery (result, p1, p2, ctx);
1154       break;
1155     case MPI_EC_EDWARDS:
1156       add_points_edwards (result, p1, p2, ctx);
1157       break;
1158     }
1159 }
1160
1161
1162 /* RESULT = P1 - P2  (Weierstrass version).*/
1163 static void
1164 sub_points_weierstrass (mpi_point_t result,
1165                         mpi_point_t p1, mpi_point_t p2,
1166                         mpi_ec_t ctx)
1167 {
1168   (void)result;
1169   (void)p1;
1170   (void)p2;
1171   (void)ctx;
1172   log_fatal ("%s: %s not yet supported\n",
1173              "_gcry_mpi_ec_sub_points", "Weierstrass");
1174 }
1175
1176
1177 /* RESULT = P1 - P2  (Montgomery version).*/
1178 static void
1179 sub_points_montgomery (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", "Montgomery");
1189 }
1190
1191
1192 /* RESULT = P1 - P2  (Twisted Edwards version).*/
1193 static void
1194 sub_points_edwards (mpi_point_t result,
1195                     mpi_point_t p1, mpi_point_t p2,
1196                     mpi_ec_t ctx)
1197 {
1198   mpi_point_t p2i = _gcry_mpi_point_new (0);
1199   point_set (p2i, p2);
1200   _gcry_mpi_neg (p2i->x, p2i->x);
1201   add_points_edwards (result, p1, p2i, ctx);
1202   _gcry_mpi_point_release (p2i);
1203 }
1204
1205
1206 /* RESULT = P1 - P2 */
1207 void
1208 _gcry_mpi_ec_sub_points (mpi_point_t result,
1209                          mpi_point_t p1, mpi_point_t p2,
1210                          mpi_ec_t ctx)
1211 {
1212   switch (ctx->model)
1213     {
1214     case MPI_EC_WEIERSTRASS:
1215       sub_points_weierstrass (result, p1, p2, ctx);
1216       break;
1217     case MPI_EC_MONTGOMERY:
1218       sub_points_montgomery (result, p1, p2, ctx);
1219       break;
1220     case MPI_EC_EDWARDS:
1221       sub_points_edwards (result, p1, p2, ctx);
1222       break;
1223     }
1224 }
1225
1226
1227 /* Scalar point multiplication - the main function for ECC.  If takes
1228    an integer SCALAR and a POINT as well as the usual context CTX.
1229    RESULT will be set to the resulting point. */
1230 void
1231 _gcry_mpi_ec_mul_point (mpi_point_t result,
1232                         gcry_mpi_t scalar, mpi_point_t point,
1233                         mpi_ec_t ctx)
1234 {
1235   gcry_mpi_t x1, y1, z1, k, h, yy;
1236   unsigned int i, loops;
1237   mpi_point_struct p1, p2, p1inv;
1238
1239   if (ctx->model == MPI_EC_EDWARDS)
1240     {
1241       /* Simple left to right binary method.  GECC Algorithm 3.27 */
1242       unsigned int nbits;
1243       int j;
1244
1245       nbits = mpi_get_nbits (scalar);
1246       mpi_set_ui (result->x, 0);
1247       mpi_set_ui (result->y, 1);
1248       mpi_set_ui (result->z, 1);
1249
1250       if (mpi_is_secure (scalar))
1251         {
1252           /* If SCALAR is in secure memory we assume that it is the
1253              secret key we use constant time operation.  */
1254           mpi_point_struct tmppnt;
1255
1256           point_init (&tmppnt);
1257           point_resize (result, ctx);
1258           point_resize (&tmppnt, ctx);
1259           for (j=nbits-1; j >= 0; j--)
1260             {
1261               _gcry_mpi_ec_dup_point (result, result, ctx);
1262               _gcry_mpi_ec_add_points (&tmppnt, result, point, ctx);
1263               point_swap_cond (result, &tmppnt, mpi_test_bit (scalar, j), ctx);
1264             }
1265           point_free (&tmppnt);
1266         }
1267       else
1268         {
1269           for (j=nbits-1; j >= 0; j--)
1270             {
1271               _gcry_mpi_ec_dup_point (result, result, ctx);
1272               if (mpi_test_bit (scalar, j))
1273                 _gcry_mpi_ec_add_points (result, result, point, ctx);
1274             }
1275         }
1276       return;
1277     }
1278   else if (ctx->model == MPI_EC_MONTGOMERY)
1279     {
1280       unsigned int nbits;
1281       int j;
1282       mpi_point_struct p1_, p2_;
1283       mpi_point_t q1, q2, prd, sum;
1284       unsigned long sw;
1285
1286       /* Compute scalar point multiplication with Montgomery Ladder.
1287          Note that we don't use Y-coordinate in the points at all.
1288          RESULT->Y will be filled by zero.  */
1289
1290       nbits = mpi_get_nbits (scalar);
1291       point_init (&p1);
1292       point_init (&p2);
1293       point_init (&p1_);
1294       point_init (&p2_);
1295       mpi_set_ui (p1.x, 1);
1296       mpi_free (p2.x);
1297       p2.x  = mpi_copy (point->x);
1298       mpi_set_ui (p2.z, 1);
1299
1300       point_resize (&p1, ctx);
1301       point_resize (&p2, ctx);
1302       point_resize (&p1_, ctx);
1303       point_resize (&p2_, ctx);
1304
1305       q1 = &p1;
1306       q2 = &p2;
1307       prd = &p1_;
1308       sum = &p2_;
1309
1310       for (j=nbits-1; j >= 0; j--)
1311         {
1312           mpi_point_t t;
1313
1314           sw = mpi_test_bit (scalar, j);
1315           point_swap_cond (q1, q2, sw, ctx);
1316           montgomery_ladder (prd, sum, q1, q2, point->x, ctx);
1317           point_swap_cond (prd, sum, sw, ctx);
1318           t = q1;  q1 = prd;  prd = t;
1319           t = q2;  q2 = sum;  sum = t;
1320         }
1321
1322       mpi_clear (result->y);
1323       sw = (nbits & 1);
1324       point_swap_cond (&p1, &p1_, sw, ctx);
1325
1326       if (p1.z->nlimbs == 0)
1327         {
1328           mpi_set_ui (result->x, 1);
1329           mpi_set_ui (result->z, 0);
1330         }
1331       else
1332         {
1333           z1 = mpi_new (0);
1334           ec_invm (z1, p1.z, ctx);
1335           ec_mulm (result->x, p1.x, z1, ctx);
1336           mpi_set_ui (result->z, 1);
1337           mpi_free (z1);
1338         }
1339
1340       point_free (&p1);
1341       point_free (&p2);
1342       point_free (&p1_);
1343       point_free (&p2_);
1344       return;
1345     }
1346
1347   x1 = mpi_alloc_like (ctx->p);
1348   y1 = mpi_alloc_like (ctx->p);
1349   h  = mpi_alloc_like (ctx->p);
1350   k  = mpi_copy (scalar);
1351   yy = mpi_copy (point->y);
1352
1353   if ( mpi_has_sign (k) )
1354     {
1355       k->sign = 0;
1356       ec_invm (yy, yy, ctx);
1357     }
1358
1359   if (!mpi_cmp_ui (point->z, 1))
1360     {
1361       mpi_set (x1, point->x);
1362       mpi_set (y1, yy);
1363     }
1364   else
1365     {
1366       gcry_mpi_t z2, z3;
1367
1368       z2 = mpi_alloc_like (ctx->p);
1369       z3 = mpi_alloc_like (ctx->p);
1370       ec_mulm (z2, point->z, point->z, ctx);
1371       ec_mulm (z3, point->z, z2, ctx);
1372       ec_invm (z2, z2, ctx);
1373       ec_mulm (x1, point->x, z2, ctx);
1374       ec_invm (z3, z3, ctx);
1375       ec_mulm (y1, yy, z3, ctx);
1376       mpi_free (z2);
1377       mpi_free (z3);
1378     }
1379   z1 = mpi_copy (mpi_const (MPI_C_ONE));
1380
1381   mpi_mul (h, k, mpi_const (MPI_C_THREE)); /* h = 3k */
1382   loops = mpi_get_nbits (h);
1383   if (loops < 2)
1384     {
1385       /* If SCALAR is zero, the above mpi_mul sets H to zero and thus
1386          LOOPs will be zero.  To avoid an underflow of I in the main
1387          loop we set LOOP to 2 and the result to (0,0,0).  */
1388       loops = 2;
1389       mpi_clear (result->x);
1390       mpi_clear (result->y);
1391       mpi_clear (result->z);
1392     }
1393   else
1394     {
1395       mpi_set (result->x, point->x);
1396       mpi_set (result->y, yy);
1397       mpi_set (result->z, point->z);
1398     }
1399   mpi_free (yy); yy = NULL;
1400
1401   p1.x = x1; x1 = NULL;
1402   p1.y = y1; y1 = NULL;
1403   p1.z = z1; z1 = NULL;
1404   point_init (&p2);
1405   point_init (&p1inv);
1406
1407   for (i=loops-2; i > 0; i--)
1408     {
1409       _gcry_mpi_ec_dup_point (result, result, ctx);
1410       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
1411         {
1412           point_set (&p2, result);
1413           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
1414         }
1415       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
1416         {
1417           point_set (&p2, result);
1418           /* Invert point: y = p - y mod p  */
1419           point_set (&p1inv, &p1);
1420           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
1421           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
1422         }
1423     }
1424
1425   point_free (&p1);
1426   point_free (&p2);
1427   point_free (&p1inv);
1428   mpi_free (h);
1429   mpi_free (k);
1430 }
1431
1432
1433 /* Return true if POINT is on the curve described by CTX.  */
1434 int
1435 _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
1436 {
1437   int res = 0;
1438   gcry_mpi_t x, y, w;
1439
1440   x = mpi_new (0);
1441   y = mpi_new (0);
1442   w = mpi_new (0);
1443
1444   switch (ctx->model)
1445     {
1446     case MPI_EC_WEIERSTRASS:
1447       {
1448         gcry_mpi_t xxx = mpi_new (0);
1449
1450         if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1451           return 0;
1452
1453         /* y^2 == x^3 + a·x + b */
1454         ec_pow2 (y, y, ctx);
1455
1456         ec_pow3 (xxx, x, ctx);
1457         ec_mulm (w, ctx->a, x, ctx);
1458         ec_addm (w, w, ctx->b, ctx);
1459         ec_addm (w, w, xxx, ctx);
1460
1461         if (!mpi_cmp (y, w))
1462           res = 1;
1463
1464         _gcry_mpi_release (xxx);
1465       }
1466       break;
1467     case MPI_EC_MONTGOMERY:
1468       {
1469 #define xx y
1470         /* With Montgomery curve, only X-coordinate is valid.  */
1471         if (_gcry_mpi_ec_get_affine (x, NULL, point, ctx))
1472           return 0;
1473
1474         /* The equation is: b * y^2 == x^3 + a · x^2 + x */
1475         /* We check if right hand is quadratic residue or not by
1476            Euler's criterion.  */
1477         /* CTX->A has (a-2)/4 and CTX->B has b^-1 */
1478         ec_mulm (w, ctx->a, mpi_const (MPI_C_FOUR), ctx);
1479         ec_addm (w, w, mpi_const (MPI_C_TWO), ctx);
1480         ec_mulm (w, w, x, ctx);
1481         ec_pow2 (xx, x, ctx);
1482         ec_addm (w, w, xx, ctx);
1483         ec_addm (w, w, mpi_const (MPI_C_ONE), ctx);
1484         ec_mulm (w, w, x, ctx);
1485         ec_mulm (w, w, ctx->b, ctx);
1486 #undef xx
1487         /* Compute Euler's criterion: w^(p-1)/2 */
1488 #define p_minus1 y
1489         ec_subm (p_minus1, ctx->p, mpi_const (MPI_C_ONE), ctx);
1490         mpi_rshift (p_minus1, p_minus1, 1);
1491         ec_powm (w, w, p_minus1, ctx);
1492
1493         res = !mpi_cmp_ui (w, 1);
1494 #undef p_minus1
1495       }
1496       break;
1497     case MPI_EC_EDWARDS:
1498       {
1499         if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1500           return 0;
1501
1502         /* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
1503         ec_pow2 (x, x, ctx);
1504         ec_pow2 (y, y, ctx);
1505         if (ctx->dialect == ECC_DIALECT_ED25519)
1506           {
1507             mpi_set (w, x);
1508             _gcry_mpi_neg (w, w);
1509           }
1510         else
1511           ec_mulm (w, ctx->a, x, ctx);
1512         ec_addm (w, w, y, ctx);
1513         ec_subm (w, w, mpi_const (MPI_C_ONE), ctx);
1514         ec_mulm (x, x, y, ctx);
1515         ec_mulm (x, x, ctx->b, ctx);
1516         ec_subm (w, w, x, ctx);
1517         if (!mpi_cmp_ui (w, 0))
1518           res = 1;
1519       }
1520       break;
1521     }
1522
1523   _gcry_mpi_release (w);
1524   _gcry_mpi_release (x);
1525   _gcry_mpi_release (y);
1526
1527   return res;
1528 }