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