xref: /openbmc/linux/lib/crypto/mpi/ec.c (revision 18afb028)
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->t.p_barrett = use_barrett > 0 ? mpi_barrett_init(ctx->p, 0) : NULL;
588 
589 	mpi_ec_get_reset(ctx);
590 
591 	if (model == MPI_EC_MONTGOMERY) {
592 		for (i = 0; i < DIM(bad_points_table); i++) {
593 			MPI p_candidate = mpi_scanval(bad_points_table[i][0]);
594 			int match_p = !mpi_cmp(ctx->p, p_candidate);
595 			int j;
596 
597 			mpi_free(p_candidate);
598 			if (!match_p)
599 				continue;
600 
601 			for (j = 0; i < DIM(ctx->t.scratch) && bad_points_table[i][j]; j++)
602 				ctx->t.scratch[j] = mpi_scanval(bad_points_table[i][j]);
603 		}
604 	} else {
605 		/* Allocate scratch variables.  */
606 		for (i = 0; i < DIM(ctx->t.scratch); i++)
607 			ctx->t.scratch[i] = mpi_alloc_like(ctx->p);
608 	}
609 
610 	ctx->addm = ec_addm;
611 	ctx->subm = ec_subm;
612 	ctx->mulm = ec_mulm;
613 	ctx->mul2 = ec_mul2;
614 	ctx->pow2 = ec_pow2;
615 
616 	for (i = 0; field_table[i].p; i++) {
617 		MPI f_p;
618 
619 		f_p = mpi_scanval(field_table[i].p);
620 		if (!f_p)
621 			break;
622 
623 		if (!mpi_cmp(p, f_p)) {
624 			ctx->addm = field_table[i].addm;
625 			ctx->subm = field_table[i].subm;
626 			ctx->mulm = field_table[i].mulm;
627 			ctx->mul2 = field_table[i].mul2;
628 			ctx->pow2 = field_table[i].pow2;
629 			mpi_free(f_p);
630 
631 			mpi_resize(ctx->a, ctx->p->nlimbs);
632 			ctx->a->nlimbs = ctx->p->nlimbs;
633 
634 			mpi_resize(ctx->b, ctx->p->nlimbs);
635 			ctx->b->nlimbs = ctx->p->nlimbs;
636 
637 			for (i = 0; i < DIM(ctx->t.scratch) && ctx->t.scratch[i]; i++)
638 				ctx->t.scratch[i]->nlimbs = ctx->p->nlimbs;
639 
640 			break;
641 		}
642 
643 		mpi_free(f_p);
644 	}
645 }
646 EXPORT_SYMBOL_GPL(mpi_ec_init);
647 
648 void mpi_ec_deinit(struct mpi_ec_ctx *ctx)
649 {
650 	int i;
651 
652 	mpi_barrett_free(ctx->t.p_barrett);
653 
654 	/* Domain parameter.  */
655 	mpi_free(ctx->p);
656 	mpi_free(ctx->a);
657 	mpi_free(ctx->b);
658 	mpi_point_release(ctx->G);
659 	mpi_free(ctx->n);
660 
661 	/* The key.  */
662 	mpi_point_release(ctx->Q);
663 	mpi_free(ctx->d);
664 
665 	/* Private data of ec.c.  */
666 	mpi_free(ctx->t.two_inv_p);
667 
668 	for (i = 0; i < DIM(ctx->t.scratch); i++)
669 		mpi_free(ctx->t.scratch[i]);
670 }
671 EXPORT_SYMBOL_GPL(mpi_ec_deinit);
672 
673 /* Compute the affine coordinates from the projective coordinates in
674  * POINT.  Set them into X and Y.  If one coordinate is not required,
675  * X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
676  * on success or !0 if POINT is at infinity.
677  */
678 int mpi_ec_get_affine(MPI x, MPI y, MPI_POINT point, struct mpi_ec_ctx *ctx)
679 {
680 	if (!mpi_cmp_ui(point->z, 0))
681 		return -1;
682 
683 	switch (ctx->model) {
684 	case MPI_EC_WEIERSTRASS: /* Using Jacobian coordinates.  */
685 		{
686 			MPI z1, z2, z3;
687 
688 			z1 = mpi_new(0);
689 			z2 = mpi_new(0);
690 			ec_invm(z1, point->z, ctx);  /* z1 = z^(-1) mod p  */
691 			ec_mulm(z2, z1, z1, ctx);    /* z2 = z^(-2) mod p  */
692 
693 			if (x)
694 				ec_mulm(x, point->x, z2, ctx);
695 
696 			if (y) {
697 				z3 = mpi_new(0);
698 				ec_mulm(z3, z2, z1, ctx);      /* z3 = z^(-3) mod p */
699 				ec_mulm(y, point->y, z3, ctx);
700 				mpi_free(z3);
701 			}
702 
703 			mpi_free(z2);
704 			mpi_free(z1);
705 		}
706 		return 0;
707 
708 	case MPI_EC_MONTGOMERY:
709 		{
710 			if (x)
711 				mpi_set(x, point->x);
712 
713 			if (y) {
714 				log_fatal("%s: Getting Y-coordinate on %s is not supported\n",
715 						"mpi_ec_get_affine", "Montgomery");
716 				return -1;
717 			}
718 		}
719 		return 0;
720 
721 	case MPI_EC_EDWARDS:
722 		{
723 			MPI z;
724 
725 			z = mpi_new(0);
726 			ec_invm(z, point->z, ctx);
727 
728 			mpi_resize(z, ctx->p->nlimbs);
729 			z->nlimbs = ctx->p->nlimbs;
730 
731 			if (x) {
732 				mpi_resize(x, ctx->p->nlimbs);
733 				x->nlimbs = ctx->p->nlimbs;
734 				ctx->mulm(x, point->x, z, ctx);
735 			}
736 			if (y) {
737 				mpi_resize(y, ctx->p->nlimbs);
738 				y->nlimbs = ctx->p->nlimbs;
739 				ctx->mulm(y, point->y, z, ctx);
740 			}
741 
742 			mpi_free(z);
743 		}
744 		return 0;
745 
746 	default:
747 		return -1;
748 	}
749 }
750 EXPORT_SYMBOL_GPL(mpi_ec_get_affine);
751 
752 /*  RESULT = 2 * POINT  (Weierstrass version). */
753 static void dup_point_weierstrass(MPI_POINT result,
754 		MPI_POINT point, struct mpi_ec_ctx *ctx)
755 {
756 #define x3 (result->x)
757 #define y3 (result->y)
758 #define z3 (result->z)
759 #define t1 (ctx->t.scratch[0])
760 #define t2 (ctx->t.scratch[1])
761 #define t3 (ctx->t.scratch[2])
762 #define l1 (ctx->t.scratch[3])
763 #define l2 (ctx->t.scratch[4])
764 #define l3 (ctx->t.scratch[5])
765 
766 	if (!mpi_cmp_ui(point->y, 0) || !mpi_cmp_ui(point->z, 0)) {
767 		/* P_y == 0 || P_z == 0 => [1:1:0] */
768 		mpi_set_ui(x3, 1);
769 		mpi_set_ui(y3, 1);
770 		mpi_set_ui(z3, 0);
771 	} else {
772 		if (ec_get_a_is_pminus3(ctx)) {
773 			/* Use the faster case.  */
774 			/* L1 = 3(X - Z^2)(X + Z^2) */
775 			/*                          T1: used for Z^2. */
776 			/*                          T2: used for the right term. */
777 			ec_pow2(t1, point->z, ctx);
778 			ec_subm(l1, point->x, t1, ctx);
779 			ec_mulm(l1, l1, mpi_const(MPI_C_THREE), ctx);
780 			ec_addm(t2, point->x, t1, ctx);
781 			ec_mulm(l1, l1, t2, ctx);
782 		} else {
783 			/* Standard case. */
784 			/* L1 = 3X^2 + aZ^4 */
785 			/*                          T1: used for aZ^4. */
786 			ec_pow2(l1, point->x, ctx);
787 			ec_mulm(l1, l1, mpi_const(MPI_C_THREE), ctx);
788 			ec_powm(t1, point->z, mpi_const(MPI_C_FOUR), ctx);
789 			ec_mulm(t1, t1, ctx->a, ctx);
790 			ec_addm(l1, l1, t1, ctx);
791 		}
792 		/* Z3 = 2YZ */
793 		ec_mulm(z3, point->y, point->z, ctx);
794 		ec_mul2(z3, z3, ctx);
795 
796 		/* L2 = 4XY^2 */
797 		/*                              T2: used for Y2; required later. */
798 		ec_pow2(t2, point->y, ctx);
799 		ec_mulm(l2, t2, point->x, ctx);
800 		ec_mulm(l2, l2, mpi_const(MPI_C_FOUR), ctx);
801 
802 		/* X3 = L1^2 - 2L2 */
803 		/*                              T1: used for L2^2. */
804 		ec_pow2(x3, l1, ctx);
805 		ec_mul2(t1, l2, ctx);
806 		ec_subm(x3, x3, t1, ctx);
807 
808 		/* L3 = 8Y^4 */
809 		/*                              T2: taken from above. */
810 		ec_pow2(t2, t2, ctx);
811 		ec_mulm(l3, t2, mpi_const(MPI_C_EIGHT), ctx);
812 
813 		/* Y3 = L1(L2 - X3) - L3 */
814 		ec_subm(y3, l2, x3, ctx);
815 		ec_mulm(y3, y3, l1, ctx);
816 		ec_subm(y3, y3, l3, ctx);
817 	}
818 
819 #undef x3
820 #undef y3
821 #undef z3
822 #undef t1
823 #undef t2
824 #undef t3
825 #undef l1
826 #undef l2
827 #undef l3
828 }
829 
830 /*  RESULT = 2 * POINT  (Montgomery version). */
831 static void dup_point_montgomery(MPI_POINT result,
832 				MPI_POINT point, struct mpi_ec_ctx *ctx)
833 {
834 	(void)result;
835 	(void)point;
836 	(void)ctx;
837 	log_fatal("%s: %s not yet supported\n",
838 			"mpi_ec_dup_point", "Montgomery");
839 }
840 
841 /*  RESULT = 2 * POINT  (Twisted Edwards version). */
842 static void dup_point_edwards(MPI_POINT result,
843 		MPI_POINT point, struct mpi_ec_ctx *ctx)
844 {
845 #define X1 (point->x)
846 #define Y1 (point->y)
847 #define Z1 (point->z)
848 #define X3 (result->x)
849 #define Y3 (result->y)
850 #define Z3 (result->z)
851 #define B (ctx->t.scratch[0])
852 #define C (ctx->t.scratch[1])
853 #define D (ctx->t.scratch[2])
854 #define E (ctx->t.scratch[3])
855 #define F (ctx->t.scratch[4])
856 #define H (ctx->t.scratch[5])
857 #define J (ctx->t.scratch[6])
858 
859 	/* Compute: (X_3 : Y_3 : Z_3) = 2( X_1 : Y_1 : Z_1 ) */
860 
861 	/* B = (X_1 + Y_1)^2  */
862 	ctx->addm(B, X1, Y1, ctx);
863 	ctx->pow2(B, B, ctx);
864 
865 	/* C = X_1^2 */
866 	/* D = Y_1^2 */
867 	ctx->pow2(C, X1, ctx);
868 	ctx->pow2(D, Y1, ctx);
869 
870 	/* E = aC */
871 	if (ctx->dialect == ECC_DIALECT_ED25519)
872 		ctx->subm(E, ctx->p, C, ctx);
873 	else
874 		ctx->mulm(E, ctx->a, C, ctx);
875 
876 	/* F = E + D */
877 	ctx->addm(F, E, D, ctx);
878 
879 	/* H = Z_1^2 */
880 	ctx->pow2(H, Z1, ctx);
881 
882 	/* J = F - 2H */
883 	ctx->mul2(J, H, ctx);
884 	ctx->subm(J, F, J, ctx);
885 
886 	/* X_3 = (B - C - D) · J */
887 	ctx->subm(X3, B, C, ctx);
888 	ctx->subm(X3, X3, D, ctx);
889 	ctx->mulm(X3, X3, J, ctx);
890 
891 	/* Y_3 = F · (E - D) */
892 	ctx->subm(Y3, E, D, ctx);
893 	ctx->mulm(Y3, Y3, F, ctx);
894 
895 	/* Z_3 = F · J */
896 	ctx->mulm(Z3, F, J, ctx);
897 
898 #undef X1
899 #undef Y1
900 #undef Z1
901 #undef X3
902 #undef Y3
903 #undef Z3
904 #undef B
905 #undef C
906 #undef D
907 #undef E
908 #undef F
909 #undef H
910 #undef J
911 }
912 
913 /*  RESULT = 2 * POINT  */
914 static void
915 mpi_ec_dup_point(MPI_POINT result, MPI_POINT point, struct mpi_ec_ctx *ctx)
916 {
917 	switch (ctx->model) {
918 	case MPI_EC_WEIERSTRASS:
919 		dup_point_weierstrass(result, point, ctx);
920 		break;
921 	case MPI_EC_MONTGOMERY:
922 		dup_point_montgomery(result, point, ctx);
923 		break;
924 	case MPI_EC_EDWARDS:
925 		dup_point_edwards(result, point, ctx);
926 		break;
927 	}
928 }
929 
930 /* RESULT = P1 + P2  (Weierstrass version).*/
931 static void add_points_weierstrass(MPI_POINT result,
932 		MPI_POINT p1, MPI_POINT p2,
933 		struct mpi_ec_ctx *ctx)
934 {
935 #define x1 (p1->x)
936 #define y1 (p1->y)
937 #define z1 (p1->z)
938 #define x2 (p2->x)
939 #define y2 (p2->y)
940 #define z2 (p2->z)
941 #define x3 (result->x)
942 #define y3 (result->y)
943 #define z3 (result->z)
944 #define l1 (ctx->t.scratch[0])
945 #define l2 (ctx->t.scratch[1])
946 #define l3 (ctx->t.scratch[2])
947 #define l4 (ctx->t.scratch[3])
948 #define l5 (ctx->t.scratch[4])
949 #define l6 (ctx->t.scratch[5])
950 #define l7 (ctx->t.scratch[6])
951 #define l8 (ctx->t.scratch[7])
952 #define l9 (ctx->t.scratch[8])
953 #define t1 (ctx->t.scratch[9])
954 #define t2 (ctx->t.scratch[10])
955 
956 	if ((!mpi_cmp(x1, x2)) && (!mpi_cmp(y1, y2)) && (!mpi_cmp(z1, z2))) {
957 		/* Same point; need to call the duplicate function.  */
958 		mpi_ec_dup_point(result, p1, ctx);
959 	} else if (!mpi_cmp_ui(z1, 0)) {
960 		/* P1 is at infinity.  */
961 		mpi_set(x3, p2->x);
962 		mpi_set(y3, p2->y);
963 		mpi_set(z3, p2->z);
964 	} else if (!mpi_cmp_ui(z2, 0)) {
965 		/* P2 is at infinity.  */
966 		mpi_set(x3, p1->x);
967 		mpi_set(y3, p1->y);
968 		mpi_set(z3, p1->z);
969 	} else {
970 		int z1_is_one = !mpi_cmp_ui(z1, 1);
971 		int z2_is_one = !mpi_cmp_ui(z2, 1);
972 
973 		/* l1 = x1 z2^2  */
974 		/* l2 = x2 z1^2  */
975 		if (z2_is_one)
976 			mpi_set(l1, x1);
977 		else {
978 			ec_pow2(l1, z2, ctx);
979 			ec_mulm(l1, l1, x1, ctx);
980 		}
981 		if (z1_is_one)
982 			mpi_set(l2, x2);
983 		else {
984 			ec_pow2(l2, z1, ctx);
985 			ec_mulm(l2, l2, x2, ctx);
986 		}
987 		/* l3 = l1 - l2 */
988 		ec_subm(l3, l1, l2, ctx);
989 		/* l4 = y1 z2^3  */
990 		ec_powm(l4, z2, mpi_const(MPI_C_THREE), ctx);
991 		ec_mulm(l4, l4, y1, ctx);
992 		/* l5 = y2 z1^3  */
993 		ec_powm(l5, z1, mpi_const(MPI_C_THREE), ctx);
994 		ec_mulm(l5, l5, y2, ctx);
995 		/* l6 = l4 - l5  */
996 		ec_subm(l6, l4, l5, ctx);
997 
998 		if (!mpi_cmp_ui(l3, 0)) {
999 			if (!mpi_cmp_ui(l6, 0)) {
1000 				/* P1 and P2 are the same - use duplicate function. */
1001 				mpi_ec_dup_point(result, p1, ctx);
1002 			} else {
1003 				/* P1 is the inverse of P2.  */
1004 				mpi_set_ui(x3, 1);
1005 				mpi_set_ui(y3, 1);
1006 				mpi_set_ui(z3, 0);
1007 			}
1008 		} else {
1009 			/* l7 = l1 + l2  */
1010 			ec_addm(l7, l1, l2, ctx);
1011 			/* l8 = l4 + l5  */
1012 			ec_addm(l8, l4, l5, ctx);
1013 			/* z3 = z1 z2 l3  */
1014 			ec_mulm(z3, z1, z2, ctx);
1015 			ec_mulm(z3, z3, l3, ctx);
1016 			/* x3 = l6^2 - l7 l3^2  */
1017 			ec_pow2(t1, l6, ctx);
1018 			ec_pow2(t2, l3, ctx);
1019 			ec_mulm(t2, t2, l7, ctx);
1020 			ec_subm(x3, t1, t2, ctx);
1021 			/* l9 = l7 l3^2 - 2 x3  */
1022 			ec_mul2(t1, x3, ctx);
1023 			ec_subm(l9, t2, t1, ctx);
1024 			/* y3 = (l9 l6 - l8 l3^3)/2  */
1025 			ec_mulm(l9, l9, l6, ctx);
1026 			ec_powm(t1, l3, mpi_const(MPI_C_THREE), ctx); /* fixme: Use saved value*/
1027 			ec_mulm(t1, t1, l8, ctx);
1028 			ec_subm(y3, l9, t1, ctx);
1029 			ec_mulm(y3, y3, ec_get_two_inv_p(ctx), ctx);
1030 		}
1031 	}
1032 
1033 #undef x1
1034 #undef y1
1035 #undef z1
1036 #undef x2
1037 #undef y2
1038 #undef z2
1039 #undef x3
1040 #undef y3
1041 #undef z3
1042 #undef l1
1043 #undef l2
1044 #undef l3
1045 #undef l4
1046 #undef l5
1047 #undef l6
1048 #undef l7
1049 #undef l8
1050 #undef l9
1051 #undef t1
1052 #undef t2
1053 }
1054 
1055 /* RESULT = P1 + P2  (Montgomery version).*/
1056 static void add_points_montgomery(MPI_POINT result,
1057 		MPI_POINT p1, MPI_POINT p2,
1058 		struct mpi_ec_ctx *ctx)
1059 {
1060 	(void)result;
1061 	(void)p1;
1062 	(void)p2;
1063 	(void)ctx;
1064 	log_fatal("%s: %s not yet supported\n",
1065 			"mpi_ec_add_points", "Montgomery");
1066 }
1067 
1068 /* RESULT = P1 + P2  (Twisted Edwards version).*/
1069 static void add_points_edwards(MPI_POINT result,
1070 		MPI_POINT p1, MPI_POINT p2,
1071 		struct mpi_ec_ctx *ctx)
1072 {
1073 #define X1 (p1->x)
1074 #define Y1 (p1->y)
1075 #define Z1 (p1->z)
1076 #define X2 (p2->x)
1077 #define Y2 (p2->y)
1078 #define Z2 (p2->z)
1079 #define X3 (result->x)
1080 #define Y3 (result->y)
1081 #define Z3 (result->z)
1082 #define A (ctx->t.scratch[0])
1083 #define B (ctx->t.scratch[1])
1084 #define C (ctx->t.scratch[2])
1085 #define D (ctx->t.scratch[3])
1086 #define E (ctx->t.scratch[4])
1087 #define F (ctx->t.scratch[5])
1088 #define G (ctx->t.scratch[6])
1089 #define tmp (ctx->t.scratch[7])
1090 
1091 	point_resize(result, ctx);
1092 
1093 	/* Compute: (X_3 : Y_3 : Z_3) = (X_1 : Y_1 : Z_1) + (X_2 : Y_2 : Z_3) */
1094 
1095 	/* A = Z1 · Z2 */
1096 	ctx->mulm(A, Z1, Z2, ctx);
1097 
1098 	/* B = A^2 */
1099 	ctx->pow2(B, A, ctx);
1100 
1101 	/* C = X1 · X2 */
1102 	ctx->mulm(C, X1, X2, ctx);
1103 
1104 	/* D = Y1 · Y2 */
1105 	ctx->mulm(D, Y1, Y2, ctx);
1106 
1107 	/* E = d · C · D */
1108 	ctx->mulm(E, ctx->b, C, ctx);
1109 	ctx->mulm(E, E, D, ctx);
1110 
1111 	/* F = B - E */
1112 	ctx->subm(F, B, E, ctx);
1113 
1114 	/* G = B + E */
1115 	ctx->addm(G, B, E, ctx);
1116 
1117 	/* X_3 = A · F · ((X_1 + Y_1) · (X_2 + Y_2) - C - D) */
1118 	ctx->addm(tmp, X1, Y1, ctx);
1119 	ctx->addm(X3, X2, Y2, ctx);
1120 	ctx->mulm(X3, X3, tmp, ctx);
1121 	ctx->subm(X3, X3, C, ctx);
1122 	ctx->subm(X3, X3, D, ctx);
1123 	ctx->mulm(X3, X3, F, ctx);
1124 	ctx->mulm(X3, X3, A, ctx);
1125 
1126 	/* Y_3 = A · G · (D - aC) */
1127 	if (ctx->dialect == ECC_DIALECT_ED25519) {
1128 		ctx->addm(Y3, D, C, ctx);
1129 	} else {
1130 		ctx->mulm(Y3, ctx->a, C, ctx);
1131 		ctx->subm(Y3, D, Y3, ctx);
1132 	}
1133 	ctx->mulm(Y3, Y3, G, ctx);
1134 	ctx->mulm(Y3, Y3, A, ctx);
1135 
1136 	/* Z_3 = F · G */
1137 	ctx->mulm(Z3, F, G, ctx);
1138 
1139 
1140 #undef X1
1141 #undef Y1
1142 #undef Z1
1143 #undef X2
1144 #undef Y2
1145 #undef Z2
1146 #undef X3
1147 #undef Y3
1148 #undef Z3
1149 #undef A
1150 #undef B
1151 #undef C
1152 #undef D
1153 #undef E
1154 #undef F
1155 #undef G
1156 #undef tmp
1157 }
1158 
1159 /* Compute a step of Montgomery Ladder (only use X and Z in the point).
1160  * Inputs:  P1, P2, and x-coordinate of DIF = P1 - P1.
1161  * Outputs: PRD = 2 * P1 and  SUM = P1 + P2.
1162  */
1163 static void montgomery_ladder(MPI_POINT prd, MPI_POINT sum,
1164 		MPI_POINT p1, MPI_POINT p2, MPI dif_x,
1165 		struct mpi_ec_ctx *ctx)
1166 {
1167 	ctx->addm(sum->x, p2->x, p2->z, ctx);
1168 	ctx->subm(p2->z, p2->x, p2->z, ctx);
1169 	ctx->addm(prd->x, p1->x, p1->z, ctx);
1170 	ctx->subm(p1->z, p1->x, p1->z, ctx);
1171 	ctx->mulm(p2->x, p1->z, sum->x, ctx);
1172 	ctx->mulm(p2->z, prd->x, p2->z, ctx);
1173 	ctx->pow2(p1->x, prd->x, ctx);
1174 	ctx->pow2(p1->z, p1->z, ctx);
1175 	ctx->addm(sum->x, p2->x, p2->z, ctx);
1176 	ctx->subm(p2->z, p2->x, p2->z, ctx);
1177 	ctx->mulm(prd->x, p1->x, p1->z, ctx);
1178 	ctx->subm(p1->z, p1->x, p1->z, ctx);
1179 	ctx->pow2(sum->x, sum->x, ctx);
1180 	ctx->pow2(sum->z, p2->z, ctx);
1181 	ctx->mulm(prd->z, p1->z, ctx->a, ctx); /* CTX->A: (a-2)/4 */
1182 	ctx->mulm(sum->z, sum->z, dif_x, ctx);
1183 	ctx->addm(prd->z, p1->x, prd->z, ctx);
1184 	ctx->mulm(prd->z, prd->z, p1->z, ctx);
1185 }
1186 
1187 /* RESULT = P1 + P2 */
1188 void mpi_ec_add_points(MPI_POINT result,
1189 		MPI_POINT p1, MPI_POINT p2,
1190 		struct mpi_ec_ctx *ctx)
1191 {
1192 	switch (ctx->model) {
1193 	case MPI_EC_WEIERSTRASS:
1194 		add_points_weierstrass(result, p1, p2, ctx);
1195 		break;
1196 	case MPI_EC_MONTGOMERY:
1197 		add_points_montgomery(result, p1, p2, ctx);
1198 		break;
1199 	case MPI_EC_EDWARDS:
1200 		add_points_edwards(result, p1, p2, ctx);
1201 		break;
1202 	}
1203 }
1204 EXPORT_SYMBOL_GPL(mpi_ec_add_points);
1205 
1206 /* Scalar point multiplication - the main function for ECC.  If takes
1207  * an integer SCALAR and a POINT as well as the usual context CTX.
1208  * RESULT will be set to the resulting point.
1209  */
1210 void mpi_ec_mul_point(MPI_POINT result,
1211 			MPI scalar, MPI_POINT point,
1212 			struct mpi_ec_ctx *ctx)
1213 {
1214 	MPI x1, y1, z1, k, h, yy;
1215 	unsigned int i, loops;
1216 	struct gcry_mpi_point p1, p2, p1inv;
1217 
1218 	if (ctx->model == MPI_EC_EDWARDS) {
1219 		/* Simple left to right binary method.  Algorithm 3.27 from
1220 		 * {author={Hankerson, Darrel and Menezes, Alfred J. and Vanstone, Scott},
1221 		 *  title = {Guide to Elliptic Curve Cryptography},
1222 		 *  year = {2003}, isbn = {038795273X},
1223 		 *  url = {http://www.cacr.math.uwaterloo.ca/ecc/},
1224 		 *  publisher = {Springer-Verlag New York, Inc.}}
1225 		 */
1226 		unsigned int nbits;
1227 		int j;
1228 
1229 		if (mpi_cmp(scalar, ctx->p) >= 0)
1230 			nbits = mpi_get_nbits(scalar);
1231 		else
1232 			nbits = mpi_get_nbits(ctx->p);
1233 
1234 		mpi_set_ui(result->x, 0);
1235 		mpi_set_ui(result->y, 1);
1236 		mpi_set_ui(result->z, 1);
1237 		point_resize(point, ctx);
1238 
1239 		point_resize(result, ctx);
1240 		point_resize(point, ctx);
1241 
1242 		for (j = nbits-1; j >= 0; j--) {
1243 			mpi_ec_dup_point(result, result, ctx);
1244 			if (mpi_test_bit(scalar, j))
1245 				mpi_ec_add_points(result, result, point, ctx);
1246 		}
1247 		return;
1248 	} else if (ctx->model == MPI_EC_MONTGOMERY) {
1249 		unsigned int nbits;
1250 		int j;
1251 		struct gcry_mpi_point p1_, p2_;
1252 		MPI_POINT q1, q2, prd, sum;
1253 		unsigned long sw;
1254 		mpi_size_t rsize;
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 
1261 		nbits = mpi_get_nbits(scalar);
1262 		point_init(&p1);
1263 		point_init(&p2);
1264 		point_init(&p1_);
1265 		point_init(&p2_);
1266 		mpi_set_ui(p1.x, 1);
1267 		mpi_free(p2.x);
1268 		p2.x = mpi_copy(point->x);
1269 		mpi_set_ui(p2.z, 1);
1270 
1271 		point_resize(&p1, ctx);
1272 		point_resize(&p2, ctx);
1273 		point_resize(&p1_, ctx);
1274 		point_resize(&p2_, ctx);
1275 
1276 		mpi_resize(point->x, ctx->p->nlimbs);
1277 		point->x->nlimbs = ctx->p->nlimbs;
1278 
1279 		q1 = &p1;
1280 		q2 = &p2;
1281 		prd = &p1_;
1282 		sum = &p2_;
1283 
1284 		for (j = nbits-1; j >= 0; j--) {
1285 			MPI_POINT t;
1286 
1287 			sw = mpi_test_bit(scalar, j);
1288 			point_swap_cond(q1, q2, sw, ctx);
1289 			montgomery_ladder(prd, sum, q1, q2, point->x, ctx);
1290 			point_swap_cond(prd, sum, sw, ctx);
1291 			t = q1;  q1 = prd;  prd = t;
1292 			t = q2;  q2 = sum;  sum = t;
1293 		}
1294 
1295 		mpi_clear(result->y);
1296 		sw = (nbits & 1);
1297 		point_swap_cond(&p1, &p1_, sw, ctx);
1298 
1299 		rsize = p1.z->nlimbs;
1300 		MPN_NORMALIZE(p1.z->d, rsize);
1301 		if (rsize == 0) {
1302 			mpi_set_ui(result->x, 1);
1303 			mpi_set_ui(result->z, 0);
1304 		} else {
1305 			z1 = mpi_new(0);
1306 			ec_invm(z1, p1.z, ctx);
1307 			ec_mulm(result->x, p1.x, z1, ctx);
1308 			mpi_set_ui(result->z, 1);
1309 			mpi_free(z1);
1310 		}
1311 
1312 		point_free(&p1);
1313 		point_free(&p2);
1314 		point_free(&p1_);
1315 		point_free(&p2_);
1316 		return;
1317 	}
1318 
1319 	x1 = mpi_alloc_like(ctx->p);
1320 	y1 = mpi_alloc_like(ctx->p);
1321 	h  = mpi_alloc_like(ctx->p);
1322 	k  = mpi_copy(scalar);
1323 	yy = mpi_copy(point->y);
1324 
1325 	if (mpi_has_sign(k)) {
1326 		k->sign = 0;
1327 		ec_invm(yy, yy, ctx);
1328 	}
1329 
1330 	if (!mpi_cmp_ui(point->z, 1)) {
1331 		mpi_set(x1, point->x);
1332 		mpi_set(y1, yy);
1333 	} else {
1334 		MPI z2, z3;
1335 
1336 		z2 = mpi_alloc_like(ctx->p);
1337 		z3 = mpi_alloc_like(ctx->p);
1338 		ec_mulm(z2, point->z, point->z, ctx);
1339 		ec_mulm(z3, point->z, z2, ctx);
1340 		ec_invm(z2, z2, ctx);
1341 		ec_mulm(x1, point->x, z2, ctx);
1342 		ec_invm(z3, z3, ctx);
1343 		ec_mulm(y1, yy, z3, ctx);
1344 		mpi_free(z2);
1345 		mpi_free(z3);
1346 	}
1347 	z1 = mpi_copy(mpi_const(MPI_C_ONE));
1348 
1349 	mpi_mul(h, k, mpi_const(MPI_C_THREE)); /* h = 3k */
1350 	loops = mpi_get_nbits(h);
1351 	if (loops < 2) {
1352 		/* If SCALAR is zero, the above mpi_mul sets H to zero and thus
1353 		 * LOOPs will be zero.  To avoid an underflow of I in the main
1354 		 * loop we set LOOP to 2 and the result to (0,0,0).
1355 		 */
1356 		loops = 2;
1357 		mpi_clear(result->x);
1358 		mpi_clear(result->y);
1359 		mpi_clear(result->z);
1360 	} else {
1361 		mpi_set(result->x, point->x);
1362 		mpi_set(result->y, yy);
1363 		mpi_set(result->z, point->z);
1364 	}
1365 	mpi_free(yy); yy = NULL;
1366 
1367 	p1.x = x1; x1 = NULL;
1368 	p1.y = y1; y1 = NULL;
1369 	p1.z = z1; z1 = NULL;
1370 	point_init(&p2);
1371 	point_init(&p1inv);
1372 
1373 	/* Invert point: y = p - y mod p  */
1374 	point_set(&p1inv, &p1);
1375 	ec_subm(p1inv.y, ctx->p, p1inv.y, ctx);
1376 
1377 	for (i = loops-2; i > 0; i--) {
1378 		mpi_ec_dup_point(result, result, ctx);
1379 		if (mpi_test_bit(h, i) == 1 && mpi_test_bit(k, i) == 0) {
1380 			point_set(&p2, result);
1381 			mpi_ec_add_points(result, &p2, &p1, ctx);
1382 		}
1383 		if (mpi_test_bit(h, i) == 0 && mpi_test_bit(k, i) == 1) {
1384 			point_set(&p2, result);
1385 			mpi_ec_add_points(result, &p2, &p1inv, ctx);
1386 		}
1387 	}
1388 
1389 	point_free(&p1);
1390 	point_free(&p2);
1391 	point_free(&p1inv);
1392 	mpi_free(h);
1393 	mpi_free(k);
1394 }
1395 EXPORT_SYMBOL_GPL(mpi_ec_mul_point);
1396 
1397 /* Return true if POINT is on the curve described by CTX.  */
1398 int mpi_ec_curve_point(MPI_POINT point, struct mpi_ec_ctx *ctx)
1399 {
1400 	int res = 0;
1401 	MPI x, y, w;
1402 
1403 	x = mpi_new(0);
1404 	y = mpi_new(0);
1405 	w = mpi_new(0);
1406 
1407 	/* Check that the point is in range.  This needs to be done here and
1408 	 * not after conversion to affine coordinates.
1409 	 */
1410 	if (mpi_cmpabs(point->x, ctx->p) >= 0)
1411 		goto leave;
1412 	if (mpi_cmpabs(point->y, ctx->p) >= 0)
1413 		goto leave;
1414 	if (mpi_cmpabs(point->z, ctx->p) >= 0)
1415 		goto leave;
1416 
1417 	switch (ctx->model) {
1418 	case MPI_EC_WEIERSTRASS:
1419 		{
1420 			MPI xxx;
1421 
1422 			if (mpi_ec_get_affine(x, y, point, ctx))
1423 				goto leave;
1424 
1425 			xxx = mpi_new(0);
1426 
1427 			/* y^2 == x^3 + a·x + b */
1428 			ec_pow2(y, y, ctx);
1429 
1430 			ec_pow3(xxx, x, ctx);
1431 			ec_mulm(w, ctx->a, x, ctx);
1432 			ec_addm(w, w, ctx->b, ctx);
1433 			ec_addm(w, w, xxx, ctx);
1434 
1435 			if (!mpi_cmp(y, w))
1436 				res = 1;
1437 
1438 			mpi_free(xxx);
1439 		}
1440 		break;
1441 
1442 	case MPI_EC_MONTGOMERY:
1443 		{
1444 #define xx y
1445 			/* With Montgomery curve, only X-coordinate is valid. */
1446 			if (mpi_ec_get_affine(x, NULL, point, ctx))
1447 				goto leave;
1448 
1449 			/* The equation is: b * y^2 == x^3 + a · x^2 + x */
1450 			/* We check if right hand is quadratic residue or not by
1451 			 * Euler's criterion.
1452 			 */
1453 			/* CTX->A has (a-2)/4 and CTX->B has b^-1 */
1454 			ec_mulm(w, ctx->a, mpi_const(MPI_C_FOUR), ctx);
1455 			ec_addm(w, w, mpi_const(MPI_C_TWO), ctx);
1456 			ec_mulm(w, w, x, ctx);
1457 			ec_pow2(xx, x, ctx);
1458 			ec_addm(w, w, xx, ctx);
1459 			ec_addm(w, w, mpi_const(MPI_C_ONE), ctx);
1460 			ec_mulm(w, w, x, ctx);
1461 			ec_mulm(w, w, ctx->b, ctx);
1462 #undef xx
1463 			/* Compute Euler's criterion: w^(p-1)/2 */
1464 #define p_minus1 y
1465 			ec_subm(p_minus1, ctx->p, mpi_const(MPI_C_ONE), ctx);
1466 			mpi_rshift(p_minus1, p_minus1, 1);
1467 			ec_powm(w, w, p_minus1, ctx);
1468 
1469 			res = !mpi_cmp_ui(w, 1);
1470 #undef p_minus1
1471 		}
1472 		break;
1473 
1474 	case MPI_EC_EDWARDS:
1475 		{
1476 			if (mpi_ec_get_affine(x, y, point, ctx))
1477 				goto leave;
1478 
1479 			mpi_resize(w, ctx->p->nlimbs);
1480 			w->nlimbs = ctx->p->nlimbs;
1481 
1482 			/* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
1483 			ctx->pow2(x, x, ctx);
1484 			ctx->pow2(y, y, ctx);
1485 			if (ctx->dialect == ECC_DIALECT_ED25519)
1486 				ctx->subm(w, ctx->p, x, ctx);
1487 			else
1488 				ctx->mulm(w, ctx->a, x, ctx);
1489 			ctx->addm(w, w, y, ctx);
1490 			ctx->mulm(x, x, y, ctx);
1491 			ctx->mulm(x, x, ctx->b, ctx);
1492 			ctx->subm(w, w, x, ctx);
1493 			if (!mpi_cmp_ui(w, 1))
1494 				res = 1;
1495 		}
1496 		break;
1497 	}
1498 
1499 leave:
1500 	mpi_free(w);
1501 	mpi_free(x);
1502 	mpi_free(y);
1503 
1504 	return res;
1505 }
1506 EXPORT_SYMBOL_GPL(mpi_ec_curve_point);
1507