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