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