Add GCRYMPI_FLAG_CONST and make use constants.
[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
31
32 #define point_init(a)  _gcry_mpi_point_init ((a))
33 #define point_free(a)  _gcry_mpi_point_free_parts ((a))
34
35
36 /* Object to represent a point in projective coordinates. */
37 /* Currently defined in mpi.h */
38
39 /* This context is used with all our EC functions. */
40 struct mpi_ec_ctx_s
41 {
42   /* Domain parameters.  */
43   gcry_mpi_t p;   /* Prime specifying the field GF(p).  */
44   gcry_mpi_t a;   /* First coefficient of the Weierstrass equation.  */
45
46   int a_is_pminus3;  /* True if A = P - 3. */
47
48   gcry_mpi_t two_inv_p;
49
50   /* Scratch variables.  */
51   gcry_mpi_t scratch[11];
52
53   /* Helper for fast reduction.  */
54 /*   int nist_nbits; /\* If this is a NIST curve, the number of bits.  *\/ */
55 /*   gcry_mpi_t s[10]; */
56 /*   gcry_mpi_t c; */
57
58 };
59
60
61 /* Create a new point option.  NBITS gives the size in bits of one
62    coordinate; it is only used to pre-allocate some resources and
63    might also be passed as 0 to use a default value.  */
64 mpi_point_t
65 gcry_mpi_point_new (unsigned int nbits)
66 {
67   mpi_point_t p;
68
69   (void)nbits;  /* Currently not used.  */
70
71   p = gcry_xmalloc (sizeof *p);
72   _gcry_mpi_point_init (p);
73   return p;
74 }
75
76
77 /* Release the point object P.  P may be NULL. */
78 void
79 gcry_mpi_point_release (mpi_point_t p)
80 {
81   if (p)
82     {
83       _gcry_mpi_point_free_parts (p);
84       gcry_free (p);
85     }
86 }
87
88
89 /* Initialize the fields of a point object.  gcry_mpi_point_free_parts
90    may be used to release the fields.  */
91 void
92 _gcry_mpi_point_init (mpi_point_t p)
93 {
94   p->x = mpi_new (0);
95   p->y = mpi_new (0);
96   p->z = mpi_new (0);
97 }
98
99
100 /* Release the parts of a point object. */
101 void
102 _gcry_mpi_point_free_parts (mpi_point_t p)
103 {
104   mpi_free (p->x); p->x = NULL;
105   mpi_free (p->y); p->y = NULL;
106   mpi_free (p->z); p->z = NULL;
107 }
108
109
110 /* Set the value from S into D.  */
111 static void
112 point_set (mpi_point_t d, mpi_point_t s)
113 {
114   mpi_set (d->x, s->x);
115   mpi_set (d->y, s->y);
116   mpi_set (d->z, s->z);
117 }
118
119 /* Set the projective coordinates from POINT into X, Y, and Z.  If a
120    coordinate is not required, X, Y, or Z may be passed as NULL.  */
121 void
122 gcry_mpi_point_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
123                     mpi_point_t point)
124 {
125   if (x)
126     mpi_set (x, point->x);
127   if (y)
128     mpi_set (y, point->y);
129   if (z)
130     mpi_set (z, point->z);
131 }
132
133
134 /* Set the projective coordinates from POINT into X, Y, and Z and
135    release POINT.  If a coordinate is not required, X, Y, or Z may be
136    passed as NULL.  */
137 void
138 gcry_mpi_point_snatch_get (gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z,
139                            mpi_point_t point)
140 {
141   mpi_snatch (x, point->x);
142   mpi_snatch (y, point->y);
143   mpi_snatch (z, point->z);
144   gcry_free (point);
145 }
146
147
148 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
149    coordinate is given as NULL, the value 0 is stored into point.  If
150    POINT is given as NULL a new point object is allocated.  Returns
151    POINT or the newly allocated point object. */
152 mpi_point_t
153 gcry_mpi_point_set (mpi_point_t point,
154                     gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
155 {
156   if (!point)
157     point = gcry_mpi_point_new (0);
158
159   if (x)
160     mpi_set (point->x, x);
161   else
162     mpi_clear (point->x);
163   if (y)
164     mpi_set (point->y, y);
165   else
166     mpi_clear (point->y);
167   if (z)
168     mpi_set (point->z, z);
169   else
170     mpi_clear (point->z);
171
172   return point;
173 }
174
175
176 /* Set the projective coordinates from X, Y, and Z into POINT.  If a
177    coordinate is given as NULL, the value 0 is stored into point.  If
178    POINT is given as NULL a new point object is allocated.  The
179    coordinates X, Y, and Z are released.  Returns POINT or the newly
180    allocated point object. */
181 mpi_point_t
182 gcry_mpi_point_snatch_set (mpi_point_t point,
183                            gcry_mpi_t x, gcry_mpi_t y, gcry_mpi_t z)
184 {
185   if (!point)
186     point = gcry_mpi_point_new (0);
187
188   if (x)
189     mpi_snatch (point->x, x);
190   else
191     mpi_clear (point->x);
192   if (y)
193     mpi_snatch (point->y, y);
194   else
195     mpi_clear (point->y);
196   if (z)
197     mpi_snatch (point->z, z);
198   else
199     mpi_clear (point->z);
200
201   return point;
202 }
203
204
205 static void
206 ec_addm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
207 {
208   mpi_addm (w, u, v, ctx->p);
209 }
210
211 static void
212 ec_subm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
213 {
214   mpi_subm (w, u, v, ctx->p);
215 }
216
217 static void
218 ec_mulm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
219 {
220 #if 0
221   /* NOTE: This code works only for limb sizes of 32 bit.  */
222   mpi_limb_t *wp, *sp;
223
224   if (ctx->nist_nbits == 192)
225     {
226       mpi_mul (w, u, v);
227       mpi_resize (w, 12);
228       wp = w->d;
229
230       sp = ctx->s[0]->d;
231       sp[0*2+0] = wp[0*2+0];
232       sp[0*2+1] = wp[0*2+1];
233       sp[1*2+0] = wp[1*2+0];
234       sp[1*2+1] = wp[1*2+1];
235       sp[2*2+0] = wp[2*2+0];
236       sp[2*2+1] = wp[2*2+1];
237
238       sp = ctx->s[1]->d;
239       sp[0*2+0] = wp[3*2+0];
240       sp[0*2+1] = wp[3*2+1];
241       sp[1*2+0] = wp[3*2+0];
242       sp[1*2+1] = wp[3*2+1];
243       sp[2*2+0] = 0;
244       sp[2*2+1] = 0;
245
246       sp = ctx->s[2]->d;
247       sp[0*2+0] = 0;
248       sp[0*2+1] = 0;
249       sp[1*2+0] = wp[4*2+0];
250       sp[1*2+1] = wp[4*2+1];
251       sp[2*2+0] = wp[4*2+0];
252       sp[2*2+1] = wp[4*2+1];
253
254       sp = ctx->s[3]->d;
255       sp[0*2+0] = wp[5*2+0];
256       sp[0*2+1] = wp[5*2+1];
257       sp[1*2+0] = wp[5*2+0];
258       sp[1*2+1] = wp[5*2+1];
259       sp[2*2+0] = wp[5*2+0];
260       sp[2*2+1] = wp[5*2+1];
261
262       ctx->s[0]->nlimbs = 6;
263       ctx->s[1]->nlimbs = 6;
264       ctx->s[2]->nlimbs = 6;
265       ctx->s[3]->nlimbs = 6;
266
267       mpi_add (ctx->c, ctx->s[0], ctx->s[1]);
268       mpi_add (ctx->c, ctx->c, ctx->s[2]);
269       mpi_add (ctx->c, ctx->c, ctx->s[3]);
270
271       while ( mpi_cmp (ctx->c, ctx->p ) >= 0 )
272         mpi_sub ( ctx->c, ctx->c, ctx->p );
273       mpi_set (w, ctx->c);
274     }
275   else if (ctx->nist_nbits == 384)
276     {
277       int i;
278       mpi_mul (w, u, v);
279       mpi_resize (w, 24);
280       wp = w->d;
281
282 #define NEXT(a) do { ctx->s[(a)]->nlimbs = 12; \
283                      sp = ctx->s[(a)]->d; \
284                      i = 0; } while (0)
285 #define X(a) do { sp[i++] = wp[(a)];} while (0)
286 #define X0(a) do { sp[i++] = 0; } while (0)
287       NEXT(0);
288       X(0);X(1);X(2);X(3);X(4);X(5);X(6);X(7);X(8);X(9);X(10);X(11);
289       NEXT(1);
290       X0();X0();X0();X0();X(21);X(22);X(23);X0();X0();X0();X0();X0();
291       NEXT(2);
292       X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);X(20);X(21);X(22);X(23);
293       NEXT(3);
294       X(21);X(22);X(23);X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);X(20);
295       NEXT(4);
296       X0();X(23);X0();X(20);X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);
297       NEXT(5);
298       X0();X0();X0();X0();X(20);X(21);X(22);X(23);X0();X0();X0();X0();
299       NEXT(6);
300       X(20);X0();X0();X(21);X(22);X(23);X0();X0();X0();X0();X0();X0();
301       NEXT(7);
302       X(23);X(12);X(13);X(14);X(15);X(16);X(17);X(18);X(19);X(20);X(21);X(22);
303       NEXT(8);
304       X0();X(20);X(21);X(22);X(23);X0();X0();X0();X0();X0();X0();X0();
305       NEXT(9);
306       X0();X0();X0();X(23);X(23);X0();X0();X0();X0();X0();X0();X0();
307 #undef X0
308 #undef X
309 #undef NEXT
310       mpi_add (ctx->c, ctx->s[0], ctx->s[1]);
311       mpi_add (ctx->c, ctx->c, ctx->s[1]);
312       mpi_add (ctx->c, ctx->c, ctx->s[2]);
313       mpi_add (ctx->c, ctx->c, ctx->s[3]);
314       mpi_add (ctx->c, ctx->c, ctx->s[4]);
315       mpi_add (ctx->c, ctx->c, ctx->s[5]);
316       mpi_add (ctx->c, ctx->c, ctx->s[6]);
317       mpi_sub (ctx->c, ctx->c, ctx->s[7]);
318       mpi_sub (ctx->c, ctx->c, ctx->s[8]);
319       mpi_sub (ctx->c, ctx->c, ctx->s[9]);
320
321       while ( mpi_cmp (ctx->c, ctx->p ) >= 0 )
322         mpi_sub ( ctx->c, ctx->c, ctx->p );
323       while ( ctx->c->sign )
324         mpi_add ( ctx->c, ctx->c, ctx->p );
325       mpi_set (w, ctx->c);
326     }
327   else
328 #endif /*0*/
329     mpi_mulm (w, u, v, ctx->p);
330 }
331
332 static void
333 ec_powm (gcry_mpi_t w, const gcry_mpi_t b, const gcry_mpi_t e,
334          mpi_ec_t ctx)
335 {
336   mpi_powm (w, b, e, ctx->p);
337 }
338
339 static void
340 ec_invm (gcry_mpi_t x, gcry_mpi_t a, mpi_ec_t ctx)
341 {
342   mpi_invm (x, a, ctx->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, gcry_mpi_t p, gcry_mpi_t a)
352 {
353   int i;
354   gcry_mpi_t tmp;
355
356   mpi_normalize (p);
357   mpi_normalize (a);
358
359   /* Fixme: Do we want to check some constraints? e.g.
360      a < p
361   */
362
363   ctx->p = mpi_copy (p);
364   ctx->a = mpi_copy (a);
365
366   tmp = mpi_alloc_like (ctx->p);
367   mpi_sub_ui (tmp, ctx->p, 3);
368   ctx->a_is_pminus3 = !mpi_cmp (ctx->a, tmp);
369   mpi_free (tmp);
370
371   ctx->two_inv_p = mpi_alloc (0);
372   ec_invm (ctx->two_inv_p, mpi_const (MPI_C_TWO), ctx);
373
374   /* Allocate scratch variables.  */
375   for (i=0; i< DIM(ctx->scratch); i++)
376     ctx->scratch[i] = mpi_alloc_like (ctx->p);
377
378   /* Prepare for fast reduction.  */
379   /* FIXME: need a test for NIST values.  However it does not gain us
380      any real advantage, for 384 bits it is actually slower than using
381      mpi_mulm.  */
382 /*   ctx->nist_nbits = mpi_get_nbits (ctx->p); */
383 /*   if (ctx->nist_nbits == 192) */
384 /*     { */
385 /*       for (i=0; i < 4; i++) */
386 /*         ctx->s[i] = mpi_new (192); */
387 /*       ctx->c    = mpi_new (192*2); */
388 /*     } */
389 /*   else if (ctx->nist_nbits == 384) */
390 /*     { */
391 /*       for (i=0; i < 10; i++) */
392 /*         ctx->s[i] = mpi_new (384); */
393 /*       ctx->c    = mpi_new (384*2); */
394 /*     } */
395 }
396
397
398 static void
399 ec_deinit (void *opaque)
400 {
401   mpi_ec_t ctx = opaque;
402   int i;
403
404   mpi_free (ctx->p);
405   mpi_free (ctx->a);
406
407   mpi_free (ctx->two_inv_p);
408
409   for (i=0; i< DIM(ctx->scratch); i++)
410     mpi_free (ctx->scratch[i]);
411
412 /*   if (ctx->nist_nbits == 192) */
413 /*     { */
414 /*       for (i=0; i < 4; i++) */
415 /*         mpi_free (ctx->s[i]); */
416 /*       mpi_free (ctx->c); */
417 /*     } */
418 /*   else if (ctx->nist_nbits == 384) */
419 /*     { */
420 /*       for (i=0; i < 10; i++) */
421 /*         mpi_free (ctx->s[i]); */
422 /*       mpi_free (ctx->c); */
423 /*     } */
424 }
425
426
427 /* This function returns a new context for elliptic curve based on the
428    field GF(p).  P is the prime specifying this field, A is the first
429    coefficient.  This function is only used within Libgcrypt and not
430    part of the public API.
431
432    This context needs to be released using _gcry_mpi_ec_free.  */
433 mpi_ec_t
434 _gcry_mpi_ec_p_internal_new (gcry_mpi_t p, gcry_mpi_t a)
435 {
436   mpi_ec_t ctx;
437
438   ctx = gcry_xcalloc (1, sizeof *ctx);
439   ec_p_init (ctx, p, a);
440
441   return ctx;
442 }
443
444
445 void
446 _gcry_mpi_ec_free (mpi_ec_t ctx)
447 {
448   if (ctx)
449     {
450       ec_deinit (ctx);
451       gcry_free (ctx);
452     }
453 }
454
455
456 /* This function returns a new context for elliptic curve operations
457    based on the field GF(p).  P is the prime specifying this field, A
458    is the first coefficient.  This function is part of the public API.
459    On error this function returns NULL and sets ERRNO.
460    The context needs to be released using gcry_ctx_release.  */
461 gcry_ctx_t
462 gcry_mpi_ec_p_new (gcry_mpi_t p, gcry_mpi_t a)
463 {
464   gcry_ctx_t ctx;
465   mpi_ec_t ec;
466
467   if (!p || !a || !mpi_cmp_ui (a, 0))
468     {
469       gpg_err_set_errno (EINVAL);
470       return NULL;
471     }
472
473   ctx = _gcry_ctx_alloc (CONTEXT_TYPE_EC, sizeof *ec, ec_deinit);
474   if (!ctx)
475     return NULL;
476   ec = _gcry_ctx_get_pointer (ctx, CONTEXT_TYPE_EC);
477   ec_p_init (ec, p, a);
478
479   return ctx;
480 }
481
482
483 /* Compute the affine coordinates from the projective coordinates in
484    POINT.  Set them into X and Y.  If one coordinate is not required,
485    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
486    on success or !0 if POINT is at infinity.  */
487 int
488 _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
489                          mpi_ec_t ctx)
490 {
491   gcry_mpi_t z1, z2, z3;
492
493   if (!mpi_cmp_ui (point->z, 0))
494     return -1;
495
496   z1 = mpi_new (0);
497   z2 = mpi_new (0);
498   ec_invm (z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
499   ec_mulm (z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
500
501   if (x)
502     ec_mulm (x, point->x, z2, ctx);
503
504   if (y)
505     {
506       z3 = mpi_new (0);
507       ec_mulm (z3, z2, z1, ctx);      /* z3 = z^(-3) mod p  */
508       ec_mulm (y, point->y, z3, ctx);
509       mpi_free (z3);
510     }
511
512   mpi_free (z2);
513   mpi_free (z1);
514   return 0;
515 }
516
517
518 \f
519 /*  RESULT = 2 * POINT  */
520 void
521 _gcry_mpi_ec_dup_point (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
522 {
523 #define x3 (result->x)
524 #define y3 (result->y)
525 #define z3 (result->z)
526 #define t1 (ctx->scratch[0])
527 #define t2 (ctx->scratch[1])
528 #define t3 (ctx->scratch[2])
529 #define l1 (ctx->scratch[3])
530 #define l2 (ctx->scratch[4])
531 #define l3 (ctx->scratch[5])
532
533   if (!mpi_cmp_ui (point->y, 0) || !mpi_cmp_ui (point->z, 0))
534     {
535       /* P_y == 0 || P_z == 0 => [1:1:0] */
536       mpi_set_ui (x3, 1);
537       mpi_set_ui (y3, 1);
538       mpi_set_ui (z3, 0);
539     }
540   else
541     {
542       if (ctx->a_is_pminus3)  /* Use the faster case.  */
543         {
544           /* L1 = 3(X - Z^2)(X + Z^2) */
545           /*                          T1: used for Z^2. */
546           /*                          T2: used for the right term.  */
547           ec_powm (t1, point->z, mpi_const (MPI_C_TWO), ctx);
548           ec_subm (l1, point->x, t1, ctx);
549           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
550           ec_addm (t2, point->x, t1, ctx);
551           ec_mulm (l1, l1, t2, ctx);
552         }
553       else /* Standard case. */
554         {
555           /* L1 = 3X^2 + aZ^4 */
556           /*                          T1: used for aZ^4. */
557           ec_powm (l1, point->x, mpi_const (MPI_C_TWO), ctx);
558           ec_mulm (l1, l1, mpi_const (MPI_C_THREE), ctx);
559           ec_powm (t1, point->z, mpi_const (MPI_C_FOUR), ctx);
560           ec_mulm (t1, t1, ctx->a, ctx);
561           ec_addm (l1, l1, t1, ctx);
562         }
563       /* Z3 = 2YZ */
564       ec_mulm (z3, point->y, point->z, ctx);
565       ec_mulm (z3, z3, mpi_const (MPI_C_TWO), ctx);
566
567       /* L2 = 4XY^2 */
568       /*                              T2: used for Y2; required later. */
569       ec_powm (t2, point->y, mpi_const (MPI_C_TWO), ctx);
570       ec_mulm (l2, t2, point->x, ctx);
571       ec_mulm (l2, l2, mpi_const (MPI_C_FOUR), ctx);
572
573       /* X3 = L1^2 - 2L2 */
574       /*                              T1: used for L2^2. */
575       ec_powm (x3, l1, mpi_const (MPI_C_TWO), ctx);
576       ec_mulm (t1, l2, mpi_const (MPI_C_TWO), ctx);
577       ec_subm (x3, x3, t1, ctx);
578
579       /* L3 = 8Y^4 */
580       /*                              T2: taken from above. */
581       ec_powm (t2, t2, mpi_const (MPI_C_TWO), ctx);
582       ec_mulm (l3, t2, mpi_const (MPI_C_EIGHT), ctx);
583
584       /* Y3 = L1(L2 - X3) - L3 */
585       ec_subm (y3, l2, x3, ctx);
586       ec_mulm (y3, y3, l1, ctx);
587       ec_subm (y3, y3, l3, ctx);
588     }
589
590 #undef x3
591 #undef y3
592 #undef z3
593 #undef t1
594 #undef t2
595 #undef t3
596 #undef l1
597 #undef l2
598 #undef l3
599 }
600
601
602
603 /* RESULT = P1 + P2 */
604 void
605 _gcry_mpi_ec_add_points (mpi_point_t result,
606                          mpi_point_t p1, mpi_point_t p2,
607                          mpi_ec_t ctx)
608 {
609 #define x1 (p1->x    )
610 #define y1 (p1->y    )
611 #define z1 (p1->z    )
612 #define x2 (p2->x    )
613 #define y2 (p2->y    )
614 #define z2 (p2->z    )
615 #define x3 (result->x)
616 #define y3 (result->y)
617 #define z3 (result->z)
618 #define l1 (ctx->scratch[0])
619 #define l2 (ctx->scratch[1])
620 #define l3 (ctx->scratch[2])
621 #define l4 (ctx->scratch[3])
622 #define l5 (ctx->scratch[4])
623 #define l6 (ctx->scratch[5])
624 #define l7 (ctx->scratch[6])
625 #define l8 (ctx->scratch[7])
626 #define l9 (ctx->scratch[8])
627 #define t1 (ctx->scratch[9])
628 #define t2 (ctx->scratch[10])
629
630   if ( (!mpi_cmp (x1, x2)) && (!mpi_cmp (y1, y2)) && (!mpi_cmp (z1, z2)) )
631     {
632       /* Same point; need to call the duplicate function.  */
633       _gcry_mpi_ec_dup_point (result, p1, ctx);
634     }
635   else if (!mpi_cmp_ui (z1, 0))
636     {
637       /* P1 is at infinity.  */
638       mpi_set (x3, p2->x);
639       mpi_set (y3, p2->y);
640       mpi_set (z3, p2->z);
641     }
642   else if (!mpi_cmp_ui (z2, 0))
643     {
644       /* P2 is at infinity.  */
645       mpi_set (x3, p1->x);
646       mpi_set (y3, p1->y);
647       mpi_set (z3, p1->z);
648     }
649   else
650     {
651       int z1_is_one = !mpi_cmp_ui (z1, 1);
652       int z2_is_one = !mpi_cmp_ui (z2, 1);
653
654       /* l1 = x1 z2^2  */
655       /* l2 = x2 z1^2  */
656       if (z2_is_one)
657         mpi_set (l1, x1);
658       else
659         {
660           ec_powm (l1, z2, mpi_const (MPI_C_TWO), ctx);
661           ec_mulm (l1, l1, x1, ctx);
662         }
663       if (z1_is_one)
664         mpi_set (l2, x1);
665       else
666         {
667           ec_powm (l2, z1, mpi_const (MPI_C_TWO), ctx);
668           ec_mulm (l2, l2, x2, ctx);
669         }
670       /* l3 = l1 - l2 */
671       ec_subm (l3, l1, l2, ctx);
672       /* l4 = y1 z2^3  */
673       ec_powm (l4, z2, mpi_const (MPI_C_THREE), ctx);
674       ec_mulm (l4, l4, y1, ctx);
675       /* l5 = y2 z1^3  */
676       ec_powm (l5, z1, mpi_const (MPI_C_THREE), ctx);
677       ec_mulm (l5, l5, y2, ctx);
678       /* l6 = l4 - l5  */
679       ec_subm (l6, l4, l5, ctx);
680
681       if (!mpi_cmp_ui (l3, 0))
682         {
683           if (!mpi_cmp_ui (l6, 0))
684             {
685               /* P1 and P2 are the same - use duplicate function.  */
686               _gcry_mpi_ec_dup_point (result, p1, ctx);
687             }
688           else
689             {
690               /* P1 is the inverse of P2.  */
691               mpi_set_ui (x3, 1);
692               mpi_set_ui (y3, 1);
693               mpi_set_ui (z3, 0);
694             }
695         }
696       else
697         {
698           /* l7 = l1 + l2  */
699           ec_addm (l7, l1, l2, ctx);
700           /* l8 = l4 + l5  */
701           ec_addm (l8, l4, l5, ctx);
702           /* z3 = z1 z2 l3  */
703           ec_mulm (z3, z1, z2, ctx);
704           ec_mulm (z3, z3, l3, ctx);
705           /* x3 = l6^2 - l7 l3^2  */
706           ec_powm (t1, l6, mpi_const (MPI_C_TWO), ctx);
707           ec_powm (t2, l3, mpi_const (MPI_C_TWO), ctx);
708           ec_mulm (t2, t2, l7, ctx);
709           ec_subm (x3, t1, t2, ctx);
710           /* l9 = l7 l3^2 - 2 x3  */
711           ec_mulm (t1, x3, mpi_const (MPI_C_TWO), ctx);
712           ec_subm (l9, t2, t1, ctx);
713           /* y3 = (l9 l6 - l8 l3^3)/2  */
714           ec_mulm (l9, l9, l6, ctx);
715           ec_powm (t1, l3, mpi_const (MPI_C_THREE), ctx); /* fixme: Use saved value*/
716           ec_mulm (t1, t1, l8, ctx);
717           ec_subm (y3, l9, t1, ctx);
718           ec_mulm (y3, y3, ctx->two_inv_p, ctx);
719         }
720     }
721
722 #undef x1
723 #undef y1
724 #undef z1
725 #undef x2
726 #undef y2
727 #undef z2
728 #undef x3
729 #undef y3
730 #undef z3
731 #undef l1
732 #undef l2
733 #undef l3
734 #undef l4
735 #undef l5
736 #undef l6
737 #undef l7
738 #undef l8
739 #undef l9
740 #undef t1
741 #undef t2
742 }
743
744
745
746 /* Scalar point multiplication - the main function for ECC.  If takes
747    an integer SCALAR and a POINT as well as the usual context CTX.
748    RESULT will be set to the resulting point. */
749 void
750 _gcry_mpi_ec_mul_point (mpi_point_t result,
751                         gcry_mpi_t scalar, mpi_point_t point,
752                         mpi_ec_t ctx)
753 {
754 #if 0
755   /* Simple left to right binary method.  GECC Algorithm 3.27 */
756   unsigned int nbits;
757   int i;
758
759   nbits = mpi_get_nbits (scalar);
760   mpi_set_ui (result->x, 1);
761   mpi_set_ui (result->y, 1);
762   mpi_set_ui (result->z, 0);
763
764   for (i=nbits-1; i >= 0; i--)
765     {
766       _gcry_mpi_ec_dup_point (result, result, ctx);
767       if (mpi_test_bit (scalar, i) == 1)
768         _gcry_mpi_ec_add_points (result, result, point, ctx);
769     }
770
771 #else
772   gcry_mpi_t x1, y1, z1, k, h, yy;
773   unsigned int i, loops;
774   mpi_point_struct p1, p2, p1inv;
775
776   x1 = mpi_alloc_like (ctx->p);
777   y1 = mpi_alloc_like (ctx->p);
778   h  = mpi_alloc_like (ctx->p);
779   k  = mpi_copy (scalar);
780   yy = mpi_copy (point->y);
781
782   if ( mpi_is_neg (k) )
783     {
784       k->sign = 0;
785       ec_invm (yy, yy, ctx);
786     }
787
788   if (!mpi_cmp_ui (point->z, 1))
789     {
790       mpi_set (x1, point->x);
791       mpi_set (y1, yy);
792     }
793   else
794     {
795       gcry_mpi_t z2, z3;
796
797       z2 = mpi_alloc_like (ctx->p);
798       z3 = mpi_alloc_like (ctx->p);
799       ec_mulm (z2, point->z, point->z, ctx);
800       ec_mulm (z3, point->z, z2, ctx);
801       ec_invm (z2, z2, ctx);
802       ec_mulm (x1, point->x, z2, ctx);
803       ec_invm (z3, z3, ctx);
804       ec_mulm (y1, yy, z3, ctx);
805       mpi_free (z2);
806       mpi_free (z3);
807     }
808   z1 = mpi_copy (mpi_const (MPI_C_ONE));
809
810   mpi_mul (h, k, mpi_const (MPI_C_THREE)); /* h = 3k */
811   loops = mpi_get_nbits (h);
812
813   mpi_set (result->x, point->x);
814   mpi_set (result->y, yy); mpi_free (yy); yy = NULL;
815   mpi_set (result->z, point->z);
816
817   p1.x = x1; x1 = NULL;
818   p1.y = y1; y1 = NULL;
819   p1.z = z1; z1 = NULL;
820   point_init (&p2);
821   point_init (&p1inv);
822
823   for (i=loops-2; i > 0; i--)
824     {
825       _gcry_mpi_ec_dup_point (result, result, ctx);
826       if (mpi_test_bit (h, i) == 1 && mpi_test_bit (k, i) == 0)
827         {
828           point_set (&p2, result);
829           _gcry_mpi_ec_add_points (result, &p2, &p1, ctx);
830         }
831       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
832         {
833           point_set (&p2, result);
834           /* Invert point: y = p - y mod p  */
835           point_set (&p1inv, &p1);
836           ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
837           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
838         }
839     }
840
841   point_free (&p1);
842   point_free (&p2);
843   point_free (&p1inv);
844   mpi_free (h);
845   mpi_free (k);
846 #endif
847 }