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