6a968a570ffdc7b9ef0ef2f398869f5f636b0147
[libgcrypt.git] / mpi / ec.c
1 /* ec.c -  Elliptic Curve functions
2  * Copyright (C) 2007 Free Software Foundation, Inc.
3  * Copyright (C) 2013 g10 Code GmbH
4  *
5  * This file is part of Libgcrypt.
6  *
7  * Libgcrypt is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU Lesser General Public License as
9  * published by the Free Software Foundation; either version 2.1 of
10  * the License, or (at your option) any later version.
11  *
12  * Libgcrypt is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with this program; if not, see <http://www.gnu.org/licenses/>.
19  */
20
21 #include <config.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <errno.h>
25
26 #include "mpi-internal.h"
27 #include "longlong.h"
28 #include "g10lib.h"
29 #include "context.h"
30 #include "ec-context.h"
31 #include "ec-internal.h"
32
33
34 #define point_init(a)  _gcry_mpi_point_init ((a))
35 #define point_free(a)  _gcry_mpi_point_free_parts ((a))
36
37
38 /* Print a point using the log fucntions.  If CTX is not NULL affine
39    coordinates will be printed.  */
40 void
41 _gcry_mpi_point_log (const char *name, mpi_point_t point, mpi_ec_t ctx)
42 {
43   gcry_mpi_t x, y;
44   char buf[100];
45
46   if (!point)
47     {
48       snprintf (buf, sizeof buf - 1, "%s.*", name);
49       log_mpidump (buf, NULL);
50       return;
51     }
52   snprintf (buf, sizeof buf - 1, "%s.X", name);
53
54   if (ctx)
55     {
56       x = mpi_new (0);
57       y = mpi_new (0);
58     }
59   if (!ctx || _gcry_mpi_ec_get_affine (x, y, point, ctx))
60     {
61       log_mpidump (buf, point->x);
62       buf[strlen(buf)-1] = 'Y';
63       log_mpidump (buf, point->y);
64       buf[strlen(buf)-1] = 'Z';
65       log_mpidump (buf, point->z);
66     }
67   else
68     {
69       buf[strlen(buf)-1] = 'x';
70       log_mpidump (buf, x);
71       buf[strlen(buf)-1] = 'y';
72       log_mpidump (buf, y);
73
74     }
75   if (ctx)
76     {
77       _gcry_mpi_release (x);
78       _gcry_mpi_release (y);
79     }
80 }
81
82
83 /* Create a new point option.  NBITS gives the size in bits of one
84    coordinate; it is only used to pre-allocate some resources and
85    might also be passed as 0 to use a default value.  */
86 mpi_point_t
87 _gcry_mpi_point_new (unsigned int nbits)
88 {
89   mpi_point_t p;
90
91   (void)nbits;  /* Currently not used.  */
92
93   p = gcry_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       gcry_free (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 /* Set the projective coordinates from POINT into X, Y, and Z.  If a
143    coordinate is not required, X, Y, or Z may be passed as NULL.  */
144 void
145 _gcry_mpi_point_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
146                      mpi_point_t point)
147 {
148   if (x)
149     mpi_set (x, point->x);
150   if (y)
151     mpi_set (y, point->y);
152   if (z)
153     mpi_set (z, point->z);
154 }
155
156
157 /* Set the projective coordinates from POINT into X, Y, and Z and
158    release POINT.  If a coordinate is not required, X, Y, or Z may be
159    passed as NULL.  */
160 void
161 _gcry_mpi_point_snatch_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
162                             mpi_point_t point)
163 {
164   mpi_snatch (x, point->x);
165   mpi_snatch (y, point->y);
166   mpi_snatch (z, point->z);
167   gcry_free (point);
168 }
169
170
171 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
172    coordinate is given as NULL, the value 0 is stored into point.  If
173    POINT is given as NULL a new point object is allocated.  Returns
174    POINT or the newly allocated point object. */
175 mpi_point_t
176 _gcry_mpi_point_set (mpi_point_t point,
177                      gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
178 {
179   if (!point)
180     point = mpi_point_new (0);
181
182   if (x)
183     mpi_set (point->x, x);
184   else
185     mpi_clear (point->x);
186   if (y)
187     mpi_set (point->y, y);
188   else
189     mpi_clear (point->y);
190   if (z)
191     mpi_set (point->z, z);
192   else
193     mpi_clear (point->z);
194
195   return point;
196 }
197
198
199 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
200    coordinate is given as NULL, the value 0 is stored into point.  If
201    POINT is given as NULL a new point object is allocated.  The
202    coordinates X, Y, and Z are released.  Returns POINT or the newly
203    allocated point object. */
204 mpi_point_t
205 _gcry_mpi_point_snatch_set (mpi_point_t point,
206                             gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
207 {
208   if (!point)
209     point = mpi_point_new (0);
210
211   if (x)
212     mpi_snatch (point->x, x);
213   else
214     mpi_clear (point->x);
215   if (y)
216     mpi_snatch (point->y, y);
217   else
218     mpi_clear (point->y);
219   if (z)
220     mpi_snatch (point->z, z);
221   else
222     mpi_clear (point->z);
223
224   return point;
225 }
226
227
228 /* W = W mod P.  */
229 static void
230 ec_mod (gcry_mpi_t w, mpi_ec_t ec)
231 {
232   if (0 && ec->dialect == ECC_DIALECT_ED25519)
233     _gcry_mpi_ec_ed25519_mod (w);
234   else if (ec->t.p_barrett)
235     _gcry_mpi_mod_barrett (w, w, ec->t.p_barrett);
236   else
237     _gcry_mpi_mod (w, w, ec->p);
238 }
239
240 static void
241 ec_addm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
242 {
243   mpi_add (w, u, v);
244   ec_mod (w, ctx);
245 }
246
247 static void
248 ec_subm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ec)
249 {
250   (void)ec;
251   mpi_sub (w, u, v);
252   /*ec_mod (w, ec);*/
253 }
254
255 static void
256 ec_mulm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
257 {
258   mpi_mul (w, u, v);
259   ec_mod (w, ctx);
260 }
261
262 /* W = 2 * U mod P.  */
263 static void
264 ec_mul2 (gcry_mpi_t w, gcry_mpi_t u, mpi_ec_t ctx)
265 {
266   mpi_lshift (w, u, 1);
267   ec_mod (w, ctx);
268 }
269
270 static void
271 ec_powm (gcry_mpi_t w, const gcry_mpi_t b, const gcry_mpi_t e,
272          mpi_ec_t ctx)
273 {
274   mpi_powm (w, b, e, ctx->p);
275   /* _gcry_mpi_abs (w); */
276 }
277
278
279 /* Shortcut for
280      ec_powm (B, B, mpi_const (MPI_C_TWO), ctx);
281    for easier optimization.  */
282 static void
283 ec_pow2 (gcry_mpi_t w, const gcry_mpi_t b, mpi_ec_t ctx)
284 {
285   /* Using mpi_mul is slightly faster (at least on amd64).  */
286   /* mpi_powm (w, b, mpi_const (MPI_C_TWO), ctx->p); */
287   ec_mulm (w, b, b, ctx);
288 }
289
290
291 /* Shortcut for
292      ec_powm (B, B, mpi_const (MPI_C_THREE), ctx);
293    for easier optimization.  */
294 static void
295 ec_pow3 (gcry_mpi_t w, const gcry_mpi_t b, mpi_ec_t ctx)
296 {
297   mpi_powm (w, b, mpi_const (MPI_C_THREE), ctx->p);
298 }
299
300
301 static void
302 ec_invm (gcry_mpi_t x, gcry_mpi_t a, mpi_ec_t ctx)
303 {
304   if (!mpi_invm (x, a, ctx->p))
305     {
306       log_error ("ec_invm: inverse does not exist:\n");
307       log_mpidump ("  a", a);
308       log_mpidump ("  p", ctx->p);
309     }
310 }
311
312
313 /* Force recomputation of all helper variables.  */
314 void
315 _gcry_mpi_ec_get_reset (mpi_ec_t ec)
316 {
317   ec->t.valid.a_is_pminus3 = 0;
318   ec->t.valid.two_inv_p = 0;
319 }
320
321
322 /* Accessor for helper variable.  */
323 static int
324 ec_get_a_is_pminus3 (mpi_ec_t ec)
325 {
326   gcry_mpi_t tmp;
327
328   if (!ec->t.valid.a_is_pminus3)
329     {
330       ec->t.valid.a_is_pminus3 = 1;
331       tmp = mpi_alloc_like (ec->p);
332       mpi_sub_ui (tmp, ec->p, 3);
333       ec->t.a_is_pminus3 = !mpi_cmp (ec->a, tmp);
334       mpi_free (tmp);
335     }
336
337   return ec->t.a_is_pminus3;
338 }
339
340
341 /* Accessor for helper variable.  */
342 static gcry_mpi_t
343 ec_get_two_inv_p (mpi_ec_t ec)
344 {
345   if (!ec->t.valid.two_inv_p)
346     {
347       ec->t.valid.two_inv_p = 1;
348       if (!ec->t.two_inv_p)
349         ec->t.two_inv_p = mpi_alloc (0);
350       ec_invm (ec->t.two_inv_p, mpi_const (MPI_C_TWO), ec);
351     }
352   return ec->t.two_inv_p;
353 }
354
355
356
357 /* This function initialized a context for elliptic curve based on the
358    field GF(p).  P is the prime specifying this field, A is the first
359    coefficient.  CTX is expected to be zeroized.  */
360 static void
361 ec_p_init (mpi_ec_t ctx, enum gcry_mpi_ec_models model,
362            enum ecc_dialects dialect,
363            int flags,
364            gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
365 {
366   int i;
367   static int use_barrett;
368
369   if (!use_barrett)
370     {
371       if (getenv ("GCRYPT_BARRETT"))
372         use_barrett = 1;
373       else
374         use_barrett = -1;
375     }
376
377   /* Fixme: Do we want to check some constraints? e.g.  a < p  */
378
379   ctx->model = model;
380   ctx->dialect = dialect;
381   ctx->flags = flags;
382   if (dialect == ECC_DIALECT_ED25519)
383     ctx->nbits = 256;
384   else
385     ctx->nbits = mpi_get_nbits (p);
386   ctx->p = mpi_copy (p);
387   ctx->a = mpi_copy (a);
388   ctx->b = mpi_copy (b);
389
390   ctx->t.p_barrett = use_barrett > 0? _gcry_mpi_barrett_init (ctx->p, 0):NULL;
391
392   _gcry_mpi_ec_get_reset (ctx);
393
394   /* Allocate scratch variables.  */
395   for (i=0; i< DIM(ctx->t.scratch); i++)
396     ctx->t.scratch[i] = mpi_alloc_like (ctx->p);
397
398   /* Prepare for fast reduction.  */
399   /* FIXME: need a test for NIST values.  However it does not gain us
400      any real advantage, for 384 bits it is actually slower than using
401      mpi_mulm.  */
402 /*   ctx->nist_nbits = mpi_get_nbits (ctx->p); */
403 /*   if (ctx->nist_nbits == 192) */
404 /*     { */
405 /*       for (i=0; i < 4; i++) */
406 /*         ctx->s[i] = mpi_new (192); */
407 /*       ctx->c    = mpi_new (192*2); */
408 /*     } */
409 /*   else if (ctx->nist_nbits == 384) */
410 /*     { */
411 /*       for (i=0; i < 10; i++) */
412 /*         ctx->s[i] = mpi_new (384); */
413 /*       ctx->c    = mpi_new (384*2); */
414 /*     } */
415 }
416
417
418 static void
419 ec_deinit (void *opaque)
420 {
421   mpi_ec_t ctx = opaque;
422   int i;
423
424   _gcry_mpi_barrett_free (ctx->t.p_barrett);
425
426   /* Domain parameter.  */
427   mpi_free (ctx->p);
428   mpi_free (ctx->a);
429   mpi_free (ctx->b);
430   _gcry_mpi_point_release (ctx->G);
431   mpi_free (ctx->n);
432
433   /* The key.  */
434   _gcry_mpi_point_release (ctx->Q);
435   mpi_free (ctx->d);
436
437   /* Private data of ec.c.  */
438   mpi_free (ctx->t.two_inv_p);
439
440   for (i=0; i< DIM(ctx->t.scratch); i++)
441     mpi_free (ctx->t.scratch[i]);
442
443 /*   if (ctx->nist_nbits == 192) */
444 /*     { */
445 /*       for (i=0; i < 4; i++) */
446 /*         mpi_free (ctx->s[i]); */
447 /*       mpi_free (ctx->c); */
448 /*     } */
449 /*   else if (ctx->nist_nbits == 384) */
450 /*     { */
451 /*       for (i=0; i < 10; i++) */
452 /*         mpi_free (ctx->s[i]); */
453 /*       mpi_free (ctx->c); */
454 /*     } */
455 }
456
457
458 /* This function returns a new context for elliptic curve based on the
459    field GF(p).  P is the prime specifying this field, A is the first
460    coefficient, B is the second coefficient, and MODEL is the model
461    for the curve.  This function is only used within Libgcrypt and not
462    part of the public API.
463
464    This context needs to be released using _gcry_mpi_ec_free.  */
465 mpi_ec_t
466 _gcry_mpi_ec_p_internal_new (enum gcry_mpi_ec_models model,
467                              enum ecc_dialects dialect,
468                              int flags,
469                              gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
470 {
471   mpi_ec_t ctx;
472
473   ctx = gcry_xcalloc (1, sizeof *ctx);
474   ec_p_init (ctx, model, dialect, flags, p, a, b);
475
476   return ctx;
477 }
478
479
480 /* This is a variant of _gcry_mpi_ec_p_internal_new which returns an
481    public context and does some error checking on the supplied
482    arguments.  On success the new context is stored at R_CTX and 0 is
483    returned; on error NULL is stored at R_CTX and an error code is
484    returned.
485
486    The context needs to be released using gcry_ctx_release.  */
487 gpg_err_code_t
488 _gcry_mpi_ec_p_new (gcry_ctx_t *r_ctx,
489                     enum gcry_mpi_ec_models model,
490                     enum ecc_dialects dialect,
491                     int flags,
492                     gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
493 {
494   gcry_ctx_t ctx;
495   mpi_ec_t ec;
496
497   *r_ctx = NULL;
498   if (!p || !a || !mpi_cmp_ui (a, 0))
499     return GPG_ERR_EINVAL;
500
501   ctx = _gcry_ctx_alloc (CONTEXT_TYPE_EC, sizeof *ec, ec_deinit);
502   if (!ctx)
503     return gpg_err_code_from_syserror ();
504   ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
505   ec_p_init (ec, model, dialect, flags, p, a, b);
506
507   *r_ctx = ctx;
508   return 0;
509 }
510
511
512 void
513 _gcry_mpi_ec_free (mpi_ec_t ctx)
514 {
515   if (ctx)
516     {
517       ec_deinit (ctx);
518       gcry_free (ctx);
519     }
520 }
521
522
523 gcry_mpi_t
524 _gcry_mpi_ec_get_mpi (const char *name, gcry_ctx_t ctx, int copy)
525 {
526   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
527
528   return _gcry_ecc_get_mpi (name, ec, copy);
529 }
530
531
532 gcry_mpi_point_t
533 _gcry_mpi_ec_get_point (const char *name, gcry_ctx_t ctx, int copy)
534 {
535   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
536
537   (void)copy;  /* Not used.  */
538
539   return _gcry_ecc_get_point (name, ec);
540 }
541
542
543 gpg_err_code_t
544 _gcry_mpi_ec_set_mpi (const char *name, gcry_mpi_t newvalue,
545                       gcry_ctx_t ctx)
546 {
547   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
548
549   return _gcry_ecc_set_mpi (name, newvalue, ec);
550 }
551
552
553 gpg_err_code_t
554 _gcry_mpi_ec_set_point (const char *name, gcry_mpi_point_t newvalue,
555                         gcry_ctx_t ctx)
556 {
557   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
558
559   return _gcry_ecc_set_point (name, newvalue, ec);
560 }
561
562
563 /* Compute the affine coordinates from the projective coordinates in
564    POINT.  Set them into X and Y.  If one coordinate is not required,
565    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
566    on success or !0 if POINT is at infinity.  */
567 int
568 _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
569                          mpi_ec_t ctx)
570 {
571   if (!mpi_cmp_ui (point->z, 0))
572     return -1;
573
574   switch (ctx->model)
575     {
576     case MPI_EC_WEIERSTRASS: /* Using Jacobian coordinates.  */
577       {
578         gcry_mpi_t z1, z2, z3;
579
580         z1 = mpi_new (0);
581         z2 = mpi_new (0);
582         ec_invm (z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
583         ec_mulm (z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
584
585         if (x)
586           ec_mulm (x, point->x, z2, ctx);
587
588         if (y)
589           {
590             z3 = mpi_new (0);
591             ec_mulm (z3, z2, z1, ctx);      /* z3 = z^(-3) mod p  */
592             ec_mulm (y, point->y, z3, ctx);
593             mpi_free (z3);
594           }
595
596         mpi_free (z2);
597         mpi_free (z1);
598       }
599       return 0;
600
601     case MPI_EC_MONTGOMERY:
602       {
603         log_fatal ("%s: %s not yet supported\n",
604                    "_gcry_mpi_ec_get_affine", "Montgomery");
605       }
606       return -1;
607
608     case MPI_EC_TWISTEDEDWARDS:
609       {
610         gcry_mpi_t z;
611
612         z = mpi_new (0);
613         ec_invm (z, point->z, ctx);
614
615         if (x)
616           ec_mulm (x, point->x, z, ctx);
617         if (y)
618           ec_mulm (y, point->y, z, ctx);
619
620         _gcry_mpi_release (z);
621       }
622       return 0;
623
624     default:
625       return -1;
626     }
627 }
628
629
630 \f
631 /*  RESULT = 2 * POINT  (Weierstrass version). */
632 static void
633 dup_point_weierstrass (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
634 {
635 #define x3 (result->x)
636 #define y3 (result->y)
637 #define z3 (result->z)
638 #define t1 (ctx->t.scratch[0])
639 #define t2 (ctx->t.scratch[1])
640 #define t3 (ctx->t.scratch[2])
641 #define l1 (ctx->t.scratch[3])
642 #define l2 (ctx->t.scratch[4])
643 #define l3 (ctx->t.scratch[5])
644
645   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
646     {
647       /* P_y == 0 || P_z == 0 => [1:1:0] */
648       mpi_set_ui (x3, 1);
649       mpi_set_ui (y3, 1);
650       mpi_set_ui (z3, 0);
651     }
652   else
653     {
654       if (ec_get_a_is_pminus3 (ctx))  /* Use the faster case.  */
655         {
656           /* L1 = 3(X - Z^2)(X + Z^2) */
657           /*                          T1: used for Z^2. */
658           /*                          T2: used for the right term.  */
659           ec_pow2 (t1, point->z, ctx);
660           ec_subm (l1, point->x, t1, ctx);
661           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
662           ec_addm (t2, point->x, t1, ctx);
663           ec_mulm (l1, l1, t2, ctx);
664         }
665       else /* Standard case. */
666         {
667           /* L1 = 3X^2 + aZ^4 */
668           /*                          T1: used for aZ^4. */
669           ec_pow2 (l1, point->x, ctx);
670           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
671           ec_powm (t1, point->z, mpi_const (MPI_C_FOUR), ctx);
672           ec_mulm (t1, t1, ctx->a, ctx);
673           ec_addm (l1, l1, t1, ctx);
674         }
675       /* Z3 = 2YZ */
676       ec_mulm (z3, point->y, point->z, ctx);
677       ec_mul2 (z3, z3, ctx);
678
679       /* L2 = 4XY^2 */
680       /*                              T2: used for Y2; required later. */
681       ec_pow2 (t2, point->y, ctx);
682       ec_mulm (l2, t2, point->x, ctx);
683       ec_mulm (l2, l2, mpi_const (MPI_C_FOUR), ctx);
684
685       /* X3 = L1^2 - 2L2 */
686       /*                              T1: used for L2^2. */
687       ec_pow2 (x3, l1, ctx);
688       ec_mul2 (t1, l2, ctx);
689       ec_subm (x3, x3, t1, ctx);
690
691       /* L3 = 8Y^4 */
692       /*                              T2: taken from above. */
693       ec_pow2 (t2, t2, ctx);
694       ec_mulm (l3, t2, mpi_const (MPI_C_EIGHT), ctx);
695
696       /* Y3 = L1(L2 - X3) - L3 */
697       ec_subm (y3, l2, x3, ctx);
698       ec_mulm (y3, y3, l1, ctx);
699       ec_subm (y3, y3, l3, ctx);
700     }
701
702 #undef x3
703 #undef y3
704 #undef z3
705 #undef t1
706 #undef t2
707 #undef t3
708 #undef l1
709 #undef l2
710 #undef l3
711 }
712
713
714 /*  RESULT = 2 * POINT  (Montgomery version). */
715 static void
716 dup_point_montgomery (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
717 {
718   (void)result;
719   (void)point;
720   (void)ctx;
721   log_fatal ("%s: %s not yet supported\n",
722              "_gcry_mpi_ec_dup_point", "Montgomery");
723 }
724
725
726 /*  RESULT = 2 * POINT  (Twisted Edwards version). */
727 static void
728 dup_point_twistededwards (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
729 {
730 #define X1 (point->x)
731 #define Y1 (point->y)
732 #define Z1 (point->z)
733 #define X3 (result->x)
734 #define Y3 (result->y)
735 #define Z3 (result->z)
736 #define B (ctx->t.scratch[0])
737 #define C (ctx->t.scratch[1])
738 #define D (ctx->t.scratch[2])
739 #define E (ctx->t.scratch[3])
740 #define F (ctx->t.scratch[4])
741 #define H (ctx->t.scratch[5])
742 #define J (ctx->t.scratch[6])
743
744   /* Compute: (X_3 : Y_3 : Z_3) = 2( X_1 : Y_1 : Z_1 ) */
745
746   /* B = (X_1 + Y_1)^2  */
747   ec_addm (B, X1, Y1, ctx);
748   ec_pow2 (B, B, ctx);
749
750   /* C = X_1^2 */
751   /* D = Y_1^2 */
752   ec_pow2 (C, X1, ctx);
753   ec_pow2 (D, Y1, ctx);
754
755   /* E = aC */
756   if (ctx->dialect == ECC_DIALECT_ED25519)
757     {
758       mpi_set (E, C);
759       _gcry_mpi_neg (E, E);
760     }
761   else
762     ec_mulm (E, ctx->a, C, ctx);
763
764   /* F = E + D */
765   ec_addm (F, E, D, ctx);
766
767   /* H = Z_1^2 */
768   ec_pow2 (H, Z1, ctx);
769
770   /* J = F - 2H */
771   ec_mul2 (J, H, ctx);
772   ec_subm (J, F, J, ctx);
773
774   /* X_3 = (B - C - D) · J */
775   ec_subm (X3, B, C, ctx);
776   ec_subm (X3, X3, D, ctx);
777   ec_mulm (X3, X3, J, ctx);
778
779   /* Y_3 = F · (E - D) */
780   ec_subm (Y3, E, D, ctx);
781   ec_mulm (Y3, Y3, F, ctx);
782
783   /* Z_3 = F · J */
784   ec_mulm (Z3, F, J, ctx);
785
786 #undef X1
787 #undef Y1
788 #undef Z1
789 #undef X3
790 #undef Y3
791 #undef Z3
792 #undef B
793 #undef C
794 #undef D
795 #undef E
796 #undef F
797 #undef H
798 #undef J
799 }
800
801
802 /*  RESULT = 2 * POINT  */
803 void
804 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
805 {
806   switch (ctx->model)
807     {
808     case MPI_EC_WEIERSTRASS:
809       dup_point_weierstrass (result, point, ctx);
810       break;
811     case MPI_EC_MONTGOMERY:
812       dup_point_montgomery (result, point, ctx);
813       break;
814     case MPI_EC_TWISTEDEDWARDS:
815       dup_point_twistededwards (result, point, ctx);
816       break;
817     }
818 }
819
820
821 /* RESULT = P1 + P2  (Weierstrass version).*/
822 static void
823 add_points_weierstrass (mpi_point_t result,
824                         mpi_point_t p1, mpi_point_t p2,
825                         mpi_ec_t ctx)
826 {
827 #define x1 (p1->x    )
828 #define y1 (p1->y    )
829 #define z1 (p1->z    )
830 #define x2 (p2->x    )
831 #define y2 (p2->y    )
832 #define z2 (p2->z    )
833 #define x3 (result->x)
834 #define y3 (result->y)
835 #define z3 (result->z)
836 #define l1 (ctx->t.scratch[0])
837 #define l2 (ctx->t.scratch[1])
838 #define l3 (ctx->t.scratch[2])
839 #define l4 (ctx->t.scratch[3])
840 #define l5 (ctx->t.scratch[4])
841 #define l6 (ctx->t.scratch[5])
842 #define l7 (ctx->t.scratch[6])
843 #define l8 (ctx->t.scratch[7])
844 #define l9 (ctx->t.scratch[8])
845 #define t1 (ctx->t.scratch[9])
846 #define t2 (ctx->t.scratch[10])
847
848   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
849     {
850       /* Same point; need to call the duplicate function.  */
851       _gcry_mpi_ec_dup_point (result, p1, ctx);
852     }
853   else if (!mpi_cmp_ui (z1, 0))
854     {
855       /* P1 is at infinity.  */
856       mpi_set (x3, p2->x);
857       mpi_set (y3, p2->y);
858       mpi_set (z3, p2->z);
859     }
860   else if (!mpi_cmp_ui (z2, 0))
861     {
862       /* P2 is at infinity.  */
863       mpi_set (x3, p1->x);
864       mpi_set (y3, p1->y);
865       mpi_set (z3, p1->z);
866     }
867   else
868     {
869       int z1_is_one = !mpi_cmp_ui (z1, 1);
870       int z2_is_one = !mpi_cmp_ui (z2, 1);
871
872       /* l1 = x1 z2^2  */
873       /* l2 = x2 z1^2  */
874       if (z2_is_one)
875         mpi_set (l1, x1);
876       else
877         {
878           ec_pow2 (l1, z2, ctx);
879           ec_mulm (l1, l1, x1, ctx);
880         }
881       if (z1_is_one)
882         mpi_set (l2, x2);
883       else
884         {
885           ec_pow2 (l2, z1, ctx);
886           ec_mulm (l2, l2, x2, ctx);
887         }
888       /* l3 = l1 - l2 */
889       ec_subm (l3, l1, l2, ctx);
890       /* l4 = y1 z2^3  */
891       ec_powm (l4, z2, mpi_const (MPI_C_THREE), ctx);
892       ec_mulm (l4, l4, y1, ctx);
893       /* l5 = y2 z1^3  */
894       ec_powm (l5, z1, mpi_const (MPI_C_THREE), ctx);
895       ec_mulm (l5, l5, y2, ctx);
896       /* l6 = l4 - l5  */
897       ec_subm (l6, l4, l5, ctx);
898
899       if (!mpi_cmp_ui (l3, 0))
900         {
901           if (!mpi_cmp_ui (l6, 0))
902             {
903               /* P1 and P2 are the same - use duplicate function.  */
904               _gcry_mpi_ec_dup_point (result, p1, ctx);
905             }
906           else
907             {
908               /* P1 is the inverse of P2.  */
909               mpi_set_ui (x3, 1);
910               mpi_set_ui (y3, 1);
911               mpi_set_ui (z3, 0);
912             }
913         }
914       else
915         {
916           /* l7 = l1 + l2  */
917           ec_addm (l7, l1, l2, ctx);
918           /* l8 = l4 + l5  */
919           ec_addm (l8, l4, l5, ctx);
920           /* z3 = z1 z2 l3  */
921           ec_mulm (z3, z1, z2, ctx);
922           ec_mulm (z3, z3, l3, ctx);
923           /* x3 = l6^2 - l7 l3^2  */
924           ec_pow2 (t1, l6, ctx);
925           ec_pow2 (t2, l3, ctx);
926           ec_mulm (t2, t2, l7, ctx);
927           ec_subm (x3, t1, t2, ctx);
928           /* l9 = l7 l3^2 - 2 x3  */
929           ec_mul2 (t1, x3, ctx);
930           ec_subm (l9, t2, t1, ctx);
931           /* y3 = (l9 l6 - l8 l3^3)/2  */
932           ec_mulm (l9, l9, l6, ctx);
933           ec_powm (t1, l3, mpi_const (MPI_C_THREE), ctx); /* fixme: Use saved value*/
934           ec_mulm (t1, t1, l8, ctx);
935           ec_subm (y3, l9, t1, ctx);
936           ec_mulm (y3, y3, ec_get_two_inv_p (ctx), ctx);
937         }
938     }
939
940 #undef x1
941 #undef y1
942 #undef z1
943 #undef x2
944 #undef y2
945 #undef z2
946 #undef x3
947 #undef y3
948 #undef z3
949 #undef l1
950 #undef l2
951 #undef l3
952 #undef l4
953 #undef l5
954 #undef l6
955 #undef l7
956 #undef l8
957 #undef l9
958 #undef t1
959 #undef t2
960 }
961
962
963 /* RESULT = P1 + P2  (Montgomery version).*/
964 static void
965 add_points_montgomery (mpi_point_t result,
966                        mpi_point_t p1, mpi_point_t p2,
967                        mpi_ec_t ctx)
968 {
969   (void)result;
970   (void)p1;
971   (void)p2;
972   (void)ctx;
973   log_fatal ("%s: %s not yet supported\n",
974              "_gcry_mpi_ec_add_points", "Montgomery");
975 }
976
977
978 /* RESULT = P1 + P2  (Twisted Edwards version).*/
979 static void
980 add_points_twistededwards (mpi_point_t result,
981                            mpi_point_t p1, mpi_point_t p2,
982                            mpi_ec_t ctx)
983 {
984 #define X1 (p1->x)
985 #define Y1 (p1->y)
986 #define Z1 (p1->z)
987 #define X2 (p2->x)
988 #define Y2 (p2->y)
989 #define Z2 (p2->z)
990 #define X3 (result->x)
991 #define Y3 (result->y)
992 #define Z3 (result->z)
993 #define A (ctx->t.scratch[0])
994 #define B (ctx->t.scratch[1])
995 #define C (ctx->t.scratch[2])
996 #define D (ctx->t.scratch[3])
997 #define E (ctx->t.scratch[4])
998 #define F (ctx->t.scratch[5])
999 #define G (ctx->t.scratch[6])
1000 #define tmp (ctx->t.scratch[7])
1001
1002   /* Compute: (X_3 : Y_3 : Z_3) = (X_1 : Y_1 : Z_1) + (X_2 : Y_2 : Z_3)  */
1003
1004   /* A = Z1 · Z2 */
1005   ec_mulm (A, Z1, Z2, ctx);
1006
1007   /* B = A^2 */
1008   ec_pow2 (B, A, ctx);
1009
1010   /* C = X1 · X2 */
1011   ec_mulm (C, X1, X2, ctx);
1012
1013   /* D = Y1 · Y2 */
1014   ec_mulm (D, Y1, Y2, ctx);
1015
1016   /* E = d · C · D */
1017   ec_mulm (E, ctx->b, C, ctx);
1018   ec_mulm (E, E, D, ctx);
1019
1020   /* F = B - E */
1021   ec_subm (F, B, E, ctx);
1022
1023   /* G = B + E */
1024   ec_addm (G, B, E, ctx);
1025
1026   /* X_3 = A · F · ((X_1 + Y_1) · (X_2 + Y_2) - C - D) */
1027   ec_addm (tmp, X1, Y1, ctx);
1028   ec_addm (X3, X2, Y2, ctx);
1029   ec_mulm (X3, X3, tmp, ctx);
1030   ec_subm (X3, X3, C, ctx);
1031   ec_subm (X3, X3, D, ctx);
1032   ec_mulm (X3, X3, F, ctx);
1033   ec_mulm (X3, X3, A, ctx);
1034
1035   /* Y_3 = A · G · (D - aC) */
1036   if (ctx->dialect == ECC_DIALECT_ED25519)
1037     {
1038       /* Using ec_addm (Y3, D, C, ctx) is possible but a litte bit
1039          slower because a subm does currently skip the mod step.  */
1040       mpi_set (Y3, C);
1041       _gcry_mpi_neg (Y3, Y3);
1042       ec_subm (Y3, D, Y3, ctx);
1043     }
1044   else
1045     {
1046       ec_mulm (Y3, ctx->a, C, ctx);
1047       ec_subm (Y3, D, Y3, ctx);
1048     }
1049   ec_mulm (Y3, Y3, G, ctx);
1050   ec_mulm (Y3, Y3, A, ctx);
1051
1052   /* Z_3 = F · G */
1053   ec_mulm (Z3, F, G, ctx);
1054
1055
1056 #undef X1
1057 #undef Y1
1058 #undef Z1
1059 #undef X2
1060 #undef Y2
1061 #undef Z2
1062 #undef X3
1063 #undef Y3
1064 #undef Z3
1065 #undef A
1066 #undef B
1067 #undef C
1068 #undef D
1069 #undef E
1070 #undef F
1071 #undef G
1072 #undef tmp
1073 }
1074
1075
1076 /* RESULT = P1 + P2 */
1077 void
1078 _gcry_mpi_ec_add_points (mpi_point_t result,
1079                          mpi_point_t p1, mpi_point_t p2,
1080                          mpi_ec_t ctx)
1081 {
1082   switch (ctx->model)
1083     {
1084     case MPI_EC_WEIERSTRASS:
1085       add_points_weierstrass (result, p1, p2, ctx);
1086       break;
1087     case MPI_EC_MONTGOMERY:
1088       add_points_montgomery (result, p1, p2, ctx);
1089       break;
1090     case MPI_EC_TWISTEDEDWARDS:
1091       add_points_twistededwards (result, p1, p2, ctx);
1092       break;
1093     }
1094 }
1095
1096
1097 /* Scalar point multiplication - the main function for ECC.  If takes
1098    an integer SCALAR and a POINT as well as the usual context CTX.
1099    RESULT will be set to the resulting point. */
1100 void
1101 _gcry_mpi_ec_mul_point (mpi_point_t result,
1102                         gcry_mpi_t scalar, mpi_point_t point,
1103                         mpi_ec_t ctx)
1104 {
1105   gcry_mpi_t x1, y1, z1, k, h, yy;
1106   unsigned int i, loops;
1107   mpi_point_struct p1, p2, p1inv;
1108
1109   if (ctx->model == MPI_EC_TWISTEDEDWARDS)
1110     {
1111       /* Simple left to right binary method.  GECC Algorithm 3.27 */
1112       unsigned int nbits;
1113       int j;
1114
1115       nbits = mpi_get_nbits (scalar);
1116       mpi_set_ui (result->x, 0);
1117       mpi_set_ui (result->y, 1);
1118       mpi_set_ui (result->z, 1);
1119
1120       if (mpi_is_secure (scalar))
1121         {
1122           /* If SCALAR is in secure memory we assume that it is the
1123              secret key we use constant time operation.  */
1124           mpi_point_struct tmppnt;
1125
1126           point_init (&tmppnt);
1127           for (j=nbits-1; j >= 0; j--)
1128             {
1129               _gcry_mpi_ec_dup_point (result, result, ctx);
1130               _gcry_mpi_ec_add_points (&tmppnt, result, point, ctx);
1131               if (mpi_test_bit (scalar, j))
1132                 point_set (result, &tmppnt);
1133             }
1134           point_free (&tmppnt);
1135         }
1136       else
1137         {
1138           for (j=nbits-1; j >= 0; j--)
1139             {
1140               _gcry_mpi_ec_dup_point (result, result, ctx);
1141               if (mpi_test_bit (scalar, j))
1142                 _gcry_mpi_ec_add_points (result, result, point, ctx);
1143             }
1144         }
1145       return;
1146     }
1147
1148   x1 = mpi_alloc_like (ctx->p);
1149   y1 = mpi_alloc_like (ctx->p);
1150   h  = mpi_alloc_like (ctx->p);
1151   k  = mpi_copy (scalar);
1152   yy = mpi_copy (point->y);
1153
1154   if ( mpi_has_sign (k) )
1155     {
1156       k->sign = 0;
1157       ec_invm (yy, yy, ctx);
1158     }
1159
1160   if (!mpi_cmp_ui (point->z, 1))
1161     {
1162       mpi_set (x1, point->x);
1163       mpi_set (y1, yy);
1164     }
1165   else
1166     {
1167       gcry_mpi_t z2, z3;
1168
1169       z2 = mpi_alloc_like (ctx->p);
1170       z3 = mpi_alloc_like (ctx->p);
1171       ec_mulm (z2, point->z, point->z, ctx);
1172       ec_mulm (z3, point->z, z2, ctx);
1173       ec_invm (z2, z2, ctx);
1174       ec_mulm (x1, point->x, z2, ctx);
1175       ec_invm (z3, z3, ctx);
1176       ec_mulm (y1, yy, z3, ctx);
1177       mpi_free (z2);
1178       mpi_free (z3);
1179     }
1180   z1 = mpi_copy (mpi_const (MPI_C_ONE));
1181
1182   mpi_mul (h, k, mpi_const (MPI_C_THREE)); /* h = 3k */
1183   loops = mpi_get_nbits (h);
1184   if (loops < 2)
1185     {
1186       /* If SCALAR is zero, the above mpi_mul sets H to zero and thus
1187          LOOPs will be zero.  To avoid an underflow of I in the main
1188          loop we set LOOP to 2 and the result to (0,0,0).  */
1189       loops = 2;
1190       mpi_clear (result->x);
1191       mpi_clear (result->y);
1192       mpi_clear (result->z);
1193     }
1194   else
1195     {
1196       mpi_set (result->x, point->x);
1197       mpi_set (result->y, yy);
1198       mpi_set (result->z, point->z);
1199     }
1200   mpi_free (yy); yy = NULL;
1201
1202   p1.x = x1; x1 = NULL;
1203   p1.y = y1; y1 = NULL;
1204   p1.z = z1; z1 = NULL;
1205   point_init (&p2);
1206   point_init (&p1inv);
1207
1208   for (i=loops-2; i > 0; i--)
1209     {
1210       _gcry_mpi_ec_dup_point (result, result, ctx);
1211       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
1212         {
1213           point_set (&p2, result);
1214           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
1215         }
1216       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
1217         {
1218           point_set (&p2, result);
1219           /* Invert point: y = p - y mod p  */
1220           point_set (&p1inv, &p1);
1221           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
1222           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
1223         }
1224     }
1225
1226   point_free (&p1);
1227   point_free (&p2);
1228   point_free (&p1inv);
1229   mpi_free (h);
1230   mpi_free (k);
1231 }
1232
1233
1234 /* Return true if POINT is on the curve described by CTX.  */
1235 int
1236 _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
1237 {
1238   int res = 0;
1239   gcry_mpi_t x, y, w;
1240
1241   x = mpi_new (0);
1242   y = mpi_new (0);
1243   w = mpi_new (0);
1244
1245   if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1246     return 0;
1247
1248   switch (ctx->model)
1249     {
1250     case MPI_EC_WEIERSTRASS:
1251       {
1252         gcry_mpi_t xxx = mpi_new (0);
1253
1254         /* y^2 == x^3 + a·x + b */
1255         ec_pow2 (y, y, ctx);
1256
1257         ec_pow3 (xxx, x, ctx);
1258         ec_mulm (w, ctx->a, x, ctx);
1259         ec_addm (w, w, ctx->b, ctx);
1260         ec_addm (w, w, xxx, ctx);
1261
1262         if (!mpi_cmp (y, w))
1263           res = 1;
1264
1265         _gcry_mpi_release (xxx);
1266       }
1267       break;
1268     case MPI_EC_MONTGOMERY:
1269       log_fatal ("%s: %s not yet supported\n",
1270                  "_gcry_mpi_ec_curve_point", "Montgomery");
1271       break;
1272     case MPI_EC_TWISTEDEDWARDS:
1273       {
1274         /* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
1275         ec_pow2 (x, x, ctx);
1276         ec_pow2 (y, y, ctx);
1277         if (ctx->dialect == ECC_DIALECT_ED25519)
1278           {
1279             mpi_set (w, x);
1280             _gcry_mpi_neg (w, w);
1281           }
1282         else
1283           ec_mulm (w, ctx->a, x, ctx);
1284         ec_addm (w, w, y, ctx);
1285         ec_subm (w, w, mpi_const (MPI_C_ONE), ctx);
1286         ec_mulm (x, x, y, ctx);
1287         ec_mulm (x, x, ctx->b, ctx);
1288         ec_subm (w, w, x, ctx);
1289         if (!mpi_cmp_ui (w, 0))
1290           res = 1;
1291       }
1292       break;
1293     }
1294
1295   _gcry_mpi_release (w);
1296   _gcry_mpi_release (x);
1297   _gcry_mpi_release (y);
1298
1299   return res;
1300 }