xref: /freebsd/crypto/libecc/src/curves/prj_pt.c (revision f0865ec9906d5a18fa2a3b61381f22ce16e606ad)
1 /*
2  *  Copyright (C) 2017 - This file is part of libecc project
3  *
4  *  Authors:
5  *      Ryad BENADJILA <ryadbenadjila@gmail.com>
6  *      Arnaud EBALARD <arnaud.ebalard@ssi.gouv.fr>
7  *      Jean-Pierre FLORI <jean-pierre.flori@ssi.gouv.fr>
8  *
9  *  Contributors:
10  *      Nicolas VIVET <nicolas.vivet@ssi.gouv.fr>
11  *      Karim KHALFALLAH <karim.khalfallah@ssi.gouv.fr>
12  *
13  *  This software is licensed under a dual BSD and GPL v2 license.
14  *  See LICENSE file at the root folder of the project.
15  */
16 #include <libecc/curves/ec_shortw.h>
17 #include <libecc/curves/prj_pt.h>
18 #include <libecc/nn/nn_logical.h>
19 #include <libecc/nn/nn_add.h>
20 #include <libecc/nn/nn_rand.h>
21 #include <libecc/fp/fp_add.h>
22 #include <libecc/fp/fp_mul.h>
23 #include <libecc/fp/fp_montgomery.h>
24 #include <libecc/fp/fp_rand.h>
25 
26 #define PRJ_PT_MAGIC ((word_t)(0xe1cd70babb1d5afeULL))
27 
28 /*
29  * Check given projective point has been correctly initialized (using
30  * prj_pt_init()). Returns 0 on success, -1 on error.
31  */
prj_pt_check_initialized(prj_pt_src_t in)32 int prj_pt_check_initialized(prj_pt_src_t in)
33 {
34 	int ret;
35 
36 	MUST_HAVE(((in != NULL) && (in->magic == PRJ_PT_MAGIC)), ret, err);
37 	ret = ec_shortw_crv_check_initialized(in->crv);
38 
39 err:
40 	return ret;
41 }
42 
43 /*
44  * Initialize the projective point structure on given curve as the point at
45  * infinity. The function returns 0 on success, -1 on error.
46  */
prj_pt_init(prj_pt_t in,ec_shortw_crv_src_t curve)47 int prj_pt_init(prj_pt_t in, ec_shortw_crv_src_t curve)
48 {
49 	int ret;
50 
51 	ret = ec_shortw_crv_check_initialized(curve); EG(ret, err);
52 
53 	MUST_HAVE((in != NULL), ret, err);
54 
55 	ret = fp_init(&(in->X), curve->a.ctx); EG(ret, err);
56 	ret = fp_init(&(in->Y), curve->a.ctx); EG(ret, err);
57 	ret = fp_init(&(in->Z), curve->a.ctx); EG(ret, err);
58 	in->crv = curve;
59 	in->magic = PRJ_PT_MAGIC;
60 
61 err:
62 	return ret;
63 }
64 
65 /*
66  * Initialize the projective point structure on given curve using given
67  * coordinates. The function returns 0 on success, -1 on error.
68  */
prj_pt_init_from_coords(prj_pt_t in,ec_shortw_crv_src_t curve,fp_src_t xcoord,fp_src_t ycoord,fp_src_t zcoord)69 int prj_pt_init_from_coords(prj_pt_t in,
70 			     ec_shortw_crv_src_t curve,
71 			     fp_src_t xcoord, fp_src_t ycoord, fp_src_t zcoord)
72 {
73 	int ret;
74 
75 	ret = prj_pt_init(in, curve); EG(ret, err);
76 	ret = fp_copy(&(in->X), xcoord); EG(ret, err);
77 	ret = fp_copy(&(in->Y), ycoord); EG(ret, err);
78 	ret = fp_copy(&(in->Z), zcoord);
79 
80 err:
81 	return ret;
82 }
83 
84 /*
85  * Uninit given projective point structure. The function returns 0 on success,
86  * -1 on error. This is an error if passed point has not already been
87  * initialized first.
88  */
prj_pt_uninit(prj_pt_t in)89 void prj_pt_uninit(prj_pt_t in)
90 {
91 	if((in != NULL) && (in->magic == PRJ_PT_MAGIC) && (in->crv != NULL)){
92 		in->crv = NULL;
93 		in->magic = WORD(0);
94 
95 		fp_uninit(&(in->X));
96 		fp_uninit(&(in->Y));
97 		fp_uninit(&(in->Z));
98 	}
99 
100 	return;
101 }
102 
103 /*
104  * Checks if projective point is the point at infinity (last coordinate is 0).
105  * In that case, 'iszero' out parameter is set to 1. It is set to 0 if the
106  * point is not the point at infinity. The function returns 0 on success, -1 on
107  * error. On error, 'iszero' is not meaningful.
108  */
prj_pt_iszero(prj_pt_src_t in,int * iszero)109 int prj_pt_iszero(prj_pt_src_t in, int *iszero)
110 {
111 	int ret;
112 
113 	ret = prj_pt_check_initialized(in); EG(ret, err);
114 	ret = fp_iszero(&(in->Z), iszero);
115 
116 err:
117 	return ret;
118 }
119 
120 /*
121  * Set given projective point 'out' to the point at infinity. The functions
122  * returns 0 on success, -1 on error.
123  */
prj_pt_zero(prj_pt_t out)124 int prj_pt_zero(prj_pt_t out)
125 {
126 	int ret;
127 
128 	ret = prj_pt_check_initialized(out); EG(ret, err);
129 
130 	ret = fp_zero(&(out->X)); EG(ret, err);
131 	ret = fp_one(&(out->Y)); EG(ret, err);
132 	ret = fp_zero(&(out->Z));
133 
134 err:
135 	return ret;
136 }
137 
138 /*
139  * Check if a projective point is indeed on its curve. The function sets
140  * 'on_curve' out parameter to 1 if the point is on the curve, 0 if not.
141  * The function returns 0 on success, -1 on error. 'on_curve' is not
142  * meaningful on error.
143  */
prj_pt_is_on_curve(prj_pt_src_t in,int * on_curve)144 int prj_pt_is_on_curve(prj_pt_src_t in,  int *on_curve)
145 {
146 	int ret, cmp;
147 
148 	/* In order to check that we are on the curve, we
149 	 * use the projective formula of the curve:
150 	 *
151 	 *   Y**2 * Z = X**3 + a * X * Z**2 + b * Z**3
152 	 *
153 	 */
154 	fp X, Y, Z;
155 	X.magic = Y.magic = Z.magic = WORD(0);
156 
157 	ret = prj_pt_check_initialized(in); EG(ret, err);
158 	ret = ec_shortw_crv_check_initialized(in->crv); EG(ret, err);
159 	MUST_HAVE((on_curve != NULL), ret, err);
160 
161 	ret = fp_init(&X, in->X.ctx); EG(ret, err);
162 	ret = fp_init(&Y, in->X.ctx); EG(ret, err);
163 	ret = fp_init(&Z, in->X.ctx); EG(ret, err);
164 
165 	/* Compute X**3 + a * X * Z**2 + b * Z**3 on one side */
166 	ret = fp_sqr(&X, &(in->X)); EG(ret, err);
167 	ret = fp_mul(&X, &X, &(in->X)); EG(ret, err);
168 	ret = fp_mul(&Z, &(in->X), &(in->crv->a)); EG(ret, err);
169 	ret = fp_mul(&Y, &(in->crv->b), &(in->Z)); EG(ret, err);
170 	ret = fp_add(&Z, &Z, &Y); EG(ret, err);
171 	ret = fp_mul(&Z, &Z, &(in->Z)); EG(ret, err);
172 	ret = fp_mul(&Z, &Z, &(in->Z)); EG(ret, err);
173 	ret = fp_add(&X, &X, &Z); EG(ret, err);
174 
175 	/* Compute Y**2 * Z on the other side */
176 	ret = fp_sqr(&Y, &(in->Y)); EG(ret, err);
177 	ret = fp_mul(&Y, &Y, &(in->Z)); EG(ret, err);
178 
179 	/* Compare the two values */
180 	ret = fp_cmp(&X, &Y, &cmp); EG(ret, err);
181 
182 	(*on_curve) = (!cmp);
183 
184 err:
185 	fp_uninit(&X);
186 	fp_uninit(&Y);
187 	fp_uninit(&Z);
188 
189 	return ret;
190 }
191 
192 /*
193  * The function copies 'in' projective point to 'out'. 'out' is initialized by
194  * the function. The function returns 0 on sucess, -1 on error.
195  */
prj_pt_copy(prj_pt_t out,prj_pt_src_t in)196 int prj_pt_copy(prj_pt_t out, prj_pt_src_t in)
197 {
198 	int ret;
199 
200 	ret = prj_pt_check_initialized(in); EG(ret, err);
201 
202 	ret = prj_pt_init(out, in->crv); EG(ret, err);
203 
204 	ret = fp_copy(&(out->X), &(in->X)); EG(ret, err);
205 	ret = fp_copy(&(out->Y), &(in->Y)); EG(ret, err);
206 	ret = fp_copy(&(out->Z), &(in->Z));
207 
208 err:
209 	return ret;
210 }
211 
212 /*
213  * Convert given projective point 'in' to affine representation in 'out'. 'out'
214  * is initialized by the function. The function returns 0 on success, -1 on
215  * error. Passing the point at infinty to the function is considered as an
216  * error.
217  */
prj_pt_to_aff(aff_pt_t out,prj_pt_src_t in)218 int prj_pt_to_aff(aff_pt_t out, prj_pt_src_t in)
219 {
220 	int ret, iszero;
221 
222 	ret = prj_pt_check_initialized(in); EG(ret, err);
223 
224 	ret = prj_pt_iszero(in, &iszero); EG(ret, err);
225 	MUST_HAVE((!iszero), ret, err);
226 
227 	ret = aff_pt_init(out, in->crv); EG(ret, err);
228 
229 	ret = fp_inv(&(out->x), &(in->Z)); EG(ret, err);
230 	ret = fp_mul(&(out->y), &(in->Y), &(out->x)); EG(ret, err);
231 	ret = fp_mul(&(out->x), &(in->X), &(out->x));
232 
233 err:
234 	return ret;
235 }
236 
237 /*
238  * Get the unique Z = 1 projective point representation ("equivalent" to affine
239  * point). The function returns 0 on success, -1 on error.
240  */
prj_pt_unique(prj_pt_t out,prj_pt_src_t in)241 int prj_pt_unique(prj_pt_t out, prj_pt_src_t in)
242 {
243 	int ret, iszero;
244 
245 	ret = prj_pt_check_initialized(in); EG(ret, err);
246 	ret = prj_pt_iszero(in, &iszero); EG(ret, err);
247 	MUST_HAVE((!iszero), ret, err);
248 
249 	if(out == in){
250 		/* Aliasing case */
251 		fp tmp;
252 		tmp.magic = WORD(0);
253 
254 		ret = fp_init(&tmp, (in->Z).ctx); EG(ret, err);
255 		ret = fp_inv(&tmp, &(in->Z)); EG(ret, err1);
256 		ret = fp_mul(&(out->Y), &(in->Y), &tmp); EG(ret, err1);
257 		ret = fp_mul(&(out->X), &(in->X), &tmp); EG(ret, err1);
258 		ret = fp_one(&(out->Z)); EG(ret, err1);
259 err1:
260 		fp_uninit(&tmp); EG(ret, err);
261 	}
262 	else{
263 	        ret = prj_pt_init(out, in->crv); EG(ret, err);
264 		ret = fp_inv(&(out->X), &(in->Z)); EG(ret, err);
265 		ret = fp_mul(&(out->Y), &(in->Y), &(out->X)); EG(ret, err);
266 		ret = fp_mul(&(out->X), &(in->X), &(out->X)); EG(ret, err);
267 		ret = fp_one(&(out->Z)); EG(ret, err);
268 	}
269 
270 
271 err:
272 	return ret;
273 }
274 
275 /*
276  * Converts affine point 'in' to projective representation in 'out'. 'out' is
277  * initialized by the function. The function returns 0 on success, -1 on error.
278  */
ec_shortw_aff_to_prj(prj_pt_t out,aff_pt_src_t in)279 int ec_shortw_aff_to_prj(prj_pt_t out, aff_pt_src_t in)
280 {
281 	int ret, on_curve;
282 
283 	ret = aff_pt_check_initialized(in); EG(ret, err);
284 
285 	/* The input affine point must be on the curve */
286 	ret = aff_pt_is_on_curve(in, &on_curve); EG(ret, err);
287 	MUST_HAVE(on_curve, ret, err);
288 
289 	ret = prj_pt_init(out, in->crv); EG(ret, err);
290 	ret = fp_copy(&(out->X), &(in->x)); EG(ret, err);
291 	ret = fp_copy(&(out->Y), &(in->y)); EG(ret, err);
292 	ret = nn_one(&(out->Z).fp_val); /* Z = 1 */
293 
294 err:
295 	return ret;
296 }
297 
298 /*
299  * Compare projective points 'in1' and 'in2'. On success, 'cmp' is set to
300  * the result of the comparison (0 if in1 == in2, !0 if in1 != in2). The
301  * function returns 0 on success, -1 on error.
302  */
prj_pt_cmp(prj_pt_src_t in1,prj_pt_src_t in2,int * cmp)303 int prj_pt_cmp(prj_pt_src_t in1, prj_pt_src_t in2, int *cmp)
304 {
305 	fp X1, X2, Y1, Y2;
306 	int ret, x_cmp, y_cmp;
307 	X1.magic = X2.magic = Y1.magic = Y2.magic = WORD(0);
308 
309 	MUST_HAVE((cmp != NULL), ret, err);
310 	ret = prj_pt_check_initialized(in1); EG(ret, err);
311 	ret = prj_pt_check_initialized(in2); EG(ret, err);
312 
313 	MUST_HAVE((in1->crv == in2->crv), ret, err);
314 
315 	ret = fp_init(&X1, (in1->X).ctx); EG(ret, err);
316 	ret = fp_init(&X2, (in2->X).ctx); EG(ret, err);
317 	ret = fp_init(&Y1, (in1->Y).ctx); EG(ret, err);
318 	ret = fp_init(&Y2, (in2->Y).ctx); EG(ret, err);
319 
320 	/*
321 	 * Montgomery multiplication is used as it is faster than
322 	 * usual multiplication and the spurious multiplicative
323 	 * factor does not matter.
324 	 */
325 	ret = fp_mul_monty(&X1, &(in1->X), &(in2->Z)); EG(ret, err);
326 	ret = fp_mul_monty(&X2, &(in2->X), &(in1->Z)); EG(ret, err);
327 	ret = fp_mul_monty(&Y1, &(in1->Y), &(in2->Z)); EG(ret, err);
328 	ret = fp_mul_monty(&Y2, &(in2->Y), &(in1->Z)); EG(ret, err);
329 
330 	ret = fp_mul_monty(&X1, &(in1->X), &(in2->Z)); EG(ret, err);
331 	ret = fp_mul_monty(&X2, &(in2->X), &(in1->Z)); EG(ret, err);
332 	ret = fp_mul_monty(&Y1, &(in1->Y), &(in2->Z)); EG(ret, err);
333 	ret = fp_mul_monty(&Y2, &(in2->Y), &(in1->Z)); EG(ret, err);
334 	ret = fp_cmp(&X1, &X2, &x_cmp); EG(ret, err);
335 	ret = fp_cmp(&Y1, &Y2, &y_cmp);
336 
337 	if (!ret) {
338 		(*cmp) = (x_cmp | y_cmp);
339 	}
340 
341 err:
342 	fp_uninit(&Y2);
343 	fp_uninit(&Y1);
344 	fp_uninit(&X2);
345 	fp_uninit(&X1);
346 
347 	return ret;
348 }
349 
350 /*
351  * NOTE: this internal functions assumes that upper layer have checked that in1 and in2
352  * are initialized, and that cmp is not NULL.
353  */
_prj_pt_eq_or_opp_X(prj_pt_src_t in1,prj_pt_src_t in2,int * cmp)354 ATTRIBUTE_WARN_UNUSED_RET static inline int _prj_pt_eq_or_opp_X(prj_pt_src_t in1, prj_pt_src_t in2, int *cmp)
355 {
356 	int ret;
357 	fp X1, X2;
358 	X1.magic = X2.magic = WORD(0);
359 
360 	/*
361 	 * Montgomery multiplication is used as it is faster than
362 	 * usual multiplication and the spurious multiplicative
363 	 * factor does not matter.
364 	 */
365 	ret = fp_init(&X1, (in1->X).ctx); EG(ret, err);
366 	ret = fp_init(&X2, (in2->X).ctx); EG(ret, err);
367 	ret = fp_mul_monty(&X1, &(in1->X), &(in2->Z)); EG(ret, err);
368 	ret = fp_mul_monty(&X2, &(in2->X), &(in1->Z)); EG(ret, err);
369 	ret = fp_cmp(&X1, &X2, cmp);
370 
371 err:
372 	fp_uninit(&X1);
373 	fp_uninit(&X2);
374 
375 	return ret;
376 }
377 
378 /*
379  * NOTE: this internal functions assumes that upper layer have checked that in1 and in2
380  * are initialized, and that eq_or_opp is not NULL.
381  */
_prj_pt_eq_or_opp_Y(prj_pt_src_t in1,prj_pt_src_t in2,int * eq_or_opp)382 ATTRIBUTE_WARN_UNUSED_RET static inline int _prj_pt_eq_or_opp_Y(prj_pt_src_t in1, prj_pt_src_t in2, int *eq_or_opp)
383 {
384 	int ret;
385 	fp Y1, Y2;
386 	Y1.magic = Y2.magic = WORD(0);
387 
388 	/*
389 	 * Montgomery multiplication is used as it is faster than
390 	 * usual multiplication and the spurious multiplicative
391 	 * factor does not matter.
392 	 */
393 	ret = fp_init(&Y1, (in1->Y).ctx); EG(ret, err);
394 	ret = fp_init(&Y2, (in2->Y).ctx); EG(ret, err);
395 	ret = fp_mul_monty(&Y1, &(in1->Y), &(in2->Z)); EG(ret, err);
396 	ret = fp_mul_monty(&Y2, &(in2->Y), &(in1->Z)); EG(ret, err);
397 	ret = fp_eq_or_opp(&Y1, &Y2, eq_or_opp);
398 
399 err:
400 	fp_uninit(&Y1);
401 	fp_uninit(&Y2);
402 
403 	return ret;
404 }
405 
406  /*
407  * The functions tests if given projective points 'in1' and 'in2' are equal or
408  * opposite. On success, the result of the comparison is given via 'eq_or_opp'
409  * out parameter (1 if equal or opposite, 0 otherwise). The function returns
410  * 0 on succes, -1 on error.
411  */
prj_pt_eq_or_opp(prj_pt_src_t in1,prj_pt_src_t in2,int * eq_or_opp)412 int prj_pt_eq_or_opp(prj_pt_src_t in1, prj_pt_src_t in2, int *eq_or_opp)
413 {
414 	int ret, cmp, _eq_or_opp;
415 
416 	ret = prj_pt_check_initialized(in1); EG(ret, err);
417 	ret = prj_pt_check_initialized(in2); EG(ret, err);
418 	MUST_HAVE((in1->crv == in2->crv), ret, err);
419 	MUST_HAVE((eq_or_opp != NULL), ret, err);
420 
421 	ret = _prj_pt_eq_or_opp_X(in1, in2, &cmp); EG(ret, err);
422 	ret = _prj_pt_eq_or_opp_Y(in1, in2, &_eq_or_opp);
423 
424 	if (!ret) {
425 		(*eq_or_opp) = ((cmp == 0) & _eq_or_opp);
426 	}
427 
428 err:
429 	return ret;
430 }
431 
432 /* Compute the opposite of a projective point. Supports aliasing.
433  * Returns 0 on success, -1 on failure.
434  */
prj_pt_neg(prj_pt_t out,prj_pt_src_t in)435 int prj_pt_neg(prj_pt_t out, prj_pt_src_t in)
436 {
437 	int ret;
438 
439 	ret = prj_pt_check_initialized(in); EG(ret, err);
440 
441 	if (out != in) { /* Copy point if not aliased */
442 		ret = prj_pt_init(out, in->crv); EG(ret, err);
443 		ret = prj_pt_copy(out, in); EG(ret, err);
444 	}
445 
446 	/* Then, negate Y */
447 	ret = fp_neg(&(out->Y), &(out->Y));
448 
449 err:
450 	return ret;
451 }
452 
453 /*
454  * Import a projective point from a buffer with the following layout; the 3
455  * coordinates (elements of Fp) are each encoded on p_len bytes, where p_len
456  * is the size of p in bytes (e.g. 66 for a prime p of 521 bits). Each
457  * coordinate is encoded in big endian. Size of buffer must exactly match
458  * 3 * p_len. The projective point is initialized by the function.
459  *
460  * The function returns 0 on success, -1 on error.
461  */
prj_pt_import_from_buf(prj_pt_t pt,const u8 * pt_buf,u16 pt_buf_len,ec_shortw_crv_src_t crv)462 int prj_pt_import_from_buf(prj_pt_t pt,
463 			   const u8 *pt_buf,
464 			   u16 pt_buf_len, ec_shortw_crv_src_t crv)
465 {
466 	int on_curve, ret;
467 	fp_ctx_src_t ctx;
468 	u16 coord_len;
469 
470 	ret = ec_shortw_crv_check_initialized(crv); EG(ret, err);
471 	MUST_HAVE((pt_buf != NULL) && (pt != NULL), ret, err);
472 
473 	ctx = crv->a.ctx;
474 	coord_len = (u16)BYTECEIL(ctx->p_bitlen);
475 	MUST_HAVE((pt_buf_len == (3 * coord_len)), ret, err);
476 
477 	ret = fp_init_from_buf(&(pt->X), ctx, pt_buf, coord_len); EG(ret, err);
478 	ret = fp_init_from_buf(&(pt->Y), ctx, pt_buf + coord_len, coord_len); EG(ret, err);
479 	ret = fp_init_from_buf(&(pt->Z), ctx, pt_buf + (2 * coord_len), coord_len); EG(ret, err);
480 
481 	/* Set the curve */
482 	pt->crv = crv;
483 
484 	/* Mark the point as initialized */
485 	pt->magic = PRJ_PT_MAGIC;
486 
487 	/* Check that the point is indeed on the provided curve, uninitialize it
488 	 * if this is not the case.
489 	 */
490 	ret = prj_pt_is_on_curve(pt, &on_curve); EG(ret, err);
491 	if (!on_curve){
492 		prj_pt_uninit(pt);
493 		ret = -1;
494 	}
495 
496 err:
497 	PTR_NULLIFY(ctx);
498 
499 	return ret;
500 }
501 
502 /*
503  * Import a projective point from an affine point buffer with the following layout; the 2
504  * coordinates (elements of Fp) are each encoded on p_len bytes, where p_len
505  * is the size of p in bytes (e.g. 66 for a prime p of 521 bits). Each
506  * coordinate is encoded in big endian. Size of buffer must exactly match
507  * 2 * p_len. The projective point is initialized by the function.
508  *
509  * The function returns 0 on success, -1 on error.
510  */
prj_pt_import_from_aff_buf(prj_pt_t pt,const u8 * pt_buf,u16 pt_buf_len,ec_shortw_crv_src_t crv)511 int prj_pt_import_from_aff_buf(prj_pt_t pt,
512 			   const u8 *pt_buf,
513 			   u16 pt_buf_len, ec_shortw_crv_src_t crv)
514 {
515 	int ret, on_curve;
516 	fp_ctx_src_t ctx;
517 	u16 coord_len;
518 
519 	ret = ec_shortw_crv_check_initialized(crv); EG(ret, err);
520 	MUST_HAVE((pt_buf != NULL) && (pt != NULL), ret, err);
521 
522 	ctx = crv->a.ctx;
523 	coord_len = (u16)BYTECEIL(ctx->p_bitlen);
524 	MUST_HAVE((pt_buf_len == (2 * coord_len)), ret, err);
525 
526 	ret = fp_init_from_buf(&(pt->X), ctx, pt_buf, coord_len); EG(ret, err);
527 	ret = fp_init_from_buf(&(pt->Y), ctx, pt_buf + coord_len, coord_len); EG(ret, err);
528 	/* Z coordinate is set to 1 */
529 	ret = fp_init(&(pt->Z), ctx); EG(ret, err);
530 	ret = fp_one(&(pt->Z)); EG(ret, err);
531 
532 	/* Set the curve */
533 	pt->crv = crv;
534 
535 	/* Mark the point as initialized */
536 	pt->magic = PRJ_PT_MAGIC;
537 
538 	/* Check that the point is indeed on the provided curve, uninitialize it
539 	 * if this is not the case.
540 	 */
541 	ret = prj_pt_is_on_curve(pt, &on_curve); EG(ret, err);
542 	if (!on_curve){
543 		prj_pt_uninit(pt);
544 		ret = -1;
545 	}
546 
547 err:
548 	PTR_NULLIFY(ctx);
549 
550 	return ret;
551 }
552 
553 
554 /* Export a projective point to a buffer with the following layout; the 3
555  * coordinates (elements of Fp) are each encoded on p_len bytes, where p_len
556  * is the size of p in bytes (e.g. 66 for a prime p of 521 bits). Each
557  * coordinate is encoded in big endian. Size of buffer must exactly match
558  * 3 * p_len.
559  *
560  * The function returns 0 on success, -1 on error.
561  */
prj_pt_export_to_buf(prj_pt_src_t pt,u8 * pt_buf,u32 pt_buf_len)562 int prj_pt_export_to_buf(prj_pt_src_t pt, u8 *pt_buf, u32 pt_buf_len)
563 {
564 	fp_ctx_src_t ctx;
565 	u16 coord_len;
566 	int ret, on_curve;
567 
568 	ret = prj_pt_check_initialized(pt); EG(ret, err);
569 
570 	MUST_HAVE((pt_buf != NULL), ret, err);
571 
572 	/* The point to be exported must be on the curve */
573 	ret = prj_pt_is_on_curve(pt, &on_curve); EG(ret, err);
574 	MUST_HAVE((on_curve), ret, err);
575 
576 	ctx = pt->crv->a.ctx;
577 	coord_len = (u16)BYTECEIL(ctx->p_bitlen);
578 	MUST_HAVE((pt_buf_len == (3 * coord_len)), ret, err);
579 
580 	/* Export the three coordinates */
581 	ret = fp_export_to_buf(pt_buf, coord_len, &(pt->X)); EG(ret, err);
582 	ret = fp_export_to_buf(pt_buf + coord_len, coord_len, &(pt->Y)); EG(ret, err);
583 	ret = fp_export_to_buf(pt_buf + (2 * coord_len), coord_len, &(pt->Z));
584 
585 err:
586 	PTR_NULLIFY(ctx);
587 
588 	return ret;
589 }
590 
591 /*
592  * Export a projective point to an affine point buffer with the following
593  * layout; the 2 coordinates (elements of Fp) are each encoded on p_len bytes,
594  * where p_len is the size of p in bytes (e.g. 66 for a prime p of 521 bits).
595  * Each coordinate is encoded in big endian. Size of buffer must exactly match
596  * 2 * p_len.
597  *
598  * The function returns 0 on success, -1 on error.
599  */
prj_pt_export_to_aff_buf(prj_pt_src_t pt,u8 * pt_buf,u32 pt_buf_len)600 int prj_pt_export_to_aff_buf(prj_pt_src_t pt, u8 *pt_buf, u32 pt_buf_len)
601 {
602 	int ret, on_curve;
603 	aff_pt tmp_aff;
604 	tmp_aff.magic = WORD(0);
605 
606 	ret = prj_pt_check_initialized(pt); EG(ret, err);
607 
608 	MUST_HAVE((pt_buf != NULL), ret, err);
609 
610 	/* The point to be exported must be on the curve */
611 	ret = prj_pt_is_on_curve(pt, &on_curve); EG(ret, err);
612 	MUST_HAVE((on_curve), ret, err);
613 
614 	/* Move to the affine unique representation */
615 	ret = prj_pt_to_aff(&tmp_aff, pt); EG(ret, err);
616 
617 	/* Export the affine point to the buffer */
618 	ret = aff_pt_export_to_buf(&tmp_aff, pt_buf, pt_buf_len);
619 
620 err:
621 	aff_pt_uninit(&tmp_aff);
622 
623 	return ret;
624 }
625 
626 
627 #ifdef NO_USE_COMPLETE_FORMULAS
628 
629 /*
630  * The function is an internal one: no check is performed on parameters,
631  * this MUST be done by the caller:
632  *
633  *  - in is initialized
634  *  - in and out must not be aliased
635  *
636  * The function will initialize 'out'. The function returns 0 on success, -1
637  * on error.
638  */
__prj_pt_dbl_monty_no_cf(prj_pt_t out,prj_pt_src_t in)639 ATTRIBUTE_WARN_UNUSED_RET static int __prj_pt_dbl_monty_no_cf(prj_pt_t out, prj_pt_src_t in)
640 {
641 	fp XX, ZZ, w, s, ss, sss, R, RR, B, h;
642 	int ret;
643 	XX.magic = ZZ.magic = w.magic = s.magic = WORD(0);
644 	ss.magic = sss.magic = R.magic = WORD(0);
645 	RR.magic = B.magic = h.magic = WORD(0);
646 
647 	ret = prj_pt_init(out, in->crv); EG(ret, err);
648 
649 	ret = fp_init(&XX, out->crv->a.ctx); EG(ret, err);
650 	ret = fp_init(&ZZ, out->crv->a.ctx); EG(ret, err);
651 	ret = fp_init(&w, out->crv->a.ctx); EG(ret, err);
652 	ret = fp_init(&s, out->crv->a.ctx); EG(ret, err);
653 	ret = fp_init(&ss, out->crv->a.ctx); EG(ret, err);
654 	ret = fp_init(&sss, out->crv->a.ctx); EG(ret, err);
655 	ret = fp_init(&R, out->crv->a.ctx); EG(ret, err);
656 	ret = fp_init(&RR, out->crv->a.ctx); EG(ret, err);
657 	ret = fp_init(&B, out->crv->a.ctx); EG(ret, err);
658 	ret = fp_init(&h, out->crv->a.ctx); EG(ret, err);
659 
660 	/* XX = X1² */
661 	ret = fp_sqr_monty(&XX, &(in->X)); EG(ret, err);
662 
663 	/* ZZ = Z1² */
664 	ret = fp_sqr_monty(&ZZ, &(in->Z)); EG(ret, err);
665 
666 	/* w = a*ZZ+3*XX */
667 	ret = fp_mul_monty(&w, &(in->crv->a_monty), &ZZ); EG(ret, err);
668 	ret = fp_add_monty(&w, &w, &XX); EG(ret, err);
669 	ret = fp_add_monty(&w, &w, &XX); EG(ret, err);
670 	ret = fp_add_monty(&w, &w, &XX); EG(ret, err);
671 
672 	/* s = 2*Y1*Z1 */
673 	ret = fp_mul_monty(&s, &(in->Y), &(in->Z)); EG(ret, err);
674 	ret = fp_add_monty(&s, &s, &s); EG(ret, err);
675 
676 	/* ss = s² */
677 	ret = fp_sqr_monty(&ss, &s); EG(ret, err);
678 
679 	/* sss = s*ss */
680 	ret = fp_mul_monty(&sss, &s, &ss); EG(ret, err);
681 
682 	/* R = Y1*s */
683 	ret = fp_mul_monty(&R, &(in->Y), &s); EG(ret, err);
684 
685 	/* RR = R² */
686 	ret = fp_sqr_monty(&RR, &R); EG(ret, err);
687 
688 	/* B = (X1+R)²-XX-RR */
689 	ret = fp_add_monty(&R, &R, &(in->X)); EG(ret, err);
690 	ret = fp_sqr_monty(&B, &R); EG(ret, err);
691 	ret = fp_sub_monty(&B, &B, &XX); EG(ret, err);
692 	ret = fp_sub_monty(&B, &B, &RR); EG(ret, err);
693 
694 	/* h = w²-2*B */
695 	ret = fp_sqr_monty(&h, &w); EG(ret, err);
696 	ret = fp_sub_monty(&h, &h, &B); EG(ret, err);
697 	ret = fp_sub_monty(&h, &h, &B); EG(ret, err);
698 
699 	/* X3 = h*s */
700 	ret = fp_mul_monty(&(out->X), &h, &s); EG(ret, err);
701 
702 	/* Y3 = w*(B-h)-2*RR */
703 	ret = fp_sub_monty(&B, &B, &h); EG(ret, err);
704 	ret = fp_mul_monty(&(out->Y), &w, &B); EG(ret, err);
705 	ret = fp_sub_monty(&(out->Y), &(out->Y), &RR); EG(ret, err);
706 	ret = fp_sub_monty(&(out->Y), &(out->Y), &RR); EG(ret, err);
707 
708 	/* Z3 = sss */
709 	ret = fp_copy(&(out->Z), &sss);
710 
711 err:
712 	fp_uninit(&XX);
713 	fp_uninit(&ZZ);
714 	fp_uninit(&w);
715 	fp_uninit(&s);
716 	fp_uninit(&ss);
717 	fp_uninit(&sss);
718 	fp_uninit(&R);
719 	fp_uninit(&RR);
720 	fp_uninit(&B);
721 	fp_uninit(&h);
722 
723 	return ret;
724 }
725 
726 /*
727  * The function is an internal one: no check is performed on parameters,
728  * this MUST be done by the caller:
729  *
730  *  - in1 and in2 are initialized
731  *  - in1 and in2 are on the same curve
732  *  - in1/in2 and out must not be aliased
733  *  - in1 and in2 must not be equal, opposite or have identical value
734  *
735  * The function will initialize 'out'. The function returns 0 on success, -1
736  * on error.
737  */
___prj_pt_add_monty_no_cf(prj_pt_t out,prj_pt_src_t in1,prj_pt_src_t in2)738 ATTRIBUTE_WARN_UNUSED_RET static int ___prj_pt_add_monty_no_cf(prj_pt_t out,
739 							       prj_pt_src_t in1,
740 							       prj_pt_src_t in2)
741 {
742 	fp Y1Z2, X1Z2, Z1Z2, u, uu, v, vv, vvv, R, A;
743 	int ret;
744 	Y1Z2.magic = X1Z2.magic = Z1Z2.magic = u.magic = uu.magic = v.magic = WORD(0);
745 	vv.magic = vvv.magic = R.magic = A.magic = WORD(0);
746 
747 	ret = prj_pt_init(out, in1->crv); EG(ret, err);
748 
749 	ret = fp_init(&Y1Z2, out->crv->a.ctx); EG(ret, err);
750 	ret = fp_init(&X1Z2, out->crv->a.ctx); EG(ret, err);
751 	ret = fp_init(&Z1Z2, out->crv->a.ctx); EG(ret, err);
752 	ret = fp_init(&u, out->crv->a.ctx); EG(ret, err);
753 	ret = fp_init(&uu, out->crv->a.ctx); EG(ret, err);
754 	ret = fp_init(&v, out->crv->a.ctx); EG(ret, err);
755 	ret = fp_init(&vv, out->crv->a.ctx); EG(ret, err);
756 	ret = fp_init(&vvv, out->crv->a.ctx); EG(ret, err);
757 	ret = fp_init(&R, out->crv->a.ctx); EG(ret, err);
758 	ret = fp_init(&A, out->crv->a.ctx); EG(ret, err);
759 
760 	/* Y1Z2 = Y1*Z2 */
761 	ret = fp_mul_monty(&Y1Z2, &(in1->Y), &(in2->Z)); EG(ret, err);
762 
763 	/* X1Z2 = X1*Z2 */
764 	ret = fp_mul_monty(&X1Z2, &(in1->X), &(in2->Z)); EG(ret, err);
765 
766 	/* Z1Z2 = Z1*Z2 */
767 	ret = fp_mul_monty(&Z1Z2, &(in1->Z), &(in2->Z)); EG(ret, err);
768 
769 	/* u = Y2*Z1-Y1Z2 */
770 	ret = fp_mul_monty(&u, &(in2->Y), &(in1->Z)); EG(ret, err);
771 	ret = fp_sub_monty(&u, &u, &Y1Z2); EG(ret, err);
772 
773 	/* uu = u² */
774 	ret = fp_sqr_monty(&uu, &u); EG(ret, err);
775 
776 	/* v = X2*Z1-X1Z2 */
777 	ret = fp_mul_monty(&v, &(in2->X), &(in1->Z)); EG(ret, err);
778 	ret = fp_sub_monty(&v, &v, &X1Z2); EG(ret, err);
779 
780 	/* vv = v² */
781 	ret = fp_sqr_monty(&vv, &v); EG(ret, err);
782 
783 	/* vvv = v*vv */
784 	ret = fp_mul_monty(&vvv, &v, &vv); EG(ret, err);
785 
786 	/* R = vv*X1Z2 */
787 	ret = fp_mul_monty(&R, &vv, &X1Z2); EG(ret, err);
788 
789 	/* A = uu*Z1Z2-vvv-2*R */
790 	ret = fp_mul_monty(&A, &uu, &Z1Z2); EG(ret, err);
791 	ret = fp_sub_monty(&A, &A, &vvv); EG(ret, err);
792 	ret = fp_sub_monty(&A, &A, &R); EG(ret, err);
793 	ret = fp_sub_monty(&A, &A, &R); EG(ret, err);
794 
795 	/* X3 = v*A */
796 	ret = fp_mul_monty(&(out->X), &v, &A); EG(ret, err);
797 
798 	/* Y3 = u*(R-A)-vvv*Y1Z2 */
799 	ret = fp_sub_monty(&R, &R, &A); EG(ret, err);
800 	ret = fp_mul_monty(&(out->Y), &u, &R); EG(ret, err);
801 	ret = fp_mul_monty(&R, &vvv, &Y1Z2); EG(ret, err);
802 	ret = fp_sub_monty(&(out->Y), &(out->Y), &R); EG(ret, err);
803 
804 	/* Z3 = vvv*Z1Z2 */
805 	ret = fp_mul_monty(&(out->Z), &vvv, &Z1Z2);
806 
807 err:
808 	fp_uninit(&Y1Z2);
809 	fp_uninit(&X1Z2);
810 	fp_uninit(&Z1Z2);
811 	fp_uninit(&u);
812 	fp_uninit(&uu);
813 	fp_uninit(&v);
814 	fp_uninit(&vv);
815 	fp_uninit(&vvv);
816 	fp_uninit(&R);
817 	fp_uninit(&A);
818 
819 	return ret;
820 }
821 
822 /*
823  * Public version of the addition w/o complete formulas to handle the case
824  * where the inputs are zero or opposite. Returns 0 on success, -1 on error.
825  */
__prj_pt_add_monty_no_cf(prj_pt_t out,prj_pt_src_t in1,prj_pt_src_t in2)826 ATTRIBUTE_WARN_UNUSED_RET static int __prj_pt_add_monty_no_cf(prj_pt_t out, prj_pt_src_t in1, prj_pt_src_t in2)
827 {
828 	int ret, iszero, eq_or_opp, cmp;
829 
830 	ret = prj_pt_check_initialized(in1); EG(ret, err);
831 	ret = prj_pt_check_initialized(in2); EG(ret, err);
832 	MUST_HAVE((in1->crv == in2->crv), ret, err);
833 
834 	ret = prj_pt_iszero(in1, &iszero); EG(ret, err);
835 	if (iszero) {
836 		/* in1 at infinity, output in2 in all cases */
837 		ret = prj_pt_init(out, in2->crv); EG(ret, err);
838 		ret = prj_pt_copy(out, in2); EG(ret, err);
839 	} else {
840 		/* in1 not at infinity, output in2 */
841 		ret = prj_pt_iszero(in2, &iszero); EG(ret, err);
842 		if (iszero) {
843 			/* in2 at infinity, output in1 */
844 			ret = prj_pt_init(out, in1->crv); EG(ret, err);
845 			ret = prj_pt_copy(out, in1); EG(ret, err);
846 		} else {
847 			/* enither in1, nor in2 at infinity */
848 
849 			/*
850 			 * The following test which guarantees in1 and in2 are not
851 			 * equal or opposite needs to be rewritten because it
852 			 * has a *HUGE* impact on perf (ec_self_tests run on
853 			 * all test vectors takes 24 times as long with this
854 			 * enabled). The same exists in non monty version.
855 			 */
856 			ret = prj_pt_eq_or_opp(in1, in2, &eq_or_opp); EG(ret, err);
857 			if (eq_or_opp) {
858 				/* in1 and in2 are either equal or opposite */
859 				ret = prj_pt_cmp(in1, in2, &cmp); EG(ret, err);
860 				if (cmp == 0) {
861 					/* in1 == in2 => doubling w/o cf */
862 					ret = __prj_pt_dbl_monty_no_cf(out, in1); EG(ret, err);
863 				} else {
864 					/* in1 == -in2 => output zero (point at infinity) */
865 					ret = prj_pt_init(out, in1->crv); EG(ret, err);
866 					ret = prj_pt_zero(out); EG(ret, err);
867 				}
868 			} else {
869 				/*
870 				 * in1 and in2 are neither 0, nor equal or
871 				 * opposite. Use the basic monty addition
872 				 * implementation w/o complete formulas.
873 				 */
874 				ret = ___prj_pt_add_monty_no_cf(out, in1, in2); EG(ret, err);
875 			}
876 		}
877 	}
878 
879 err:
880 	return ret;
881 }
882 
883 
884 #else /* NO_USE_COMPLETE_FORMULAS */
885 
886 
887 /*
888  * If NO_USE_COMPLETE_FORMULAS flag is not defined addition formulas from Algorithm 3
889  * of https://joostrenes.nl/publications/complete.pdf are used, otherwise
890  * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#doubling-dbl-2007-bl
891  */
__prj_pt_dbl_monty_cf(prj_pt_t out,prj_pt_src_t in)892 ATTRIBUTE_WARN_UNUSED_RET static int __prj_pt_dbl_monty_cf(prj_pt_t out, prj_pt_src_t in)
893 {
894 	fp t0, t1, t2, t3;
895 	int ret;
896 	t0.magic = t1.magic = t2.magic = t3.magic = WORD(0);
897 
898 	ret = prj_pt_init(out, in->crv); EG(ret, err);
899 
900 	ret = fp_init(&t0, out->crv->a.ctx); EG(ret, err);
901 	ret = fp_init(&t1, out->crv->a.ctx); EG(ret, err);
902 	ret = fp_init(&t2, out->crv->a.ctx); EG(ret, err);
903 	ret = fp_init(&t3, out->crv->a.ctx); EG(ret, err);
904 
905 	ret = fp_mul_monty(&t0, &in->X, &in->X); EG(ret, err);
906 	ret = fp_mul_monty(&t1, &in->Y, &in->Y); EG(ret, err);
907 	ret = fp_mul_monty(&t2, &in->Z, &in->Z); EG(ret, err);
908 	ret = fp_mul_monty(&t3, &in->X, &in->Y); EG(ret, err);
909 	ret = fp_add_monty(&t3, &t3, &t3); EG(ret, err);
910 
911 	ret = fp_mul_monty(&out->Z, &in->X, &in->Z); EG(ret, err);
912 	ret = fp_add_monty(&out->Z, &out->Z, &out->Z); EG(ret, err);
913 	ret = fp_mul_monty(&out->X, &in->crv->a_monty, &out->Z); EG(ret, err);
914 	ret = fp_mul_monty(&out->Y, &in->crv->b3_monty, &t2); EG(ret, err);
915 	ret = fp_add_monty(&out->Y, &out->X, &out->Y); EG(ret, err);
916 
917 	ret = fp_sub_monty(&out->X, &t1, &out->Y); EG(ret, err);
918 	ret = fp_add_monty(&out->Y, &t1, &out->Y); EG(ret, err);
919 	ret = fp_mul_monty(&out->Y, &out->X, &out->Y); EG(ret, err);
920 	ret = fp_mul_monty(&out->X, &t3, &out->X); EG(ret, err);
921 	ret = fp_mul_monty(&out->Z, &in->crv->b3_monty, &out->Z); EG(ret, err);
922 
923 	ret = fp_mul_monty(&t2, &in->crv->a_monty, &t2); EG(ret, err);
924 	ret = fp_sub_monty(&t3, &t0, &t2); EG(ret, err);
925 	ret = fp_mul_monty(&t3, &in->crv->a_monty, &t3); EG(ret, err);
926 	ret = fp_add_monty(&t3, &t3, &out->Z); EG(ret, err);
927 	ret = fp_add_monty(&out->Z, &t0, &t0); EG(ret, err);
928 
929 	ret = fp_add_monty(&t0, &out->Z, &t0); EG(ret, err);
930 	ret = fp_add_monty(&t0, &t0, &t2); EG(ret, err);
931 	ret = fp_mul_monty(&t0, &t0, &t3); EG(ret, err);
932 	ret = fp_add_monty(&out->Y, &out->Y, &t0); EG(ret, err);
933 	ret = fp_mul_monty(&t2, &in->Y, &in->Z); EG(ret, err);
934 
935 	ret = fp_add_monty(&t2, &t2, &t2); EG(ret, err);
936 	ret = fp_mul_monty(&t0, &t2, &t3); EG(ret, err);
937 	ret = fp_sub_monty(&out->X, &out->X, &t0); EG(ret, err);
938 	ret = fp_mul_monty(&out->Z, &t2, &t1); EG(ret, err);
939 	ret = fp_add_monty(&out->Z, &out->Z, &out->Z); EG(ret, err);
940 
941 	ret = fp_add_monty(&out->Z, &out->Z, &out->Z);
942 
943 err:
944 	fp_uninit(&t0);
945 	fp_uninit(&t1);
946 	fp_uninit(&t2);
947 	fp_uninit(&t3);
948 
949 	return ret;
950 }
951 
952 /*
953  * If NO_USE_COMPLETE_FORMULAS flag is not defined addition formulas from Algorithm 1
954  * of https://joostrenes.nl/publications/complete.pdf are used, otherwise
955  * http://www.hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-1998-cmo-2
956  */
957 
958 /*
959  * The function is an internal one: no check is performed on parameters,
960  * this MUST be done by the caller:
961  *
962  *  - in1 and in2 are initialized
963  *  - in1 and in2 are on the same curve
964  *  - in1/in2 and out must not be aliased
965  *  - in1 and in2 must not be an "exceptional" pair, i.e. (in1-in2) is not a point
966  *  of order exactly 2
967  *
968  * The function will initialize 'out'. The function returns 0 on success, -1
969  * on error.
970  */
__prj_pt_add_monty_cf(prj_pt_t out,prj_pt_src_t in1,prj_pt_src_t in2)971 ATTRIBUTE_WARN_UNUSED_RET static int __prj_pt_add_monty_cf(prj_pt_t out,
972 							   prj_pt_src_t in1,
973 							   prj_pt_src_t in2)
974 {
975 	int cmp1, cmp2;
976 	fp t0, t1, t2, t3, t4, t5;
977 	int ret;
978 	t0.magic = t1.magic = t2.magic = WORD(0);
979 	t3.magic = t4.magic = t5.magic = WORD(0);
980 
981 	ret = prj_pt_init(out, in1->crv); EG(ret, err);
982 
983 	ret = fp_init(&t0, out->crv->a.ctx); EG(ret, err);
984 	ret = fp_init(&t1, out->crv->a.ctx); EG(ret, err);
985 	ret = fp_init(&t2, out->crv->a.ctx); EG(ret, err);
986 	ret = fp_init(&t3, out->crv->a.ctx); EG(ret, err);
987 	ret = fp_init(&t4, out->crv->a.ctx); EG(ret, err);
988 	ret = fp_init(&t5, out->crv->a.ctx); EG(ret, err);
989 
990 	ret = fp_mul_monty(&t0, &in1->X, &in2->X); EG(ret, err);
991 	ret = fp_mul_monty(&t1, &in1->Y, &in2->Y); EG(ret, err);
992 	ret = fp_mul_monty(&t2, &in1->Z, &in2->Z); EG(ret, err);
993 	ret = fp_add_monty(&t3, &in1->X, &in1->Y); EG(ret, err);
994 	ret = fp_add_monty(&t4, &in2->X, &in2->Y); EG(ret, err);
995 
996 	ret = fp_mul_monty(&t3, &t3, &t4); EG(ret, err);
997 	ret = fp_add_monty(&t4, &t0, &t1); EG(ret, err);
998 	ret = fp_sub_monty(&t3, &t3, &t4); EG(ret, err);
999 	ret = fp_add_monty(&t4, &in1->X, &in1->Z); EG(ret, err);
1000 	ret = fp_add_monty(&t5, &in2->X, &in2->Z); EG(ret, err);
1001 
1002 	ret = fp_mul_monty(&t4, &t4, &t5); EG(ret, err);
1003 	ret = fp_add_monty(&t5, &t0, &t2); EG(ret, err);
1004 	ret = fp_sub_monty(&t4, &t4, &t5); EG(ret, err);
1005 	ret = fp_add_monty(&t5, &in1->Y, &in1->Z); EG(ret, err);
1006 	ret = fp_add_monty(&out->X, &in2->Y, &in2->Z); EG(ret, err);
1007 
1008 	ret = fp_mul_monty(&t5, &t5, &out->X); EG(ret, err);
1009 	ret = fp_add_monty(&out->X, &t1, &t2); EG(ret, err);
1010 	ret = fp_sub_monty(&t5, &t5, &out->X); EG(ret, err);
1011 	ret = fp_mul_monty(&out->Z, &in1->crv->a_monty, &t4); EG(ret, err);
1012 	ret = fp_mul_monty(&out->X, &in1->crv->b3_monty, &t2); EG(ret, err);
1013 
1014 	ret = fp_add_monty(&out->Z, &out->X, &out->Z); EG(ret, err);
1015 	ret = fp_sub_monty(&out->X, &t1, &out->Z); EG(ret, err);
1016 	ret = fp_add_monty(&out->Z, &t1, &out->Z); EG(ret, err);
1017 	ret = fp_mul_monty(&out->Y, &out->X, &out->Z); EG(ret, err);
1018 	ret = fp_add_monty(&t1, &t0, &t0); EG(ret, err);
1019 
1020 	ret = fp_add_monty(&t1, &t1, &t0); EG(ret, err);
1021 	ret = fp_mul_monty(&t2, &in1->crv->a_monty, &t2); EG(ret, err);
1022 	ret = fp_mul_monty(&t4, &in1->crv->b3_monty, &t4); EG(ret, err);
1023 	ret = fp_add_monty(&t1, &t1, &t2); EG(ret, err);
1024 	ret = fp_sub_monty(&t2, &t0, &t2); EG(ret, err);
1025 
1026 	ret = fp_mul_monty(&t2, &in1->crv->a_monty, &t2); EG(ret, err);
1027 	ret = fp_add_monty(&t4, &t4, &t2); EG(ret, err);
1028 	ret = fp_mul_monty(&t0, &t1, &t4); EG(ret, err);
1029 	ret = fp_add_monty(&out->Y, &out->Y, &t0); EG(ret, err);
1030 	ret = fp_mul_monty(&t0, &t5, &t4); EG(ret, err);
1031 
1032 	ret = fp_mul_monty(&out->X, &t3, &out->X); EG(ret, err);
1033 	ret = fp_sub_monty(&out->X, &out->X, &t0); EG(ret, err);
1034 	ret = fp_mul_monty(&t0, &t3, &t1); EG(ret, err);
1035 	ret = fp_mul_monty(&out->Z, &t5, &out->Z); EG(ret, err);
1036 	ret = fp_add_monty(&out->Z, &out->Z, &t0);
1037 
1038 	/* Check for "exceptional" pairs of input points with
1039 	 * checking if Y = Z = 0 as output (see the Bosma-Lenstra
1040 	 * article "Complete Systems of Two Addition Laws for
1041 	 * Elliptic Curves"). This should only happen on composite
1042 	 * order curves (i.e. not on prime order curves).
1043 	 *
1044 	 * In this case, we raise an error as the result is
1045 	 * not sound. This should not happen in our nominal
1046 	 * cases with regular signature and protocols, and if
1047 	 * it happens this usually means that bad points have
1048 	 * been injected.
1049 	 *
1050 	 * NOTE: if for some reasons you need to deal with
1051 	 * all the possible pairs of points including these
1052 	 * exceptional pairs of inputs with an order 2 difference,
1053 	 * you should fallback to the incomplete formulas using the
1054 	 * COMPLETE=0 compilation toggle. Beware that in this
1055 	 * case, the library will be more sensitive to
1056 	 * side-channel attacks.
1057 	 */
1058 	ret = fp_iszero(&(out->Z), &cmp1); EG(ret, err);
1059 	ret = fp_iszero(&(out->Y), &cmp2); EG(ret, err);
1060 	MUST_HAVE(!((cmp1) && (cmp2)), ret, err);
1061 
1062 err:
1063 	fp_uninit(&t0);
1064 	fp_uninit(&t1);
1065 	fp_uninit(&t2);
1066 	fp_uninit(&t3);
1067 	fp_uninit(&t4);
1068 	fp_uninit(&t5);
1069 
1070 	return ret;
1071 }
1072 #endif  /* NO_USE_COMPLETE_FORMULAS */
1073 
1074 /*
1075  * Internal function:
1076  *
1077  *  - not supporting aliasing,
1078  *  - requiring caller to check in parameter is initialized
1079  *
1080  * Based on library configuration, the function either use complete formulas
1081  * or not.
1082  */
_prj_pt_dbl_monty(prj_pt_t out,prj_pt_src_t in)1083 static int _prj_pt_dbl_monty(prj_pt_t out, prj_pt_src_t in)
1084 {
1085 	int ret;
1086 
1087 #ifdef NO_USE_COMPLETE_FORMULAS
1088 	int iszero;
1089 	ret = prj_pt_iszero(in, &iszero); EG(ret, err);
1090 	if (iszero) {
1091 		ret = prj_pt_init(out, in->crv); EG(ret, err);
1092 		ret = prj_pt_zero(out);
1093 	} else {
1094 		ret = __prj_pt_dbl_monty_no_cf(out, in);
1095 	}
1096 #else
1097 	ret = __prj_pt_dbl_monty_cf(out, in); EG(ret, err);
1098 #endif
1099 
1100 err:
1101 	return ret;
1102 }
1103 
1104 /*
1105  * Internal version that peform in place doubling of given val,
1106  * by using a temporary copy. Sanity checks on parameters must
1107  * be done by caller.
1108  */
_prj_pt_dbl_monty_aliased(prj_pt_t val)1109 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_dbl_monty_aliased(prj_pt_t val)
1110 {
1111 	prj_pt out_cpy;
1112 	int ret;
1113 	out_cpy.magic = WORD(0);
1114 
1115 	ret = _prj_pt_dbl_monty(&out_cpy, val); EG(ret, err);
1116 	ret = prj_pt_copy(val, &out_cpy);
1117 
1118 err:
1119 	prj_pt_uninit(&out_cpy);
1120 
1121 	return ret;
1122 }
1123 
1124 /*
1125  * Public function for projective point doubling. The function handles the init
1126  * check of 'in' parameter which must be guaranteed for internal functions.
1127  * 'out' parameter need not be initialized and can be aliased with 'in'
1128  * parameter.
1129  *
1130  * The function returns 0 on success, -1 on error.
1131  */
prj_pt_dbl(prj_pt_t out,prj_pt_src_t in)1132 ATTRIBUTE_WARN_UNUSED_RET int prj_pt_dbl(prj_pt_t out, prj_pt_src_t in)
1133 {
1134 	int ret;
1135 
1136 	ret = prj_pt_check_initialized(in); EG(ret, err);
1137 
1138 	if (out == in) {
1139 		ret = _prj_pt_dbl_monty_aliased(out);
1140 	} else {
1141 		ret = _prj_pt_dbl_monty(out, in);
1142 	}
1143 
1144 err:
1145 	return ret;
1146 }
1147 
1148 /*
1149  * Internal function:
1150  *
1151  *  - not supporting aliasing,
1152  *  - requiring caller to check in1 and in2 parameter
1153  *
1154  * Based on library configuration, the function either use complete formulas
1155  * or not.
1156  */
_prj_pt_add_monty(prj_pt_t out,prj_pt_src_t in1,prj_pt_src_t in2)1157 ATTRIBUTE_WARN_UNUSED_RET static inline int _prj_pt_add_monty(prj_pt_t out,
1158 							      prj_pt_src_t in1,
1159 							      prj_pt_src_t in2)
1160 {
1161 #ifndef NO_USE_COMPLETE_FORMULAS
1162 	return __prj_pt_add_monty_cf(out, in1, in2);
1163 #else
1164 	return __prj_pt_add_monty_no_cf(out, in1, in2);
1165 #endif
1166 }
1167 
1168 /*
1169  * The function is an internal one that specifically handles aliasing. No check
1170  * is performed on parameters, this MUST be done by the caller:
1171  *
1172  *  - in1 and in2 are initialized
1173  *  - in1 and in2 are on the same curve
1174  *
1175  * The function will initialize 'out'. The function returns 0 on success, -1
1176  * on error.
1177  */
_prj_pt_add_monty_aliased(prj_pt_t out,prj_pt_src_t in1,prj_pt_src_t in2)1178 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_add_monty_aliased(prj_pt_t out,
1179 								prj_pt_src_t in1,
1180 								prj_pt_src_t in2)
1181 {
1182 	int ret;
1183 	prj_pt out_cpy;
1184 	out_cpy.magic = WORD(0);
1185 
1186 	ret = _prj_pt_add_monty(&out_cpy, in1, in2); EG(ret, err);
1187 	ret = prj_pt_copy(out, &out_cpy); EG(ret, err);
1188 
1189 err:
1190 	prj_pt_uninit(&out_cpy);
1191 
1192 	return ret;
1193 }
1194 
1195 /*
1196  * Public function for projective point addition. The function handles the
1197  * init checks of 'in1' and 'in2' parameters, along with the check they
1198  * use the same curve. This must be guaranteed for internal functions.
1199  * 'out' parameter need not be initialized and can be aliased with either
1200  * 'in1' or 'in2' parameter.
1201  *
1202  * The function returns 0 on success, -1 on error.
1203  */
prj_pt_add(prj_pt_t out,prj_pt_src_t in1,prj_pt_src_t in2)1204 int prj_pt_add(prj_pt_t out, prj_pt_src_t in1, prj_pt_src_t in2)
1205 {
1206 	int ret;
1207 
1208 	ret = prj_pt_check_initialized(in1); EG(ret, err);
1209 	ret = prj_pt_check_initialized(in2); EG(ret, err);
1210 	MUST_HAVE((in1->crv == in2->crv), ret, err);
1211 
1212 	if ((out == in1) || (out == in2)) {
1213 		ret = _prj_pt_add_monty_aliased(out, in1, in2);
1214 	} else {
1215 		ret = _prj_pt_add_monty(out, in1, in2);
1216 	}
1217 
1218 err:
1219 	return ret;
1220 }
1221 
1222 /*******************************************************************************/
1223 /****** Scalar multiplication algorithms ***************************************/
1224 /***********/
1225 /*
1226  * The description below summarizes the following algorithms.
1227  *
1228  * Double-and-Add-Always and Montgomery Ladder masked using Itoh et al. anti-ADPA
1229  * (Address-bit DPA) countermeasure.
1230  * See "A Practical Countermeasure against Address-Bit Differential Power Analysis"
1231  * by Itoh, Izu and Takenaka for more information.
1232  *
1233  * NOTE: these masked variants of the Double-and-Add-Always and Montgomery Ladder algorithms
1234  * are used by default as Itoh et al. countermeasure has a very small impact on performance
1235  * and is inherently more robust againt DPA. The only case where we use another variant is
1236  * for devices with low memory as Itoh requires many temporary variables that consume many
1237  * temporary stack space.
1238  *
1239  * NOTE: the algorithms inherently depend on the MSB of the
1240  * scalar. In order to avoid leaking this MSB and fall into HNP (Hidden Number
1241  * Problem) issues, we use the trick described in https://eprint.iacr.org/2011/232.pdf
1242  * to have the MSB always set. However, since the scalar m might be less or bigger than
1243  * the order q of the curve, we distinguish three situations:
1244  *     - The scalar m is < q (the order), in this case we compute:
1245  *         -
1246  *        | m' = m + (2 * q) if [log(k + q)] == [log(q)],
1247  *        | m' = m + q otherwise.
1248  *         -
1249  *     - The scalar m is >= q and < q**2, in this case we compute:
1250  *         -
1251  *        | m' = m + (2 * (q**2)) if [log(k + (q**2))] == [log(q**2)],
1252  *        | m' = m + (q**2) otherwise.
1253  *         -
1254  *     - The scalar m is >= (q**2), in this case m == m'
1255  *
1256  *   => We only deal with 0 <= m < (q**2) using the countermeasure. When m >= (q**2),
1257  *      we stick with m' = m, accepting MSB issues (not much can be done in this case
1258  *      anyways). In the two first cases, Double-and-Add-Always and Montgomery Ladder are
1259  *      performed in constant time wrt the size of the scalar m.
1260  */
1261 /***********/
1262 /*
1263  * Internal point blinding function: as it is internal, in is supposed to be initialized and
1264  * aliasing is NOT supported.
1265  */
_blind_projective_point(prj_pt_t out,prj_pt_src_t in)1266 ATTRIBUTE_WARN_UNUSED_RET static int _blind_projective_point(prj_pt_t out, prj_pt_src_t in)
1267 {
1268 	int ret;
1269 
1270 	/* Random for projective coordinates masking */
1271 	/* NOTE: to limit stack usage, we reuse out->Z as a temporary
1272 	 * variable. This does not work if in == out, hence the check.
1273 	 */
1274 	MUST_HAVE((in != out), ret, err);
1275 
1276 	ret = prj_pt_init(out, in->crv); EG(ret, err);
1277 
1278 	/* Get a random value l in Fp */
1279 	ret = fp_get_random(&(out->Z), in->X.ctx); EG(ret, err);
1280 
1281 	/*
1282 	 * Blind the point with projective coordinates
1283 	 * (X, Y, Z) => (l*X, l*Y, l*Z)
1284 	 */
1285 	ret = fp_mul_monty(&(out->X), &(in->X), &(out->Z)); EG(ret, err);
1286 	ret = fp_mul_monty(&(out->Y), &(in->Y), &(out->Z)); EG(ret, err);
1287 	ret = fp_mul_monty(&(out->Z), &(in->Z), &(out->Z));
1288 
1289 err:
1290 	return ret;
1291 }
1292 
1293 /* If nothing is specified regarding the scalar multiplication algorithm, we use
1294  * the Montgomery Ladder. For the specific case of small stack devices, we release
1295  * some pressure on the stack by explicitly using double and always WITHOUT the Itoh
1296  * et al. countermeasure against A-DPA as it is quite consuming.
1297  */
1298 #if defined(USE_SMALL_STACK) && defined(USE_MONTY_LADDER)
1299 #error "Small stack is only compatible with USE_DOUBLE_ADD_ALWAYS while USE_MONTY_LADDER has been explicitly asked!"
1300 #endif
1301 
1302 #if defined(USE_SMALL_STACK)
1303 #ifndef USE_DOUBLE_ADD_ALWAYS
1304 #define USE_DOUBLE_ADD_ALWAYS
1305 #endif
1306 #endif
1307 
1308 #if !defined(USE_DOUBLE_ADD_ALWAYS) && !defined(USE_MONTY_LADDER)
1309 #define USE_MONTY_LADDER
1310 #endif
1311 
1312 #if defined(USE_DOUBLE_ADD_ALWAYS) && defined(USE_MONTY_LADDER)
1313 #error "You can either choose USE_DOUBLE_ADD_ALWAYS or USE_MONTY_LADDER, not both!"
1314 #endif
1315 
1316 #if defined(USE_DOUBLE_ADD_ALWAYS) && !defined(USE_SMALL_STACK)
_prj_pt_mul_ltr_monty_dbl_add_always(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1317 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_mul_ltr_monty_dbl_add_always(prj_pt_t out, nn_src_t m, prj_pt_src_t in)
1318 {
1319 	/* We use Itoh et al. notations here for T and the random r */
1320 	prj_pt T[3];
1321 	bitcnt_t mlen;
1322 	u8 mbit, rbit;
1323 	/* Random for masking the Double and Add Always algorithm */
1324 	nn r;
1325 	/* The new scalar we will use with MSB fixed to 1 (noted m' above).
1326 	 * This helps dealing with constant time.
1327 	 */
1328 	nn m_msb_fixed;
1329 	nn_src_t curve_order;
1330 	nn curve_order_square;
1331 	int ret, ret_ops, cmp;
1332 	r.magic = m_msb_fixed.magic = curve_order_square.magic = WORD(0);
1333 	T[0].magic = T[1].magic = T[2].magic = WORD(0);
1334 
1335 	/* Compute m' from m depending on the rule described above */
1336 	curve_order = &(in->crv->order);
1337 	/* First compute q**2 */
1338 	ret = nn_sqr(&curve_order_square, curve_order); EG(ret, err);
1339 	/* Then compute m' depending on m size */
1340 	ret = nn_cmp(m, curve_order, &cmp); EG(ret, err);
1341 	if (cmp < 0){
1342 		bitcnt_t msb_bit_len, order_bitlen;
1343 
1344 		/* Case where m < q */
1345 		ret = nn_add(&m_msb_fixed, m, curve_order); EG(ret, err);
1346 		ret = nn_bitlen(&m_msb_fixed, &msb_bit_len); EG(ret, err);
1347 		ret = nn_bitlen(curve_order, &order_bitlen); EG(ret, err);
1348 		ret = nn_cnd_add((msb_bit_len == order_bitlen), &m_msb_fixed,
1349 				  &m_msb_fixed, curve_order); EG(ret, err);
1350 	} else {
1351 		ret = nn_cmp(m, &curve_order_square, &cmp); EG(ret, err);
1352 		if (cmp < 0) {
1353 			bitcnt_t msb_bit_len, curve_order_square_bitlen;
1354 
1355 			/* Case where m >= q and m < (q**2) */
1356 			ret = nn_add(&m_msb_fixed, m, &curve_order_square); EG(ret, err);
1357 			ret = nn_bitlen(&m_msb_fixed, &msb_bit_len); EG(ret, err);
1358 			ret = nn_bitlen(&curve_order_square, &curve_order_square_bitlen); EG(ret, err);
1359 			ret = nn_cnd_add((msb_bit_len == curve_order_square_bitlen),
1360 					&m_msb_fixed, &m_msb_fixed, &curve_order_square); EG(ret, err);
1361 		} else {
1362 			/* Case where m >= (q**2) */
1363 			ret = nn_copy(&m_msb_fixed, m); EG(ret, err);
1364 		}
1365 	}
1366 	ret = nn_bitlen(&m_msb_fixed, &mlen); EG(ret, err);
1367 	MUST_HAVE(mlen != 0, ret, err);
1368 	mlen--;
1369 
1370 	/* Hide possible internal failures for double and add
1371 	 * operations and perform the operation in constant time.
1372 	 */
1373 	ret_ops = 0;
1374 
1375 	/* Get a random r with the same size of m_msb_fixed */
1376 	ret = nn_get_random_len(&r, m_msb_fixed.wlen * WORD_BYTES); EG(ret, err);
1377 
1378 	ret = nn_getbit(&r, mlen, &rbit); EG(ret, err);
1379 
1380 	/* Initialize points */
1381 	ret = prj_pt_init(&T[0], in->crv); EG(ret, err);
1382 	ret = prj_pt_init(&T[1], in->crv); EG(ret, err);
1383 
1384 	/*
1385 	 * T[2] = R(P)
1386 	 * Blind the point with projective coordinates
1387 	 * (X, Y, Z) => (l*X, l*Y, l*Z)
1388 	 */
1389 	ret = _blind_projective_point(&T[2], in); EG(ret, err);
1390 
1391 	/*  T[r[n-1]] = T[2] */
1392 	ret = prj_pt_copy(&T[rbit], &T[2]); EG(ret, err);
1393 
1394 	/* Main loop of Double and Add Always */
1395 	while (mlen > 0) {
1396 		u8 rbit_next;
1397 		--mlen;
1398 		/* rbit is r[i+1], and rbit_next is r[i] */
1399 		ret = nn_getbit(&r, mlen, &rbit_next); EG(ret, err);
1400 
1401 		/* mbit is m[i] */
1402 		ret = nn_getbit(&m_msb_fixed, mlen, &mbit); EG(ret, err);
1403 
1404 		/* Double: T[r[i+1]] = ECDBL(T[r[i+1]]) */
1405 #ifndef NO_USE_COMPLETE_FORMULAS
1406 		/*
1407 		 * NOTE: in case of complete formulas, we use the
1408 		 * addition for doubling, incurring a small performance hit
1409 		 * for better SCA resistance.
1410 		 */
1411 		ret_ops |= prj_pt_add(&T[rbit], &T[rbit], &T[rbit]);
1412 #else
1413 		ret_ops |= prj_pt_dbl(&T[rbit], &T[rbit]);
1414 #endif
1415 		/* Add:  T[1-r[i+1]] = ECADD(T[r[i+1]],T[2]) */
1416 		ret_ops |= prj_pt_add(&T[1-rbit], &T[rbit], &T[2]);
1417 
1418 		/*
1419 		 * T[r[i]] = T[d[i] ^ r[i+1]]
1420 		 * NOTE: we use the low level nn_copy function here to avoid
1421 		 * any possible leakage on operands with prj_pt_copy
1422 		 */
1423 		ret = nn_copy(&(T[rbit_next].X.fp_val), &(T[mbit ^ rbit].X.fp_val)); EG(ret, err);
1424 		ret = nn_copy(&(T[rbit_next].Y.fp_val), &(T[mbit ^ rbit].Y.fp_val)); EG(ret, err);
1425 		ret = nn_copy(&(T[rbit_next].Z.fp_val), &(T[mbit ^ rbit].Z.fp_val)); EG(ret, err);
1426 
1427 		/* Update rbit */
1428 		rbit = rbit_next;
1429 	}
1430 	/* Output: T[r[0]] */
1431 	ret = prj_pt_copy(out, &T[rbit]); EG(ret, err);
1432 
1433 	/* Take into consideration our double and add errors */
1434 	ret |= ret_ops;
1435 
1436 err:
1437 	prj_pt_uninit(&T[0]);
1438 	prj_pt_uninit(&T[1]);
1439 	prj_pt_uninit(&T[2]);
1440 	nn_uninit(&r);
1441 	nn_uninit(&m_msb_fixed);
1442 	nn_uninit(&curve_order_square);
1443 
1444 	PTR_NULLIFY(curve_order);
1445 
1446 	return ret;
1447 }
1448 #endif
1449 
1450 #if defined(USE_DOUBLE_ADD_ALWAYS) && defined(USE_SMALL_STACK)
1451 /* NOTE: in small stack case where we compile for low memory devices, we do not use Itoh et al. countermeasure
1452  * as it requires too much temporary space on the stack.
1453  */
_prj_pt_mul_ltr_monty_dbl_add_always(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1454 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_mul_ltr_monty_dbl_add_always(prj_pt_t out, nn_src_t m, prj_pt_src_t in)
1455 {
1456 	int ret, ret_ops;
1457 
1458 	/* Hide possible internal failures for double and add
1459 	 * operations and perform the operation in constant time.
1460 	 */
1461 	ret_ops = 0;
1462 
1463 	/* Blind the input point projective coordinates */
1464 	ret = _blind_projective_point(out, in); EG(ret, err);
1465 
1466 	/*******************/
1467 	{
1468 		bitcnt_t mlen;
1469 		u8 mbit;
1470 		/* The new scalar we will use with MSB fixed to 1 (noted m' above).
1471 		 * This helps dealing with constant time.
1472 		 */
1473 		nn m_msb_fixed;
1474 		nn_src_t curve_order;
1475 		int cmp;
1476 		m_msb_fixed.magic = WORD(0);
1477 
1478 		{
1479 			nn curve_order_square;
1480 			curve_order_square.magic = WORD(0);
1481 
1482 			/* Compute m' from m depending on the rule described above */
1483 			curve_order = &(in->crv->order);
1484 			/* First compute q**2 */
1485 			ret = nn_sqr(&curve_order_square, curve_order); EG(ret, err1);
1486 			/* Then compute m' depending on m size */
1487 			ret = nn_cmp(m, curve_order, &cmp); EG(ret, err1);
1488 			if (cmp < 0){
1489 				bitcnt_t msb_bit_len, order_bitlen;
1490 
1491 				/* Case where m < q */
1492 				ret = nn_add(&m_msb_fixed, m, curve_order); EG(ret, err1);
1493 				ret = nn_bitlen(&m_msb_fixed, &msb_bit_len); EG(ret, err1);
1494 				ret = nn_bitlen(curve_order, &order_bitlen); EG(ret, err1);
1495 				ret = nn_cnd_add((msb_bit_len == order_bitlen), &m_msb_fixed,
1496 					  &m_msb_fixed, curve_order); EG(ret, err1);
1497 			} else {
1498 				ret = nn_cmp(m, &curve_order_square, &cmp); EG(ret, err1);
1499 				if (cmp < 0) {
1500 					bitcnt_t msb_bit_len, curve_order_square_bitlen;
1501 
1502 					/* Case where m >= q and m < (q**2) */
1503 					ret = nn_add(&m_msb_fixed, m, &curve_order_square); EG(ret, err1);
1504 					ret = nn_bitlen(&m_msb_fixed, &msb_bit_len); EG(ret, err1);
1505 					ret = nn_bitlen(&curve_order_square, &curve_order_square_bitlen); EG(ret, err1);
1506 					ret = nn_cnd_add((msb_bit_len == curve_order_square_bitlen),
1507 							&m_msb_fixed, &m_msb_fixed, &curve_order_square); EG(ret, err1);
1508 				} else {
1509 					/* Case where m >= (q**2) */
1510 					ret = nn_copy(&m_msb_fixed, m); EG(ret, err1);
1511 				}
1512 			}
1513 err1:
1514 			nn_uninit(&curve_order_square); EG(ret, err);
1515 		}
1516 
1517 		ret = nn_bitlen(&m_msb_fixed, &mlen); EG(ret, err);
1518 		MUST_HAVE((mlen != 0), ret, err);
1519 		mlen--;
1520 
1521 		{
1522 			prj_pt dbl;
1523 			dbl.magic = WORD(0);
1524 
1525 			/* Initialize temporary point */
1526 			ret = prj_pt_init(&dbl, in->crv); EG(ret, err2);
1527 
1528 			/* Main loop of Double and Add Always */
1529 			while (mlen > 0) {
1530 				--mlen;
1531 				/* mbit is m[i] */
1532 				ret = nn_getbit(&m_msb_fixed, mlen, &mbit); EG(ret, err2);
1533 
1534 #ifndef NO_USE_COMPLETE_FORMULAS
1535 				/*
1536 				 * NOTE: in case of complete formulas, we use the
1537 				 * addition for doubling, incurring a small performance hit
1538 				 * for better SCA resistance.
1539 				 */
1540 				ret_ops |= prj_pt_add(&dbl, out, out);
1541 #else
1542 				ret_ops |= prj_pt_dbl(&dbl, out);
1543 #endif
1544 				ret_ops |= prj_pt_add(out, &dbl, in);
1545 				/* Swap */
1546 				ret = nn_cnd_swap(!mbit, &(out->X.fp_val), &(dbl.X.fp_val)); EG(ret, err2);
1547 				ret = nn_cnd_swap(!mbit, &(out->Y.fp_val), &(dbl.Y.fp_val)); EG(ret, err2);
1548 				ret = nn_cnd_swap(!mbit, &(out->Z.fp_val), &(dbl.Z.fp_val)); EG(ret, err2);
1549 			}
1550 err2:
1551 			prj_pt_uninit(&dbl); EG(ret, err);
1552 		}
1553 
1554 err:
1555 		nn_uninit(&m_msb_fixed);
1556 
1557 		PTR_NULLIFY(curve_order);
1558 	}
1559 
1560 	/* Take into consideration our double and add errors */
1561 	ret |= ret_ops;
1562 
1563 	return ret;
1564 }
1565 #endif
1566 
1567 
1568 #ifdef USE_MONTY_LADDER
_prj_pt_mul_ltr_monty_ladder(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1569 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_mul_ltr_monty_ladder(prj_pt_t out, nn_src_t m, prj_pt_src_t in)
1570 {
1571 	/* We use Itoh et al. notations here for T and the random r */
1572 	prj_pt T[3];
1573 	bitcnt_t mlen;
1574 	u8 mbit, rbit;
1575 	/* Random for masking the Montgomery Ladder algorithm */
1576 	nn r;
1577 	/* The new scalar we will use with MSB fixed to 1 (noted m' above).
1578 	 * This helps dealing with constant time.
1579 	 */
1580 	nn m_msb_fixed;
1581 	nn_src_t curve_order;
1582 	nn curve_order_square;
1583 	int ret, ret_ops, cmp;
1584 	r.magic = m_msb_fixed.magic = curve_order_square.magic = WORD(0);
1585 	T[0].magic = T[1].magic = T[2].magic = WORD(0);
1586 
1587 	/* Compute m' from m depending on the rule described above */
1588 	curve_order = &(in->crv->order);
1589 
1590 	/* First compute q**2 */
1591 	ret = nn_sqr(&curve_order_square, curve_order); EG(ret, err);
1592 
1593 	/* Then compute m' depending on m size */
1594 	ret = nn_cmp(m, curve_order, &cmp); EG(ret, err);
1595 	if (cmp < 0) {
1596 		bitcnt_t msb_bit_len, order_bitlen;
1597 
1598 		/* Case where m < q */
1599 		ret = nn_add(&m_msb_fixed, m, curve_order); EG(ret, err);
1600 		ret = nn_bitlen(&m_msb_fixed, &msb_bit_len); EG(ret, err);
1601 		ret = nn_bitlen(curve_order, &order_bitlen); EG(ret, err);
1602 		ret = nn_cnd_add((msb_bit_len == order_bitlen), &m_msb_fixed,
1603 				&m_msb_fixed, curve_order); EG(ret, err);
1604 	} else {
1605 		ret = nn_cmp(m, &curve_order_square, &cmp); EG(ret, err);
1606 		if (cmp < 0) {
1607 			bitcnt_t msb_bit_len, curve_order_square_bitlen;
1608 
1609 			/* Case where m >= q and m < (q**2) */
1610 			ret = nn_add(&m_msb_fixed, m, &curve_order_square); EG(ret, err);
1611 			ret = nn_bitlen(&m_msb_fixed, &msb_bit_len); EG(ret, err);
1612 			ret = nn_bitlen(&curve_order_square, &curve_order_square_bitlen); EG(ret, err);
1613 			ret = nn_cnd_add((msb_bit_len == curve_order_square_bitlen),
1614 					 &m_msb_fixed, &m_msb_fixed, &curve_order_square); EG(ret, err);
1615 		} else {
1616 			/* Case where m >= (q**2) */
1617 			ret = nn_copy(&m_msb_fixed, m); EG(ret, err);
1618 		}
1619 	}
1620 
1621 	ret = nn_bitlen(&m_msb_fixed, &mlen); EG(ret, err);
1622 	MUST_HAVE((mlen != 0), ret, err);
1623 	mlen--;
1624 
1625 	/* Hide possible internal failures for double and add
1626 	 * operations and perform the operation in constant time.
1627 	 */
1628 	ret_ops = 0;
1629 
1630 	/* Get a random r with the same size of m_msb_fixed */
1631 	ret = nn_get_random_len(&r, (u16)(m_msb_fixed.wlen * WORD_BYTES)); EG(ret, err);
1632 
1633 	ret = nn_getbit(&r, mlen, &rbit); EG(ret, err);
1634 
1635 	/* Initialize points */
1636 	ret = prj_pt_init(&T[0], in->crv); EG(ret, err);
1637 	ret = prj_pt_init(&T[1], in->crv); EG(ret, err);
1638 	ret = prj_pt_init(&T[2], in->crv); EG(ret, err);
1639 
1640 	/* Initialize T[r[n-1]] to input point */
1641 	/*
1642 	 * Blind the point with projective coordinates
1643 	 * (X, Y, Z) => (l*X, l*Y, l*Z)
1644 	 */
1645 	ret = _blind_projective_point(&T[rbit], in); EG(ret, err);
1646 
1647 	/* Initialize T[1-r[n-1]] with ECDBL(T[r[n-1]])) */
1648 #ifndef NO_USE_COMPLETE_FORMULAS
1649 	/*
1650 	 * NOTE: in case of complete formulas, we use the
1651 	 * addition for doubling, incurring a small performance hit
1652 	 * for better SCA resistance.
1653 	 */
1654 	ret_ops |= prj_pt_add(&T[1-rbit], &T[rbit], &T[rbit]);
1655 #else
1656 	ret_ops |= prj_pt_dbl(&T[1-rbit], &T[rbit]);
1657 #endif
1658 
1659 	/* Main loop of the Montgomery Ladder */
1660 	while (mlen > 0) {
1661 		u8 rbit_next;
1662 		--mlen;
1663 		/* rbit is r[i+1], and rbit_next is r[i] */
1664 		ret = nn_getbit(&r, mlen, &rbit_next); EG(ret, err);
1665 
1666 		/* mbit is m[i] */
1667 		ret = nn_getbit(&m_msb_fixed, mlen, &mbit); EG(ret, err);
1668 		/* Double: T[2] = ECDBL(T[d[i] ^ r[i+1]]) */
1669 
1670 #ifndef NO_USE_COMPLETE_FORMULAS
1671 		/* NOTE: in case of complete formulas, we use the
1672 		 * addition for doubling, incurring a small performance hit
1673 		 * for better SCA resistance.
1674 		 */
1675 		ret_ops |= prj_pt_add(&T[2], &T[mbit ^ rbit], &T[mbit ^ rbit]);
1676 #else
1677 		ret_ops |= prj_pt_dbl(&T[2], &T[mbit ^ rbit]);
1678 #endif
1679 
1680 		/* Add: T[1] = ECADD(T[0],T[1]) */
1681 		ret_ops |= prj_pt_add(&T[1], &T[0], &T[1]);
1682 
1683 		/* T[0] = T[2-(d[i] ^ r[i])] */
1684 		/*
1685 		 * NOTE: we use the low level nn_copy function here to avoid
1686 		 * any possible leakage on operands with prj_pt_copy
1687 		 */
1688 		ret = nn_copy(&(T[0].X.fp_val), &(T[2-(mbit ^ rbit_next)].X.fp_val)); EG(ret, err);
1689 		ret = nn_copy(&(T[0].Y.fp_val), &(T[2-(mbit ^ rbit_next)].Y.fp_val)); EG(ret, err);
1690 		ret = nn_copy(&(T[0].Z.fp_val), &(T[2-(mbit ^ rbit_next)].Z.fp_val)); EG(ret, err);
1691 
1692 		/* T[1] = T[1+(d[i] ^ r[i])] */
1693 		/* NOTE: we use the low level nn_copy function here to avoid
1694 		 * any possible leakage on operands with prj_pt_copy
1695 		 */
1696 		ret = nn_copy(&(T[1].X.fp_val), &(T[1+(mbit ^ rbit_next)].X.fp_val)); EG(ret, err);
1697 		ret = nn_copy(&(T[1].Y.fp_val), &(T[1+(mbit ^ rbit_next)].Y.fp_val)); EG(ret, err);
1698 		ret = nn_copy(&(T[1].Z.fp_val), &(T[1+(mbit ^ rbit_next)].Z.fp_val)); EG(ret, err);
1699 
1700 		/* Update rbit */
1701 		rbit = rbit_next;
1702 	}
1703 	/* Output: T[r[0]] */
1704 	ret = prj_pt_copy(out, &T[rbit]); EG(ret, err);
1705 
1706 	/* Take into consideration our double and add errors */
1707 	ret |= ret_ops;
1708 
1709 err:
1710 	prj_pt_uninit(&T[0]);
1711 	prj_pt_uninit(&T[1]);
1712 	prj_pt_uninit(&T[2]);
1713 	nn_uninit(&r);
1714 	nn_uninit(&m_msb_fixed);
1715 	nn_uninit(&curve_order_square);
1716 
1717 	PTR_NULLIFY(curve_order);
1718 
1719 	return ret;
1720 }
1721 #endif
1722 
1723 /* Main projective scalar multiplication function.
1724  * Depending on the preprocessing options, we use either the
1725  * Double and Add Always algorithm, or the Montgomery Ladder one.
1726  */
_prj_pt_mul_ltr_monty(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1727 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_mul_ltr_monty(prj_pt_t out, nn_src_t m, prj_pt_src_t in){
1728 #if defined(USE_DOUBLE_ADD_ALWAYS)
1729 	return _prj_pt_mul_ltr_monty_dbl_add_always(out, m, in);
1730 #elif defined(USE_MONTY_LADDER)
1731 	return _prj_pt_mul_ltr_monty_ladder(out, m, in);
1732 #else
1733 #error "Error: neither Double and Add Always nor Montgomery Ladder has been selected!"
1734 #endif
1735 }
1736 
1737 /* version with 'm' passed via 'out'. */
_prj_pt_mul_ltr_monty_aliased(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1738 ATTRIBUTE_WARN_UNUSED_RET static int _prj_pt_mul_ltr_monty_aliased(prj_pt_t out, nn_src_t m, prj_pt_src_t in)
1739 {
1740 	prj_pt out_cpy;
1741 	int ret;
1742 	out_cpy.magic = WORD(0);
1743 
1744 	ret = prj_pt_init(&out_cpy, in->crv); EG(ret, err);
1745 	ret = _prj_pt_mul_ltr_monty(&out_cpy, m, in); EG(ret, err);
1746 	ret = prj_pt_copy(out, &out_cpy);
1747 
1748 err:
1749 	prj_pt_uninit(&out_cpy);
1750 	return ret;
1751 }
1752 
1753 /* Aliased version. This is the public main interface of our
1754  * scalar multiplication algorithm. Checks that the input point
1755  * and that the output point are on the curve are performed here
1756  * (before and after calling the core algorithm, albeit Double and
1757  * Add Always or Montgomery Ladder).
1758  */
prj_pt_mul(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1759 int prj_pt_mul(prj_pt_t out, nn_src_t m, prj_pt_src_t in)
1760 {
1761 	int ret, on_curve;
1762 
1763 	ret = prj_pt_check_initialized(in); EG(ret, err);
1764 	ret = nn_check_initialized(m); EG(ret, err);
1765 
1766 	/* Check that the input is on the curve */
1767 	MUST_HAVE((!prj_pt_is_on_curve(in, &on_curve)) && on_curve, ret, err);
1768 
1769 	if (out == in) {
1770 		ret = _prj_pt_mul_ltr_monty_aliased(out, m, in); EG(ret, err);
1771 	} else {
1772 		ret = _prj_pt_mul_ltr_monty(out, m, in); EG(ret, err);
1773 	}
1774 
1775 	/* Check that the output is on the curve */
1776 	MUST_HAVE((!prj_pt_is_on_curve(out, &on_curve)) && on_curve, ret, err);
1777 
1778 err:
1779 	return ret;
1780 }
1781 
prj_pt_mul_blind(prj_pt_t out,nn_src_t m,prj_pt_src_t in)1782 int prj_pt_mul_blind(prj_pt_t out, nn_src_t m, prj_pt_src_t in)
1783 {
1784 	/* Blind the scalar m with (b*q) where q is the curve order.
1785 	 * NOTE: the curve order and the "generator" order are
1786 	 * usually the same (i.e. cofactor = 1) for the classical
1787 	 * prime fields curves. However some exceptions exist
1788 	 * (e.g. Wei25519 and Wei448), and in this case it is
1789 	 * curcial to use the curve order for a generic blinding
1790 	 * working on any point on the curve.
1791 	 */
1792 	nn b;
1793 	nn_src_t q;
1794 	int ret;
1795 	b.magic = WORD(0);
1796 
1797 	ret = prj_pt_check_initialized(in); EG(ret, err);
1798 
1799 	q = &(in->crv->order);
1800 
1801 	ret = nn_init(&b, 0); EG(ret, err);
1802 
1803 	ret = nn_get_random_mod(&b, q); EG(ret, err);
1804 
1805 	ret = nn_mul(&b, &b, q); EG(ret, err);
1806 	ret = nn_add(&b, &b, m); EG(ret, err);
1807 
1808 	/* NOTE: point blinding is performed in the lower functions */
1809 	/* NOTE: check that input and output points are on the curve are
1810 	 * performed in the lower functions.
1811 	 */
1812 
1813 	/* Perform the scalar multiplication */
1814 	ret = prj_pt_mul(out, &b, in);
1815 
1816 err:
1817 	nn_uninit(&b);
1818 
1819 	PTR_NULLIFY(q);
1820 
1821 	return ret;
1822 }
1823 
1824 /* Naive double and add scalar multiplication.
1825  *
1826  * This scalar multiplication is used on public values and is optimized with no
1827  * countermeasures, and it is usually faster as scalar can be small with few bits
1828  * to process (e.g. cofactors, etc.).
1829  *
1830  * out is initialized by the function.
1831  *
1832  * XXX: WARNING: this function must only be used on public points!
1833  *
1834  */
__prj_pt_unprotected_mult(prj_pt_t out,nn_src_t scalar,prj_pt_src_t public_in)1835 static int __prj_pt_unprotected_mult(prj_pt_t out, nn_src_t scalar, prj_pt_src_t public_in)
1836 {
1837         u8 expbit;
1838         bitcnt_t explen;
1839         int ret, iszero, on_curve;
1840 
1841         ret = prj_pt_check_initialized(public_in); EG(ret, err);
1842         ret = nn_check_initialized(scalar); EG(ret, err);
1843 
1844 	/* This function does not support aliasing */
1845 	MUST_HAVE((out != public_in), ret, err);
1846 
1847 	/* Check that the input is on the curve */
1848 	MUST_HAVE((!prj_pt_is_on_curve(public_in, &on_curve)) && on_curve, ret, err);
1849 
1850         ret = nn_iszero(scalar, &iszero); EG(ret, err);
1851 	/* Multiplication by zero is the point at infinity */
1852 	if(iszero){
1853 		ret = prj_pt_zero(out); EG(ret, err);
1854 		goto err;
1855 	}
1856 
1857         ret = nn_bitlen(scalar, &explen); EG(ret, err);
1858         /* Sanity check */
1859         MUST_HAVE((explen > 0), ret, err);
1860         explen = (bitcnt_t)(explen - 1);
1861 	ret = prj_pt_copy(out, public_in); EG(ret, err);
1862 
1863         while (explen > 0) {
1864                 explen = (bitcnt_t)(explen - 1);
1865                 ret = nn_getbit(scalar, explen, &expbit); EG(ret, err);
1866                 ret = prj_pt_dbl(out, out); EG(ret, err);
1867                 if(expbit){
1868                         ret = prj_pt_add(out, out, public_in); EG(ret, err);
1869                 }
1870         }
1871 
1872 	/* Check that the output is on the curve */
1873 	MUST_HAVE((!prj_pt_is_on_curve(out, &on_curve)) && on_curve, ret, err);
1874 
1875 err:
1876         VAR_ZEROIFY(expbit);
1877         VAR_ZEROIFY(explen);
1878 
1879         return ret;
1880 }
1881 
1882 /* Aliased version of __prj_pt_unprotected_mult */
_prj_pt_unprotected_mult(prj_pt_t out,nn_src_t scalar,prj_pt_src_t public_in)1883 int _prj_pt_unprotected_mult(prj_pt_t out, nn_src_t scalar, prj_pt_src_t public_in)
1884 {
1885 	int ret;
1886 
1887 	if(out == public_in){
1888                 prj_pt A;
1889                 A.magic = WORD(0);
1890 
1891                 ret = prj_pt_copy(&A, public_in); EG(ret, err1);
1892 		ret = __prj_pt_unprotected_mult(out, scalar, &A);
1893 err1:
1894 		prj_pt_uninit(&A);
1895 		goto err;
1896 	}
1897 	else{
1898 		ret = __prj_pt_unprotected_mult(out, scalar, public_in);
1899 	}
1900 err:
1901 	return ret;
1902 }
1903 /*
1904  * Check if an integer is (a multiple of) a projective point order.
1905  *
1906  * The function returns 0 on success, -1 on error. The value check is set to 1 if the projective
1907  * point has order in_isorder, 0 otherwise. The value is meaningless on error.
1908  */
check_prj_pt_order(prj_pt_src_t in_shortw,nn_src_t in_isorder,prj_pt_sensitivity s,int * check)1909 int check_prj_pt_order(prj_pt_src_t in_shortw, nn_src_t in_isorder, prj_pt_sensitivity s, int *check)
1910 {
1911 	int ret, iszero;
1912 	prj_pt res;
1913 	res.magic = WORD(0);
1914 
1915 	/* First sanity checks */
1916 	ret = prj_pt_check_initialized(in_shortw); EG(ret, err);
1917 	ret = nn_check_initialized(in_isorder); EG(ret, err);
1918 	MUST_HAVE((check != NULL), ret, err);
1919 
1920 	/* Then, perform the scalar multiplication */
1921 	if(s == PUBLIC_PT){
1922 		/* If this is a public point, we can use the naive scalar multiplication */
1923 		ret = _prj_pt_unprotected_mult(&res, in_isorder, in_shortw); EG(ret, err);
1924 	}
1925 	else{
1926 		/* If the point is private, it is sensitive and we proceed with the secure
1927 		 * scalar blind multiplication.
1928 		 */
1929 		ret = prj_pt_mul_blind(&res, in_isorder, in_shortw); EG(ret, err);
1930 	}
1931 
1932 	/* Check if we have the point at infinity */
1933 	ret = prj_pt_iszero(&res, &iszero); EG(ret, err);
1934 	(*check) = iszero;
1935 
1936 err:
1937 	prj_pt_uninit(&res);
1938 
1939 	return ret;
1940 }
1941 
1942 /*****************************************************************************/
1943 
1944 /*
1945  * Map points from Edwards to short Weierstrass projective points through Montgomery (composition mapping).
1946  *     Point at infinity (0, 1) -> (0, 1, 0) is treated as an exception, which is trivially not constant time.
1947  *     This is OK since our mapping functions should be used at the non sensitive input and output
1948  *     interfaces.
1949  *
1950  * The function returns 0 on success, -1 on error.
1951  */
aff_pt_edwards_to_prj_pt_shortw(aff_pt_edwards_src_t in_edwards,ec_shortw_crv_src_t shortw_crv,prj_pt_t out_shortw,fp_src_t alpha_edwards)1952 int aff_pt_edwards_to_prj_pt_shortw(aff_pt_edwards_src_t in_edwards,
1953 				    ec_shortw_crv_src_t shortw_crv,
1954 				    prj_pt_t out_shortw,
1955 				    fp_src_t alpha_edwards)
1956 {
1957 	int ret, iszero, cmp;
1958 	aff_pt out_shortw_aff;
1959 	fp one;
1960 	out_shortw_aff.magic = one.magic = WORD(0);
1961 
1962 	/* Check the curves compatibility */
1963 	ret = aff_pt_edwards_check_initialized(in_edwards); EG(ret, err);
1964 	ret = curve_edwards_shortw_check(in_edwards->crv, shortw_crv, alpha_edwards); EG(ret, err);
1965 
1966 	/* Initialize output point with curve */
1967 	ret = prj_pt_init(out_shortw, shortw_crv); EG(ret, err);
1968 
1969 	ret = fp_init(&one, in_edwards->x.ctx); EG(ret, err);
1970 	ret = fp_one(&one); EG(ret, err);
1971 
1972 	/* Check if we are the point at infinity
1973 	 * This check induces a non consant time exception, but the current function must be called on
1974 	 * public data anyways.
1975 	 */
1976 	ret = fp_iszero(&(in_edwards->x), &iszero); EG(ret, err);
1977 	ret = fp_cmp(&(in_edwards->y), &one, &cmp); EG(ret, err);
1978 	if(iszero && (cmp == 0)){
1979 		ret = prj_pt_zero(out_shortw); EG(ret, err);
1980 		ret = 0;
1981 		goto err;
1982 	}
1983 
1984 	/* Use the affine mapping */
1985 	ret = aff_pt_edwards_to_shortw(in_edwards, shortw_crv, &out_shortw_aff, alpha_edwards); EG(ret, err);
1986 	/* And then map the short Weierstrass affine to projective coordinates */
1987 	ret = ec_shortw_aff_to_prj(out_shortw, &out_shortw_aff);
1988 
1989 err:
1990 	fp_uninit(&one);
1991 	aff_pt_uninit(&out_shortw_aff);
1992 
1993 	return ret;
1994 }
1995 
1996 /*
1997  * Map points from short Weierstrass projective points to Edwards through Montgomery (composition mapping).
1998  *     Point at infinity with Z=0 (in projective coordinates) -> (0, 1) is treated as an exception, which is trivially not constant time.
1999  *     This is OK since our mapping functions should be used at the non sensitive input and output
2000  *     interfaces.
2001  *
2002  * The function returns 0 on success, -1 on error.
2003  */
prj_pt_shortw_to_aff_pt_edwards(prj_pt_src_t in_shortw,ec_edwards_crv_src_t edwards_crv,aff_pt_edwards_t out_edwards,fp_src_t alpha_edwards)2004 int prj_pt_shortw_to_aff_pt_edwards(prj_pt_src_t in_shortw,
2005 				    ec_edwards_crv_src_t edwards_crv,
2006 				    aff_pt_edwards_t out_edwards,
2007 				    fp_src_t alpha_edwards)
2008 {
2009 	int ret, iszero;
2010 	aff_pt in_shortw_aff;
2011 	in_shortw_aff.magic = WORD(0);
2012 
2013 	/* Check the curves compatibility */
2014 	ret = prj_pt_check_initialized(in_shortw); EG(ret, err);
2015 	ret = curve_edwards_shortw_check(edwards_crv, in_shortw->crv, alpha_edwards); EG(ret, err);
2016 
2017 	/* Initialize output point with curve */
2018 	ret = aff_pt_init(&in_shortw_aff, in_shortw->crv); EG(ret, err);
2019 
2020 	/* Check if we are the point at infinity
2021 	 * This check induces a non consant time exception, but the current function must be called on
2022 	 * public data anyways.
2023 	 */
2024 	ret = prj_pt_iszero(in_shortw, &iszero); EG(ret, err);
2025 	if(iszero){
2026 		fp zero, one;
2027 		zero.magic = one.magic = WORD(0);
2028 
2029 		ret = fp_init(&zero, in_shortw->X.ctx); EG(ret, err1);
2030 		ret = fp_init(&one, in_shortw->X.ctx); EG(ret, err1);
2031 
2032 		ret = fp_zero(&zero); EG(ret, err1);
2033 		ret = fp_one(&one); EG(ret, err1);
2034 
2035 		ret = aff_pt_edwards_init_from_coords(out_edwards, edwards_crv, &zero, &one);
2036 
2037 err1:
2038 		fp_uninit(&zero);
2039 		fp_uninit(&one);
2040 
2041 		goto err;
2042 	}
2043 
2044 	/* Map projective to affine on the short Weierstrass */
2045 	ret = prj_pt_to_aff(&in_shortw_aff, in_shortw); EG(ret, err);
2046 	/* Use the affine mapping */
2047 	ret = aff_pt_shortw_to_edwards(&in_shortw_aff, edwards_crv, out_edwards, alpha_edwards);
2048 
2049 err:
2050 	aff_pt_uninit(&in_shortw_aff);
2051 
2052 	return ret;
2053 }
2054 
2055 /*
2056  * Map points from Montgomery to short Weierstrass projective points.
2057  *
2058  * The function returns 0 on success, -1 on error.
2059  */
aff_pt_montgomery_to_prj_pt_shortw(aff_pt_montgomery_src_t in_montgomery,ec_shortw_crv_src_t shortw_crv,prj_pt_t out_shortw)2060 int aff_pt_montgomery_to_prj_pt_shortw(aff_pt_montgomery_src_t in_montgomery,
2061 				       ec_shortw_crv_src_t shortw_crv,
2062 				       prj_pt_t out_shortw)
2063 {
2064 	int ret;
2065 	aff_pt out_shortw_aff;
2066 	out_shortw_aff.magic = WORD(0);
2067 
2068 	/* Check the curves compatibility */
2069 	ret = aff_pt_montgomery_check_initialized(in_montgomery); EG(ret, err);
2070 	ret = curve_montgomery_shortw_check(in_montgomery->crv, shortw_crv); EG(ret, err);
2071 
2072 	/* Initialize output point with curve */
2073 	ret = prj_pt_init(out_shortw, shortw_crv); EG(ret, err);
2074 
2075 	/* Use the affine mapping */
2076 	ret = aff_pt_montgomery_to_shortw(in_montgomery, shortw_crv, &out_shortw_aff); EG(ret, err);
2077 	/* And then map the short Weierstrass affine to projective coordinates */
2078 	ret = ec_shortw_aff_to_prj(out_shortw, &out_shortw_aff);
2079 
2080 err:
2081 	aff_pt_uninit(&out_shortw_aff);
2082 
2083 	return ret;
2084 }
2085 
2086 /*
2087  * Map points from short Weierstrass projective points to Montgomery.
2088  *
2089  * The function returns 0 on success, -1 on error.
2090  */
prj_pt_shortw_to_aff_pt_montgomery(prj_pt_src_t in_shortw,ec_montgomery_crv_src_t montgomery_crv,aff_pt_montgomery_t out_montgomery)2091 int prj_pt_shortw_to_aff_pt_montgomery(prj_pt_src_t in_shortw, ec_montgomery_crv_src_t montgomery_crv, aff_pt_montgomery_t out_montgomery)
2092 {
2093 	int ret;
2094 	aff_pt in_shortw_aff;
2095 	in_shortw_aff.magic = WORD(0);
2096 
2097 	/* Check the curves compatibility */
2098 	ret = prj_pt_check_initialized(in_shortw); EG(ret, err);
2099 	ret = curve_montgomery_shortw_check(montgomery_crv, in_shortw->crv); EG(ret, err);
2100 
2101 	/* Initialize output point with curve */
2102 	ret = aff_pt_init(&in_shortw_aff, in_shortw->crv); EG(ret, err);
2103 
2104 	/* Map projective to affine on the short Weierstrass */
2105 	ret = prj_pt_to_aff(&in_shortw_aff, in_shortw); EG(ret, err);
2106 	/* Use the affine mapping */
2107 	ret = aff_pt_shortw_to_montgomery(&in_shortw_aff, montgomery_crv, out_montgomery);
2108 
2109 err:
2110 	aff_pt_uninit(&in_shortw_aff);
2111 
2112 	return ret;
2113 }
2114