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