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