mpi: Add gcry_mpi_ec_sub.
[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       unsigned long sw;
1255
1256       /* Compute scalar point multiplication with Montgomery Ladder.
1257          Note that we don't use Y-coordinate in the points at all.
1258          RESULT->Y will be filled by zero.  */
1259
1260       nbits = mpi_get_nbits (scalar);
1261       point_init (&p1);
1262       point_init (&p2);
1263       point_init (&p1_);
1264       point_init (&p2_);
1265       mpi_set_ui (p1.x, 1);
1266       mpi_free (p2.x);
1267       p2.x  = mpi_copy (point->x);
1268       mpi_set_ui (p2.z, 1);
1269
1270       for (j=nbits-1; j >= 0; j--)
1271         {
1272           sw = mpi_test_bit (scalar, j);
1273           mpi_swap_cond (p1.x, p2.x, sw);
1274           mpi_swap_cond (p1.z, p2.z, sw);
1275           montgomery_ladder (&p1_, &p2_, &p1, &p2, point->x, ctx);
1276           mpi_swap_cond (p1_.x, p2_.x, sw);
1277           mpi_swap_cond (p1_.z, p2_.z, sw);
1278
1279           if (--j < 0)
1280             break;
1281
1282           sw = mpi_test_bit (scalar, j);
1283           mpi_swap_cond (p1_.x, p2_.x, sw);
1284           mpi_swap_cond (p1_.z, p2_.z, sw);
1285           montgomery_ladder (&p1, &p2, &p1_, &p2_, point->x, ctx);
1286           mpi_swap_cond (p1.x, p2.x, sw);
1287           mpi_swap_cond (p1.z, p2.z, sw);
1288         }
1289
1290       z1 = mpi_new (0);
1291       mpi_clear (result->y);
1292       sw = (nbits & 1);
1293       mpi_swap_cond (p1.x, p1_.x, sw);
1294       mpi_swap_cond (p1.z, p1_.z, sw);
1295
1296       if (p1.z->nlimbs == 0)
1297         {
1298           mpi_set_ui (result->x, 1);
1299           mpi_set_ui (result->z, 0);
1300         }
1301       else
1302         {
1303           ec_invm (z1, p1.z, ctx);
1304           ec_mulm (result->x, p1.x, z1, ctx);
1305           mpi_set_ui (result->z, 1);
1306         }
1307
1308       mpi_free (z1);
1309       point_free (&p1);
1310       point_free (&p2);
1311       point_free (&p1_);
1312       point_free (&p2_);
1313       return;
1314     }
1315
1316   x1 = mpi_alloc_like (ctx->p);
1317   y1 = mpi_alloc_like (ctx->p);
1318   h  = mpi_alloc_like (ctx->p);
1319   k  = mpi_copy (scalar);
1320   yy = mpi_copy (point->y);
1321
1322   if ( mpi_has_sign (k) )
1323     {
1324       k->sign = 0;
1325       ec_invm (yy, yy, ctx);
1326     }
1327
1328   if (!mpi_cmp_ui (point->z, 1))
1329     {
1330       mpi_set (x1, point->x);
1331       mpi_set (y1, yy);
1332     }
1333   else
1334     {
1335       gcry_mpi_t z2, z3;
1336
1337       z2 = mpi_alloc_like (ctx->p);
1338       z3 = mpi_alloc_like (ctx->p);
1339       ec_mulm (z2, point->z, point->z, ctx);
1340       ec_mulm (z3, point->z, z2, ctx);
1341       ec_invm (z2, z2, ctx);
1342       ec_mulm (x1, point->x, z2, ctx);
1343       ec_invm (z3, z3, ctx);
1344       ec_mulm (y1, yy, z3, ctx);
1345       mpi_free (z2);
1346       mpi_free (z3);
1347     }
1348   z1 = mpi_copy (mpi_const (MPI_C_ONE));
1349
1350   mpi_mul (h, k, mpi_const (MPI_C_THREE)); /* h = 3k */
1351   loops = mpi_get_nbits (h);
1352   if (loops < 2)
1353     {
1354       /* If SCALAR is zero, the above mpi_mul sets H to zero and thus
1355          LOOPs will be zero.  To avoid an underflow of I in the main
1356          loop we set LOOP to 2 and the result to (0,0,0).  */
1357       loops = 2;
1358       mpi_clear (result->x);
1359       mpi_clear (result->y);
1360       mpi_clear (result->z);
1361     }
1362   else
1363     {
1364       mpi_set (result->x, point->x);
1365       mpi_set (result->y, yy);
1366       mpi_set (result->z, point->z);
1367     }
1368   mpi_free (yy); yy = NULL;
1369
1370   p1.x = x1; x1 = NULL;
1371   p1.y = y1; y1 = NULL;
1372   p1.z = z1; z1 = NULL;
1373   point_init (&p2);
1374   point_init (&p1inv);
1375
1376   for (i=loops-2; i > 0; i--)
1377     {
1378       _gcry_mpi_ec_dup_point (result, result, ctx);
1379       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
1380         {
1381           point_set (&p2, result);
1382           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
1383         }
1384       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
1385         {
1386           point_set (&p2, result);
1387           /* Invert point: y = p - y mod p  */
1388           point_set (&p1inv, &p1);
1389           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
1390           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
1391         }
1392     }
1393
1394   point_free (&p1);
1395   point_free (&p2);
1396   point_free (&p1inv);
1397   mpi_free (h);
1398   mpi_free (k);
1399 }
1400
1401
1402 /* Return true if POINT is on the curve described by CTX.  */
1403 int
1404 _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
1405 {
1406   int res = 0;
1407   gcry_mpi_t x, y, w;
1408
1409   x = mpi_new (0);
1410   y = mpi_new (0);
1411   w = mpi_new (0);
1412
1413   switch (ctx->model)
1414     {
1415     case MPI_EC_WEIERSTRASS:
1416       {
1417         gcry_mpi_t xxx = mpi_new (0);
1418
1419         if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1420           return 0;
1421
1422         /* y^2 == x^3 + a·x + b */
1423         ec_pow2 (y, y, ctx);
1424
1425         ec_pow3 (xxx, x, ctx);
1426         ec_mulm (w, ctx->a, x, ctx);
1427         ec_addm (w, w, ctx->b, ctx);
1428         ec_addm (w, w, xxx, ctx);
1429
1430         if (!mpi_cmp (y, w))
1431           res = 1;
1432
1433         _gcry_mpi_release (xxx);
1434       }
1435       break;
1436     case MPI_EC_MONTGOMERY:
1437       {
1438 #define xx y
1439         /* With Montgomery curve, only X-coordinate is valid.  */
1440         if (_gcry_mpi_ec_get_affine (x, NULL, point, ctx))
1441           return 0;
1442
1443         /* The equation is: b * y^2 == x^3 + a · x^2 + x */
1444         /* We check if right hand is quadratic residue or not by
1445            Euler's criterion.  */
1446         /* CTX->A has (a-2)/4 and CTX->B has b^-1 */
1447         ec_mulm (w, ctx->a, mpi_const (MPI_C_FOUR), ctx);
1448         ec_addm (w, w, mpi_const (MPI_C_TWO), ctx);
1449         ec_mulm (w, w, x, ctx);
1450         ec_pow2 (xx, x, ctx);
1451         ec_addm (w, w, xx, ctx);
1452         ec_addm (w, w, mpi_const (MPI_C_ONE), ctx);
1453         ec_mulm (w, w, x, ctx);
1454         ec_mulm (w, w, ctx->b, ctx);
1455 #undef xx
1456         /* Compute Euler's criterion: w^(p-1)/2 */
1457 #define p_minus1 y
1458         ec_subm (p_minus1, ctx->p, mpi_const (MPI_C_ONE), ctx);
1459         mpi_rshift (p_minus1, p_minus1, 1);
1460         ec_powm (w, w, p_minus1, ctx);
1461
1462         res = mpi_cmp_ui (w, 1);
1463 #undef p_minus1
1464       }
1465       break;
1466     case MPI_EC_EDWARDS:
1467       {
1468         if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
1469           return 0;
1470
1471         /* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
1472         ec_pow2 (x, x, ctx);
1473         ec_pow2 (y, y, ctx);
1474         if (ctx->dialect == ECC_DIALECT_ED25519)
1475           {
1476             mpi_set (w, x);
1477             _gcry_mpi_neg (w, w);
1478           }
1479         else
1480           ec_mulm (w, ctx->a, x, ctx);
1481         ec_addm (w, w, y, ctx);
1482         ec_subm (w, w, mpi_const (MPI_C_ONE), ctx);
1483         ec_mulm (x, x, y, ctx);
1484         ec_mulm (x, x, ctx->b, ctx);
1485         ec_subm (w, w, x, ctx);
1486         if (!mpi_cmp_ui (w, 0))
1487           res = 1;
1488       }
1489       break;
1490     }
1491
1492   _gcry_mpi_release (w);
1493   _gcry_mpi_release (x);
1494   _gcry_mpi_release (y);
1495
1496   return res;
1497 }