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