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