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