ecc: Add cofactor to domain parameters.
[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 = 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 /* 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   xfree (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   mpi_free (ctx->h);
433
434   /* The key.  */
435   _gcry_mpi_point_release (ctx->Q);
436   mpi_free (ctx->d);
437
438   /* Private data of ec.c.  */
439   mpi_free (ctx->t.two_inv_p);
440
441   for (i=0; i< DIM(ctx->t.scratch); i++)
442     mpi_free (ctx->t.scratch[i]);
443
444 /*   if (ctx->nist_nbits == 192) */
445 /*     { */
446 /*       for (i=0; i < 4; i++) */
447 /*         mpi_free (ctx->s[i]); */
448 /*       mpi_free (ctx->c); */
449 /*     } */
450 /*   else if (ctx->nist_nbits == 384) */
451 /*     { */
452 /*       for (i=0; i < 10; i++) */
453 /*         mpi_free (ctx->s[i]); */
454 /*       mpi_free (ctx->c); */
455 /*     } */
456 }
457
458
459 /* This function returns a new context for elliptic curve based on the
460    field GF(p).  P is the prime specifying this field, A is the first
461    coefficient, B is the second coefficient, and MODEL is the model
462    for the curve.  This function is only used within Libgcrypt and not
463    part of the public API.
464
465    This context needs to be released using _gcry_mpi_ec_free.  */
466 mpi_ec_t
467 _gcry_mpi_ec_p_internal_new (enum gcry_mpi_ec_models model,
468                              enum ecc_dialects dialect,
469                              int flags,
470                              gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
471 {
472   mpi_ec_t ctx;
473
474   ctx = xcalloc (1, sizeof *ctx);
475   ec_p_init (ctx, model, dialect, flags, p, a, b);
476
477   return ctx;
478 }
479
480
481 /* This is a variant of _gcry_mpi_ec_p_internal_new which returns an
482    public context and does some error checking on the supplied
483    arguments.  On success the new context is stored at R_CTX and 0 is
484    returned; on error NULL is stored at R_CTX and an error code is
485    returned.
486
487    The context needs to be released using gcry_ctx_release.  */
488 gpg_err_code_t
489 _gcry_mpi_ec_p_new (gcry_ctx_t *r_ctx,
490                     enum gcry_mpi_ec_models model,
491                     enum ecc_dialects dialect,
492                     int flags,
493                     gcry_mpi_t p, gcry_mpi_t a, gcry_mpi_t b)
494 {
495   gcry_ctx_t ctx;
496   mpi_ec_t ec;
497
498   *r_ctx = NULL;
499   if (!p || !a)
500     return GPG_ERR_EINVAL;
501
502   ctx = _gcry_ctx_alloc (CONTEXT_TYPE_EC, sizeof *ec, ec_deinit);
503   if (!ctx)
504     return gpg_err_code_from_syserror ();
505   ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
506   ec_p_init (ec, model, dialect, flags, p, a, b);
507
508   *r_ctx = ctx;
509   return 0;
510 }
511
512
513 void
514 _gcry_mpi_ec_free (mpi_ec_t ctx)
515 {
516   if (ctx)
517     {
518       ec_deinit (ctx);
519       xfree (ctx);
520     }
521 }
522
523
524 gcry_mpi_t
525 _gcry_mpi_ec_get_mpi (const char *name, gcry_ctx_t ctx, int copy)
526 {
527   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
528
529   return _gcry_ecc_get_mpi (name, ec, copy);
530 }
531
532
533 gcry_mpi_point_t
534 _gcry_mpi_ec_get_point (const char *name, gcry_ctx_t ctx, int copy)
535 {
536   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
537
538   (void)copy;  /* Not used.  */
539
540   return _gcry_ecc_get_point (name, ec);
541 }
542
543
544 gpg_err_code_t
545 _gcry_mpi_ec_set_mpi (const char *name, gcry_mpi_t newvalue,
546                       gcry_ctx_t ctx)
547 {
548   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
549
550   return _gcry_ecc_set_mpi (name, newvalue, ec);
551 }
552
553
554 gpg_err_code_t
555 _gcry_mpi_ec_set_point (const char *name, gcry_mpi_point_t newvalue,
556                         gcry_ctx_t ctx)
557 {
558   mpi_ec_t ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
559
560   return _gcry_ecc_set_point (name, newvalue, ec);
561 }
562
563
564 /* Compute the affine coordinates from the projective coordinates in
565    POINT.  Set them into X and Y.  If one coordinate is not required,
566    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
567    on success or !0 if POINT is at infinity.  */
568 int
569 _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
570                          mpi_ec_t ctx)
571 {
572   if (!mpi_cmp_ui (point->z, 0))
573     return -1;
574
575   switch (ctx->model)
576     {
577     case MPI_EC_WEIERSTRASS: /* Using Jacobian coordinates.  */
578       {
579         gcry_mpi_t z1, z2, z3;
580
581         z1 = mpi_new (0);
582         z2 = mpi_new (0);
583         ec_invm (z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
584         ec_mulm (z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
585
586         if (x)
587           ec_mulm (x, point->x, z2, ctx);
588
589         if (y)
590           {
591             z3 = mpi_new (0);
592             ec_mulm (z3, z2, z1, ctx);      /* z3 = z^(-3) mod p  */
593             ec_mulm (y, point->y, z3, ctx);
594             mpi_free (z3);
595           }
596
597         mpi_free (z2);
598         mpi_free (z1);
599       }
600       return 0;
601
602     case MPI_EC_MONTGOMERY:
603       {
604         log_fatal ("%s: %s not yet supported\n",
605                    "_gcry_mpi_ec_get_affine", "Montgomery");
606       }
607       return -1;
608
609     case MPI_EC_EDWARDS:
610       {
611         gcry_mpi_t z;
612
613         z = mpi_new (0);
614         ec_invm (z, point->z, ctx);
615
616         if (x)
617           ec_mulm (x, point->x, z, ctx);
618         if (y)
619           ec_mulm (y, point->y, z, ctx);
620
621         _gcry_mpi_release (z);
622       }
623       return 0;
624
625     default:
626       return -1;
627     }
628 }
629
630
631 \f
632 /*  RESULT = 2 * POINT  (Weierstrass version). */
633 static void
634 dup_point_weierstrass (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
635 {
636 #define x3 (result->x)
637 #define y3 (result->y)
638 #define z3 (result->z)
639 #define t1 (ctx->t.scratch[0])
640 #define t2 (ctx->t.scratch[1])
641 #define t3 (ctx->t.scratch[2])
642 #define l1 (ctx->t.scratch[3])
643 #define l2 (ctx->t.scratch[4])
644 #define l3 (ctx->t.scratch[5])
645
646   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
647     {
648       /* P_y == 0 || P_z == 0 => [1:1:0] */
649       mpi_set_ui (x3, 1);
650       mpi_set_ui (y3, 1);
651       mpi_set_ui (z3, 0);
652     }
653   else
654     {
655       if (ec_get_a_is_pminus3 (ctx))  /* Use the faster case.  */
656         {
657           /* L1 = 3(X - Z^2)(X + Z^2) */
658           /*                          T1: used for Z^2. */
659           /*                          T2: used for the right term.  */
660           ec_pow2 (t1, point->z, ctx);
661           ec_subm (l1, point->x, t1, ctx);
662           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
663           ec_addm (t2, point->x, t1, ctx);
664           ec_mulm (l1, l1, t2, ctx);
665         }
666       else /* Standard case. */
667         {
668           /* L1 = 3X^2 + aZ^4 */
669           /*                          T1: used for aZ^4. */
670           ec_pow2 (l1, point->x, ctx);
671           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
672           ec_powm (t1, point->z, mpi_const (MPI_C_FOUR), ctx);
673           ec_mulm (t1, t1, ctx->a, ctx);
674           ec_addm (l1, l1, t1, ctx);
675         }
676       /* Z3 = 2YZ */
677       ec_mulm (z3, point->y, point->z, ctx);
678       ec_mul2 (z3, z3, ctx);
679
680       /* L2 = 4XY^2 */
681       /*                              T2: used for Y2; required later. */
682       ec_pow2 (t2, point->y, ctx);
683       ec_mulm (l2, t2, point->x, ctx);
684       ec_mulm (l2, l2, mpi_const (MPI_C_FOUR), ctx);
685
686       /* X3 = L1^2 - 2L2 */
687       /*                              T1: used for L2^2. */
688       ec_pow2 (x3, l1, ctx);
689       ec_mul2 (t1, l2, ctx);
690       ec_subm (x3, x3, t1, ctx);
691
692       /* L3 = 8Y^4 */
693       /*                              T2: taken from above. */
694       ec_pow2 (t2, t2, ctx);
695       ec_mulm (l3, t2, mpi_const (MPI_C_EIGHT), ctx);
696
697       /* Y3 = L1(L2 - X3) - L3 */
698       ec_subm (y3, l2, x3, ctx);
699       ec_mulm (y3, y3, l1, ctx);
700       ec_subm (y3, y3, l3, ctx);
701     }
702
703 #undef x3
704 #undef y3
705 #undef z3
706 #undef t1
707 #undef t2
708 #undef t3
709 #undef l1
710 #undef l2
711 #undef l3
712 }
713
714
715 /*  RESULT = 2 * POINT  (Montgomery version). */
716 static void
717 dup_point_montgomery (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
718 {
719   (void)result;
720   (void)point;
721   (void)ctx;
722   log_fatal ("%s: %s not yet supported\n",
723              "_gcry_mpi_ec_dup_point", "Montgomery");
724 }
725
726
727 /*  RESULT = 2 * POINT  (Twisted Edwards version). */
728 static void
729 dup_point_edwards (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
730 {
731 #define X1 (point->x)
732 #define Y1 (point->y)
733 #define Z1 (point->z)
734 #define X3 (result->x)
735 #define Y3 (result->y)
736 #define Z3 (result->z)
737 #define B (ctx->t.scratch[0])
738 #define C (ctx->t.scratch[1])
739 #define D (ctx->t.scratch[2])
740 #define E (ctx->t.scratch[3])
741 #define F (ctx->t.scratch[4])
742 #define H (ctx->t.scratch[5])
743 #define J (ctx->t.scratch[6])
744
745   /* Compute: (X_3 : Y_3 : Z_3) = 2( X_1 : Y_1 : Z_1 ) */
746
747   /* B = (X_1 + Y_1)^2  */
748   ec_addm (B, X1, Y1, ctx);
749   ec_pow2 (B, B, ctx);
750
751   /* C = X_1^2 */
752   /* D = Y_1^2 */
753   ec_pow2 (C, X1, ctx);
754   ec_pow2 (D, Y1, ctx);
755
756   /* E = aC */
757   if (ctx->dialect == ECC_DIALECT_ED25519)
758     {
759       mpi_set (E, C);
760       _gcry_mpi_neg (E, E);
761     }
762   else
763     ec_mulm (E, ctx->a, C, ctx);
764
765   /* F = E + D */
766   ec_addm (F, E, D, ctx);
767
768   /* H = Z_1^2 */
769   ec_pow2 (H, Z1, ctx);
770
771   /* J = F - 2H */
772   ec_mul2 (J, H, ctx);
773   ec_subm (J, F, J, ctx);
774
775   /* X_3 = (B - C - D) · J */
776   ec_subm (X3, B, C, ctx);
777   ec_subm (X3, X3, D, ctx);
778   ec_mulm (X3, X3, J, ctx);
779
780   /* Y_3 = F · (E - D) */
781   ec_subm (Y3, E, D, ctx);
782   ec_mulm (Y3, Y3, F, ctx);
783
784   /* Z_3 = F · J */
785   ec_mulm (Z3, F, J, ctx);
786
787 #undef X1
788 #undef Y1
789 #undef Z1
790 #undef X3
791 #undef Y3
792 #undef Z3
793 #undef B
794 #undef C
795 #undef D
796 #undef E
797 #undef F
798 #undef H
799 #undef J
800 }
801
802
803 /*  RESULT = 2 * POINT  */
804 void
805 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
806 {
807   switch (ctx->model)
808     {
809     case MPI_EC_WEIERSTRASS:
810       dup_point_weierstrass (result, point, ctx);
811       break;
812     case MPI_EC_MONTGOMERY:
813       dup_point_montgomery (result, point, ctx);
814       break;
815     case MPI_EC_EDWARDS:
816       dup_point_edwards (result, point, ctx);
817       break;
818     }
819 }
820
821
822 /* RESULT = P1 + P2  (Weierstrass version).*/
823 static void
824 add_points_weierstrass (mpi_point_t result,
825                         mpi_point_t p1, mpi_point_t p2,
826                         mpi_ec_t ctx)
827 {
828 #define x1 (p1->x    )
829 #define y1 (p1->y    )
830 #define z1 (p1->z    )
831 #define x2 (p2->x    )
832 #define y2 (p2->y    )
833 #define z2 (p2->z    )
834 #define x3 (result->x)
835 #define y3 (result->y)
836 #define z3 (result->z)
837 #define l1 (ctx->t.scratch[0])
838 #define l2 (ctx->t.scratch[1])
839 #define l3 (ctx->t.scratch[2])
840 #define l4 (ctx->t.scratch[3])
841 #define l5 (ctx->t.scratch[4])
842 #define l6 (ctx->t.scratch[5])
843 #define l7 (ctx->t.scratch[6])
844 #define l8 (ctx->t.scratch[7])
845 #define l9 (ctx->t.scratch[8])
846 #define t1 (ctx->t.scratch[9])
847 #define t2 (ctx->t.scratch[10])
848
849   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
850     {
851       /* Same point; need to call the duplicate function.  */
852       _gcry_mpi_ec_dup_point (result, p1, ctx);
853     }
854   else if (!mpi_cmp_ui (z1, 0))
855     {
856       /* P1 is at infinity.  */
857       mpi_set (x3, p2->x);
858       mpi_set (y3, p2->y);
859       mpi_set (z3, p2->z);
860     }
861   else if (!mpi_cmp_ui (z2, 0))
862     {
863       /* P2 is at infinity.  */
864       mpi_set (x3, p1->x);
865       mpi_set (y3, p1->y);
866       mpi_set (z3, p1->z);
867     }
868   else
869     {
870       int z1_is_one = !mpi_cmp_ui (z1, 1);
871       int z2_is_one = !mpi_cmp_ui (z2, 1);
872
873       /* l1 = x1 z2^2  */
874       /* l2 = x2 z1^2  */
875       if (z2_is_one)
876         mpi_set (l1, x1);
877       else
878         {
879           ec_pow2 (l1, z2, ctx);
880           ec_mulm (l1, l1, x1, ctx);
881         }
882       if (z1_is_one)
883         mpi_set (l2, x2);
884       else
885         {
886           ec_pow2 (l2, z1, ctx);
887           ec_mulm (l2, l2, x2, ctx);
888         }
889       /* l3 = l1 - l2 */
890       ec_subm (l3, l1, l2, ctx);
891       /* l4 = y1 z2^3  */
892       ec_powm (l4, z2, mpi_const (MPI_C_THREE), ctx);
893       ec_mulm (l4, l4, y1, ctx);
894       /* l5 = y2 z1^3  */
895       ec_powm (l5, z1, mpi_const (MPI_C_THREE), ctx);
896       ec_mulm (l5, l5, y2, ctx);
897       /* l6 = l4 - l5  */
898       ec_subm (l6, l4, l5, ctx);
899
900       if (!mpi_cmp_ui (l3, 0))
901         {
902           if (!mpi_cmp_ui (l6, 0))
903             {
904               /* P1 and P2 are the same - use duplicate function.  */
905               _gcry_mpi_ec_dup_point (result, p1, ctx);
906             }
907           else
908             {
909               /* P1 is the inverse of P2.  */
910               mpi_set_ui (x3, 1);
911               mpi_set_ui (y3, 1);
912               mpi_set_ui (z3, 0);
913             }
914         }
915       else
916         {
917           /* l7 = l1 + l2  */
918           ec_addm (l7, l1, l2, ctx);
919           /* l8 = l4 + l5  */
920           ec_addm (l8, l4, l5, ctx);
921           /* z3 = z1 z2 l3  */
922           ec_mulm (z3, z1, z2, ctx);
923           ec_mulm (z3, z3, l3, ctx);
924           /* x3 = l6^2 - l7 l3^2  */
925           ec_pow2 (t1, l6, ctx);
926           ec_pow2 (t2, l3, ctx);
927           ec_mulm (t2, t2, l7, ctx);
928           ec_subm (x3, t1, t2, ctx);
929           /* l9 = l7 l3^2 - 2 x3  */
930           ec_mul2 (t1, x3, ctx);
931           ec_subm (l9, t2, t1, ctx);
932           /* y3 = (l9 l6 - l8 l3^3)/2  */
933           ec_mulm (l9, l9, l6, ctx);
934           ec_powm (t1, l3, mpi_const (MPI_C_THREE), ctx); /* fixme: Use saved value*/
935           ec_mulm (t1, t1, l8, ctx);
936           ec_subm (y3, l9, t1, ctx);
937           ec_mulm (y3, y3, ec_get_two_inv_p (ctx), ctx);
938         }
939     }
940
941 #undef x1
942 #undef y1
943 #undef z1
944 #undef x2
945 #undef y2
946 #undef z2
947 #undef x3
948 #undef y3
949 #undef z3
950 #undef l1
951 #undef l2
952 #undef l3
953 #undef l4
954 #undef l5
955 #undef l6
956 #undef l7
957 #undef l8
958 #undef l9
959 #undef t1
960 #undef t2
961 }
962
963
964 /* RESULT = P1 + P2  (Montgomery version).*/
965 static void
966 add_points_montgomery (mpi_point_t result,
967                        mpi_point_t p1, mpi_point_t p2,
968                        mpi_ec_t ctx)
969 {
970   (void)result;
971   (void)p1;
972   (void)p2;
973   (void)ctx;
974   log_fatal ("%s: %s not yet supported\n",
975              "_gcry_mpi_ec_add_points", "Montgomery");
976 }
977
978
979 /* RESULT = P1 + P2  (Twisted Edwards version).*/
980 static void
981 add_points_edwards (mpi_point_t result,
982                     mpi_point_t p1, mpi_point_t p2,
983                     mpi_ec_t ctx)
984 {
985 #define X1 (p1->x)
986 #define Y1 (p1->y)
987 #define Z1 (p1->z)
988 #define X2 (p2->x)
989 #define Y2 (p2->y)
990 #define Z2 (p2->z)
991 #define X3 (result->x)
992 #define Y3 (result->y)
993 #define Z3 (result->z)
994 #define A (ctx->t.scratch[0])
995 #define B (ctx->t.scratch[1])
996 #define C (ctx->t.scratch[2])
997 #define D (ctx->t.scratch[3])
998 #define E (ctx->t.scratch[4])
999 #define F (ctx->t.scratch[5])
1000 #define G (ctx->t.scratch[6])
1001 #define tmp (ctx->t.scratch[7])
1002
1003   /* Compute: (X_3 : Y_3 : Z_3) = (X_1 : Y_1 : Z_1) + (X_2 : Y_2 : Z_3)  */
1004
1005   /* A = Z1 · Z2 */
1006   ec_mulm (A, Z1, Z2, ctx);
1007
1008   /* B = A^2 */
1009   ec_pow2 (B, A, ctx);
1010
1011   /* C = X1 · X2 */
1012   ec_mulm (C, X1, X2, ctx);
1013
1014   /* D = Y1 · Y2 */
1015   ec_mulm (D, Y1, Y2, ctx);
1016
1017   /* E = d · C · D */
1018   ec_mulm (E, ctx->b, C, ctx);
1019   ec_mulm (E, E, D, ctx);
1020
1021   /* F = B - E */
1022   ec_subm (F, B, E, ctx);
1023
1024   /* G = B + E */
1025   ec_addm (G, B, E, ctx);
1026
1027   /* X_3 = A · F · ((X_1 + Y_1) · (X_2 + Y_2) - C - D) */
1028   ec_addm (tmp, X1, Y1, ctx);
1029   ec_addm (X3, X2, Y2, ctx);
1030   ec_mulm (X3, X3, tmp, ctx);
1031   ec_subm (X3, X3, C, ctx);
1032   ec_subm (X3, X3, D, ctx);
1033   ec_mulm (X3, X3, F, ctx);
1034   ec_mulm (X3, X3, A, ctx);
1035
1036   /* Y_3 = A · G · (D - aC) */
1037   if (ctx->dialect == ECC_DIALECT_ED25519)
1038     {
1039       /* Using ec_addm (Y3, D, C, ctx) is possible but a litte bit
1040          slower because a subm does currently skip the mod step.  */
1041       mpi_set (Y3, C);
1042       _gcry_mpi_neg (Y3, Y3);
1043       ec_subm (Y3, D, Y3, ctx);
1044     }
1045   else
1046     {
1047       ec_mulm (Y3, ctx->a, C, ctx);
1048       ec_subm (Y3, D, Y3, ctx);
1049     }
1050   ec_mulm (Y3, Y3, G, ctx);
1051   ec_mulm (Y3, Y3, A, ctx);
1052
1053   /* Z_3 = F · G */
1054   ec_mulm (Z3, F, G, ctx);
1055
1056
1057 #undef X1
1058 #undef Y1
1059 #undef Z1
1060 #undef X2
1061 #undef Y2
1062 #undef Z2
1063 #undef X3
1064 #undef Y3
1065 #undef Z3
1066 #undef A
1067 #undef B
1068 #undef C
1069 #undef D
1070 #undef E
1071 #undef F
1072 #undef G
1073 #undef tmp
1074 }
1075
1076
1077 /* RESULT = P1 + P2 */
1078 void
1079 _gcry_mpi_ec_add_points (mpi_point_t result,
1080                          mpi_point_t p1, mpi_point_t p2,
1081                          mpi_ec_t ctx)
1082 {
1083   switch (ctx->model)
1084     {
1085     case MPI_EC_WEIERSTRASS:
1086       add_points_weierstrass (result, p1, p2, ctx);
1087       break;
1088     case MPI_EC_MONTGOMERY:
1089       add_points_montgomery (result, p1, p2, ctx);
1090       break;
1091     case MPI_EC_EDWARDS:
1092       add_points_edwards (result, p1, p2, ctx);
1093       break;
1094     }
1095 }
1096
1097
1098 /* Scalar point multiplication - the main function for ECC.  If takes
1099    an integer SCALAR and a POINT as well as the usual context CTX.
1100    RESULT will be set to the resulting point. */
1101 void
1102 _gcry_mpi_ec_mul_point (mpi_point_t result,
1103                         gcry_mpi_t scalar, mpi_point_t point,
1104                         mpi_ec_t ctx)
1105 {
1106   gcry_mpi_t x1, y1, z1, k, h, yy;
1107   unsigned int i, loops;
1108   mpi_point_struct p1, p2, p1inv;
1109
1110   if (ctx->model == MPI_EC_EDWARDS)
1111     {
1112       /* Simple left to right binary method.  GECC Algorithm 3.27 */
1113       unsigned int nbits;
1114       int j;
1115
1116       nbits = mpi_get_nbits (scalar);
1117       mpi_set_ui (result->x, 0);
1118       mpi_set_ui (result->y, 1);
1119       mpi_set_ui (result->z, 1);
1120
1121       if (mpi_is_secure (scalar))
1122         {
1123           /* If SCALAR is in secure memory we assume that it is the
1124              secret key we use constant time operation.  */
1125           mpi_point_struct tmppnt;
1126
1127           point_init (&tmppnt);
1128           for (j=nbits-1; j >= 0; j--)
1129             {
1130               _gcry_mpi_ec_dup_point (result, result, ctx);
1131               _gcry_mpi_ec_add_points (&tmppnt, result, point, ctx);
1132               if (mpi_test_bit (scalar, j))
1133                 point_set (result, &tmppnt);
1134             }
1135           point_free (&tmppnt);
1136         }
1137       else
1138         {
1139           for (j=nbits-1; j >= 0; j--)
1140             {
1141               _gcry_mpi_ec_dup_point (result, result, ctx);
1142               if (mpi_test_bit (scalar, j))
1143                 _gcry_mpi_ec_add_points (result, result, point, ctx);
1144             }
1145         }
1146       return;
1147     }
1148
1149   x1 = mpi_alloc_like (ctx->p);
1150   y1 = mpi_alloc_like (ctx->p);
1151   h  = mpi_alloc_like (ctx->p);
1152   k  = mpi_copy (scalar);
1153   yy = mpi_copy (point->y);
1154
1155   if ( mpi_has_sign (k) )
1156     {
1157       k->sign = 0;
1158       ec_invm (yy, yy, ctx);
1159     }
1160
1161   if (!mpi_cmp_ui (point->z, 1))
1162     {
1163       mpi_set (x1, point->x);
1164       mpi_set (y1, yy);
1165     }
1166   else
1167     {
1168       gcry_mpi_t z2, z3;
1169
1170       z2 = mpi_alloc_like (ctx->p);
1171       z3 = mpi_alloc_like (ctx->p);
1172       ec_mulm (z2, point->z, point->z, ctx);
1173       ec_mulm (z3, point->z, z2, ctx);
1174       ec_invm (z2, z2, ctx);
1175       ec_mulm (x1, point->x, z2, ctx);
1176       ec_invm (z3, z3, ctx);
1177       ec_mulm (y1, yy, z3, ctx);
1178       mpi_free (z2);
1179       mpi_free (z3);
1180     }
1181   z1 = mpi_copy (mpi_const (MPI_C_ONE));
1182
1183   mpi_mul (h, k, mpi_const (MPI_C_THREE)); /* h = 3k */
1184   loops = mpi_get_nbits (h);
1185   if (loops < 2)
1186     {
1187       /* If SCALAR is zero, the above mpi_mul sets H to zero and thus
1188          LOOPs will be zero.  To avoid an underflow of I in the main
1189          loop we set LOOP to 2 and the result to (0,0,0).  */
1190       loops = 2;
1191       mpi_clear (result->x);
1192       mpi_clear (result->y);
1193       mpi_clear (result->z);
1194     }
1195   else
1196     {
1197       mpi_set (result->x, point->x);
1198       mpi_set (result->y, yy);
1199       mpi_set (result->z, point->z);
1200     }
1201   mpi_free (yy); yy = NULL;
1202
1203   p1.x = x1; x1 = NULL;
1204   p1.y = y1; y1 = NULL;
1205   p1.z = z1; z1 = NULL;
1206   point_init (&p2);
1207   point_init (&p1inv);
1208
1209   for (i=loops-2; i > 0; i--)
1210     {
1211       _gcry_mpi_ec_dup_point (result, result, ctx);
1212       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
1213         {
1214           point_set (&p2, result);
1215           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
1216         }
1217       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
1218         {
1219           point_set (&p2, result);
1220           /* Invert point: y = p - y mod p  */
1221           point_set (&p1inv, &p1);
1222           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
1223           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
1224         }
1225     }
1226
1227   point_free (&p1);
1228   point_free (&p2);
1229   point_free (&p1inv);
1230   mpi_free (h);
1231   mpi_free (k);
1232 }
1233
1234
1235 /* Return true if POINT is on the curve described by CTX.  */
1236 int
1237 _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
1238 {
1239   int res = 0;
1240   gcry_mpi_t x, y, w;
1241
1242   x = mpi_new (0);
1243   y = mpi_new (0);
1244   w = mpi_new (0);
1245
1246   if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1247     return 0;
1248
1249   switch (ctx->model)
1250     {
1251     case MPI_EC_WEIERSTRASS:
1252       {
1253         gcry_mpi_t xxx = mpi_new (0);
1254
1255         /* y^2 == x^3 + a·x + b */
1256         ec_pow2 (y, y, ctx);
1257
1258         ec_pow3 (xxx, x, ctx);
1259         ec_mulm (w, ctx->a, x, ctx);
1260         ec_addm (w, w, ctx->b, ctx);
1261         ec_addm (w, w, xxx, ctx);
1262
1263         if (!mpi_cmp (y, w))
1264           res = 1;
1265
1266         _gcry_mpi_release (xxx);
1267       }
1268       break;
1269     case MPI_EC_MONTGOMERY:
1270       log_fatal ("%s: %s not yet supported\n",
1271                  "_gcry_mpi_ec_curve_point", "Montgomery");
1272       break;
1273     case MPI_EC_EDWARDS:
1274       {
1275         /* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
1276         ec_pow2 (x, x, ctx);
1277         ec_pow2 (y, y, ctx);
1278         if (ctx->dialect == ECC_DIALECT_ED25519)
1279           {
1280             mpi_set (w, x);
1281             _gcry_mpi_neg (w, w);
1282           }
1283         else
1284           ec_mulm (w, ctx->a, x, ctx);
1285         ec_addm (w, w, y, ctx);
1286         ec_subm (w, w, mpi_const (MPI_C_ONE), ctx);
1287         ec_mulm (x, x, y, ctx);
1288         ec_mulm (x, x, ctx->b, ctx);
1289         ec_subm (w, w, x, ctx);
1290         if (!mpi_cmp_ui (w, 0))
1291           res = 1;
1292       }
1293       break;
1294     }
1295
1296   _gcry_mpi_release (w);
1297   _gcry_mpi_release (x);
1298   _gcry_mpi_release (y);
1299
1300   return res;
1301 }