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