xref: /freebsd/contrib/bearssl/src/ec/ec_p256_m62.c (revision cfd6422a5217410fbd66f7a7a8a64d9d85e61229)
1 /*
2  * Copyright (c) 2018 Thomas Pornin <pornin@bolet.org>
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining
5  * a copy of this software and associated documentation files (the
6  * "Software"), to deal in the Software without restriction, including
7  * without limitation the rights to use, copy, modify, merge, publish,
8  * distribute, sublicense, and/or sell copies of the Software, and to
9  * permit persons to whom the Software is furnished to do so, subject to
10  * the following conditions:
11  *
12  * The above copyright notice and this permission notice shall be
13  * included in all copies or substantial portions of the Software.
14  *
15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
18  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
19  * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
20  * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
21  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 
25 #include "inner.h"
26 
27 #if BR_INT128 || BR_UMUL128
28 
29 #if BR_UMUL128
30 #include <intrin.h>
31 #endif
32 
33 static const unsigned char P256_G[] = {
34 	0x04, 0x6B, 0x17, 0xD1, 0xF2, 0xE1, 0x2C, 0x42, 0x47, 0xF8,
35 	0xBC, 0xE6, 0xE5, 0x63, 0xA4, 0x40, 0xF2, 0x77, 0x03, 0x7D,
36 	0x81, 0x2D, 0xEB, 0x33, 0xA0, 0xF4, 0xA1, 0x39, 0x45, 0xD8,
37 	0x98, 0xC2, 0x96, 0x4F, 0xE3, 0x42, 0xE2, 0xFE, 0x1A, 0x7F,
38 	0x9B, 0x8E, 0xE7, 0xEB, 0x4A, 0x7C, 0x0F, 0x9E, 0x16, 0x2B,
39 	0xCE, 0x33, 0x57, 0x6B, 0x31, 0x5E, 0xCE, 0xCB, 0xB6, 0x40,
40 	0x68, 0x37, 0xBF, 0x51, 0xF5
41 };
42 
43 static const unsigned char P256_N[] = {
44 	0xFF, 0xFF, 0xFF, 0xFF, 0x00, 0x00, 0x00, 0x00, 0xFF, 0xFF,
45 	0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xBC, 0xE6, 0xFA, 0xAD,
46 	0xA7, 0x17, 0x9E, 0x84, 0xF3, 0xB9, 0xCA, 0xC2, 0xFC, 0x63,
47 	0x25, 0x51
48 };
49 
50 static const unsigned char *
51 api_generator(int curve, size_t *len)
52 {
53 	(void)curve;
54 	*len = sizeof P256_G;
55 	return P256_G;
56 }
57 
58 static const unsigned char *
59 api_order(int curve, size_t *len)
60 {
61 	(void)curve;
62 	*len = sizeof P256_N;
63 	return P256_N;
64 }
65 
66 static size_t
67 api_xoff(int curve, size_t *len)
68 {
69 	(void)curve;
70 	*len = 32;
71 	return 1;
72 }
73 
74 /*
75  * A field element is encoded as five 64-bit integers, in basis 2^52.
76  * Limbs may occasionally exceed 2^52.
77  *
78  * A _partially reduced_ value is such that the following hold:
79  *   - top limb is less than 2^48 + 2^30
80  *   - the other limbs fit on 53 bits each
81  * In particular, such a value is less than twice the modulus p.
82  */
83 
84 #define BIT(n)   ((uint64_t)1 << (n))
85 #define MASK48   (BIT(48) - BIT(0))
86 #define MASK52   (BIT(52) - BIT(0))
87 
88 /* R = 2^260 mod p */
89 static const uint64_t F256_R[] = {
90 	0x0000000000010, 0xF000000000000, 0xFFFFFFFFFFFFF,
91 	0xFFEFFFFFFFFFF, 0x00000000FFFFF
92 };
93 
94 /* Curve equation is y^2 = x^3 - 3*x + B. This constant is B*R mod p
95    (Montgomery representation of B). */
96 static const uint64_t P256_B_MONTY[] = {
97 	0xDF6229C4BDDFD, 0xCA8843090D89C, 0x212ED6ACF005C,
98 	0x83415A220ABF7, 0x0C30061DD4874
99 };
100 
101 /*
102  * Addition in the field. Carry propagation is not performed.
103  * On input, limbs may be up to 63 bits each; on output, they will
104  * be up to one bit more than on input.
105  */
106 static inline void
107 f256_add(uint64_t *d, const uint64_t *a, const uint64_t *b)
108 {
109 	d[0] = a[0] + b[0];
110 	d[1] = a[1] + b[1];
111 	d[2] = a[2] + b[2];
112 	d[3] = a[3] + b[3];
113 	d[4] = a[4] + b[4];
114 }
115 
116 /*
117  * Partially reduce the provided value.
118  * Input: limbs can go up to 61 bits each.
119  * Output: partially reduced.
120  */
121 static inline void
122 f256_partial_reduce(uint64_t *a)
123 {
124 	uint64_t w, cc, s;
125 
126 	/*
127 	 * Propagate carries.
128 	 */
129 	w = a[0];
130 	a[0] = w & MASK52;
131 	cc = w >> 52;
132 	w = a[1] + cc;
133 	a[1] = w & MASK52;
134 	cc = w >> 52;
135 	w = a[2] + cc;
136 	a[2] = w & MASK52;
137 	cc = w >> 52;
138 	w = a[3] + cc;
139 	a[3] = w & MASK52;
140 	cc = w >> 52;
141 	a[4] += cc;
142 
143 	s = a[4] >> 48;             /* s < 2^14 */
144 	a[0] += s;                  /* a[0] < 2^52 + 2^14 */
145 	w = a[1] - (s << 44);
146 	a[1] = w & MASK52;          /* a[1] < 2^52 */
147 	cc = -(w >> 52) & 0xFFF;    /* cc < 16 */
148 	w = a[2] - cc;
149 	a[2] = w & MASK52;          /* a[2] < 2^52 */
150 	cc = w >> 63;               /* cc = 0 or 1 */
151 	w = a[3] - cc - (s << 36);
152 	a[3] = w & MASK52;          /* a[3] < 2^52 */
153 	cc = w >> 63;               /* cc = 0 or 1 */
154 	w = a[4] & MASK48;
155 	a[4] = w + (s << 16) - cc;  /* a[4] < 2^48 + 2^30 */
156 }
157 
158 /*
159  * Subtraction in the field.
160  * Input: limbs must fit on 60 bits each; in particular, the complete
161  * integer will be less than 2^268 + 2^217.
162  * Output: partially reduced.
163  */
164 static inline void
165 f256_sub(uint64_t *d, const uint64_t *a, const uint64_t *b)
166 {
167 	uint64_t t[5], w, s, cc;
168 
169 	/*
170 	 * We compute d = 2^13*p + a - b; this ensures a positive
171 	 * intermediate value.
172 	 *
173 	 * Each individual addition/subtraction may yield a positive or
174 	 * negative result; thus, we need to handle a signed carry, thus
175 	 * with sign extension. We prefer not to use signed types (int64_t)
176 	 * because conversion from unsigned to signed is cumbersome (a
177 	 * direct cast with the top bit set is undefined behavior; instead,
178 	 * we have to use pointer aliasing, using the guaranteed properties
179 	 * of exact-width types, but this requires the compiler to optimize
180 	 * away the writes and reads from RAM), and right-shifting a
181 	 * signed negative value is implementation-defined. Therefore,
182 	 * we use a custom sign extension.
183 	 */
184 
185 	w = a[0] - b[0] - BIT(13);
186 	t[0] = w & MASK52;
187 	cc = w >> 52;
188 	cc |= -(cc & BIT(11));
189 	w = a[1] - b[1] + cc;
190 	t[1] = w & MASK52;
191 	cc = w >> 52;
192 	cc |= -(cc & BIT(11));
193 	w = a[2] - b[2] + cc;
194 	t[2] = (w & MASK52) + BIT(5);
195 	cc = w >> 52;
196 	cc |= -(cc & BIT(11));
197 	w = a[3] - b[3] + cc;
198 	t[3] = (w & MASK52) + BIT(49);
199 	cc = w >> 52;
200 	cc |= -(cc & BIT(11));
201 	t[4] = (BIT(61) - BIT(29)) + a[4] - b[4] + cc;
202 
203 	/*
204 	 * Perform partial reduction. Rule is:
205 	 *  2^256 = 2^224 - 2^192 - 2^96 + 1 mod p
206 	 *
207 	 * At that point:
208 	 *    0 <= t[0] <= 2^52 - 1
209 	 *    0 <= t[1] <= 2^52 - 1
210 	 *    2^5 <= t[2] <= 2^52 + 2^5 - 1
211 	 *    2^49 <= t[3] <= 2^52 + 2^49 - 1
212 	 *    2^59 < t[4] <= 2^61 + 2^60 - 2^29
213 	 *
214 	 * Thus, the value 's' (t[4] / 2^48) will be necessarily
215 	 * greater than 2048, and less than 12288.
216 	 */
217 	s = t[4] >> 48;
218 
219 	d[0] = t[0] + s;             /* d[0] <= 2^52 + 12287 */
220 	w = t[1] - (s << 44);
221 	d[1] = w & MASK52;           /* d[1] <= 2^52 - 1 */
222 	cc = -(w >> 52) & 0xFFF;     /* cc <= 48 */
223 	w = t[2] - cc;
224 	cc = w >> 63;                /* cc = 0 or 1 */
225 	d[2] = w + (cc << 52);       /* d[2] <= 2^52 + 31 */
226 	w = t[3] - cc - (s << 36);
227 	cc = w >> 63;                /* cc = 0 or 1 */
228 	d[3] = w + (cc << 52);       /* t[3] <= 2^52 + 2^49 - 1 */
229 	d[4] = (t[4] & MASK48) + (s << 16) - cc;  /* d[4] < 2^48 + 2^30 */
230 
231 	/*
232 	 * If s = 0, then none of the limbs is modified, and there cannot
233 	 * be an overflow; if s != 0, then (s << 16) > cc, and there is
234 	 * no overflow either.
235 	 */
236 }
237 
238 /*
239  * Montgomery multiplication in the field.
240  * Input: limbs must fit on 56 bits each.
241  * Output: partially reduced.
242  */
243 static void
244 f256_montymul(uint64_t *d, const uint64_t *a, const uint64_t *b)
245 {
246 #if BR_INT128
247 
248 	int i;
249 	uint64_t t[5];
250 
251 	t[0] = 0;
252 	t[1] = 0;
253 	t[2] = 0;
254 	t[3] = 0;
255 	t[4] = 0;
256 	for (i = 0; i < 5; i ++) {
257 		uint64_t x, f, cc, w, s;
258 		unsigned __int128 z;
259 
260 		/*
261 		 * Since limbs of a[] and b[] fit on 56 bits each,
262 		 * each individual product fits on 112 bits. Also,
263 		 * the factor f fits on 52 bits, so f<<48 fits on
264 		 * 112 bits too. This guarantees that carries (cc)
265 		 * will fit on 62 bits, thus no overflow.
266 		 *
267 		 * The operations below compute:
268 		 *   t <- (t + x*b + f*p) / 2^64
269 		 */
270 		x = a[i];
271 		z = (unsigned __int128)b[0] * (unsigned __int128)x
272 			+ (unsigned __int128)t[0];
273 		f = (uint64_t)z & MASK52;
274 		cc = (uint64_t)(z >> 52);
275 		z = (unsigned __int128)b[1] * (unsigned __int128)x
276 			+ (unsigned __int128)t[1] + cc
277 			+ ((unsigned __int128)f << 44);
278 		t[0] = (uint64_t)z & MASK52;
279 		cc = (uint64_t)(z >> 52);
280 		z = (unsigned __int128)b[2] * (unsigned __int128)x
281 			+ (unsigned __int128)t[2] + cc;
282 		t[1] = (uint64_t)z & MASK52;
283 		cc = (uint64_t)(z >> 52);
284 		z = (unsigned __int128)b[3] * (unsigned __int128)x
285 			+ (unsigned __int128)t[3] + cc
286 			+ ((unsigned __int128)f << 36);
287 		t[2] = (uint64_t)z & MASK52;
288 		cc = (uint64_t)(z >> 52);
289 		z = (unsigned __int128)b[4] * (unsigned __int128)x
290 			+ (unsigned __int128)t[4] + cc
291 			+ ((unsigned __int128)f << 48)
292 			- ((unsigned __int128)f << 16);
293 		t[3] = (uint64_t)z & MASK52;
294 		t[4] = (uint64_t)(z >> 52);
295 
296 		/*
297 		 * t[4] may be up to 62 bits here; we need to do a
298 		 * partial reduction. Note that limbs t[0] to t[3]
299 		 * fit on 52 bits each.
300 		 */
301 		s = t[4] >> 48;             /* s < 2^14 */
302 		t[0] += s;                  /* t[0] < 2^52 + 2^14 */
303 		w = t[1] - (s << 44);
304 		t[1] = w & MASK52;          /* t[1] < 2^52 */
305 		cc = -(w >> 52) & 0xFFF;    /* cc < 16 */
306 		w = t[2] - cc;
307 		t[2] = w & MASK52;          /* t[2] < 2^52 */
308 		cc = w >> 63;               /* cc = 0 or 1 */
309 		w = t[3] - cc - (s << 36);
310 		t[3] = w & MASK52;          /* t[3] < 2^52 */
311 		cc = w >> 63;               /* cc = 0 or 1 */
312 		w = t[4] & MASK48;
313 		t[4] = w + (s << 16) - cc;  /* t[4] < 2^48 + 2^30 */
314 
315 		/*
316 		 * The final t[4] cannot overflow because cc is 0 or 1,
317 		 * and cc can be 1 only if s != 0.
318 		 */
319 	}
320 
321 	d[0] = t[0];
322 	d[1] = t[1];
323 	d[2] = t[2];
324 	d[3] = t[3];
325 	d[4] = t[4];
326 
327 #elif BR_UMUL128
328 
329 	int i;
330 	uint64_t t[5];
331 
332 	t[0] = 0;
333 	t[1] = 0;
334 	t[2] = 0;
335 	t[3] = 0;
336 	t[4] = 0;
337 	for (i = 0; i < 5; i ++) {
338 		uint64_t x, f, cc, w, s, zh, zl;
339 		unsigned char k;
340 
341 		/*
342 		 * Since limbs of a[] and b[] fit on 56 bits each,
343 		 * each individual product fits on 112 bits. Also,
344 		 * the factor f fits on 52 bits, so f<<48 fits on
345 		 * 112 bits too. This guarantees that carries (cc)
346 		 * will fit on 62 bits, thus no overflow.
347 		 *
348 		 * The operations below compute:
349 		 *   t <- (t + x*b + f*p) / 2^64
350 		 */
351 		x = a[i];
352 		zl = _umul128(b[0], x, &zh);
353 		k = _addcarry_u64(0, t[0], zl, &zl);
354 		(void)_addcarry_u64(k, 0, zh, &zh);
355 		f = zl & MASK52;
356 		cc = (zl >> 52) | (zh << 12);
357 
358 		zl = _umul128(b[1], x, &zh);
359 		k = _addcarry_u64(0, t[1], zl, &zl);
360 		(void)_addcarry_u64(k, 0, zh, &zh);
361 		k = _addcarry_u64(0, cc, zl, &zl);
362 		(void)_addcarry_u64(k, 0, zh, &zh);
363 		k = _addcarry_u64(0, f << 44, zl, &zl);
364 		(void)_addcarry_u64(k, f >> 20, zh, &zh);
365 		t[0] = zl & MASK52;
366 		cc = (zl >> 52) | (zh << 12);
367 
368 		zl = _umul128(b[2], x, &zh);
369 		k = _addcarry_u64(0, t[2], zl, &zl);
370 		(void)_addcarry_u64(k, 0, zh, &zh);
371 		k = _addcarry_u64(0, cc, zl, &zl);
372 		(void)_addcarry_u64(k, 0, zh, &zh);
373 		t[1] = zl & MASK52;
374 		cc = (zl >> 52) | (zh << 12);
375 
376 		zl = _umul128(b[3], x, &zh);
377 		k = _addcarry_u64(0, t[3], zl, &zl);
378 		(void)_addcarry_u64(k, 0, zh, &zh);
379 		k = _addcarry_u64(0, cc, zl, &zl);
380 		(void)_addcarry_u64(k, 0, zh, &zh);
381 		k = _addcarry_u64(0, f << 36, zl, &zl);
382 		(void)_addcarry_u64(k, f >> 28, zh, &zh);
383 		t[2] = zl & MASK52;
384 		cc = (zl >> 52) | (zh << 12);
385 
386 		zl = _umul128(b[4], x, &zh);
387 		k = _addcarry_u64(0, t[4], zl, &zl);
388 		(void)_addcarry_u64(k, 0, zh, &zh);
389 		k = _addcarry_u64(0, cc, zl, &zl);
390 		(void)_addcarry_u64(k, 0, zh, &zh);
391 		k = _addcarry_u64(0, f << 48, zl, &zl);
392 		(void)_addcarry_u64(k, f >> 16, zh, &zh);
393 		k = _subborrow_u64(0, zl, f << 16, &zl);
394 		(void)_subborrow_u64(k, zh, f >> 48, &zh);
395 		t[3] = zl & MASK52;
396 		t[4] = (zl >> 52) | (zh << 12);
397 
398 		/*
399 		 * t[4] may be up to 62 bits here; we need to do a
400 		 * partial reduction. Note that limbs t[0] to t[3]
401 		 * fit on 52 bits each.
402 		 */
403 		s = t[4] >> 48;             /* s < 2^14 */
404 		t[0] += s;                  /* t[0] < 2^52 + 2^14 */
405 		w = t[1] - (s << 44);
406 		t[1] = w & MASK52;          /* t[1] < 2^52 */
407 		cc = -(w >> 52) & 0xFFF;    /* cc < 16 */
408 		w = t[2] - cc;
409 		t[2] = w & MASK52;          /* t[2] < 2^52 */
410 		cc = w >> 63;               /* cc = 0 or 1 */
411 		w = t[3] - cc - (s << 36);
412 		t[3] = w & MASK52;          /* t[3] < 2^52 */
413 		cc = w >> 63;               /* cc = 0 or 1 */
414 		w = t[4] & MASK48;
415 		t[4] = w + (s << 16) - cc;  /* t[4] < 2^48 + 2^30 */
416 
417 		/*
418 		 * The final t[4] cannot overflow because cc is 0 or 1,
419 		 * and cc can be 1 only if s != 0.
420 		 */
421 	}
422 
423 	d[0] = t[0];
424 	d[1] = t[1];
425 	d[2] = t[2];
426 	d[3] = t[3];
427 	d[4] = t[4];
428 
429 #endif
430 }
431 
432 /*
433  * Montgomery squaring in the field; currently a basic wrapper around
434  * multiplication (inline, should be optimized away).
435  * TODO: see if some extra speed can be gained here.
436  */
437 static inline void
438 f256_montysquare(uint64_t *d, const uint64_t *a)
439 {
440 	f256_montymul(d, a, a);
441 }
442 
443 /*
444  * Convert to Montgomery representation.
445  */
446 static void
447 f256_tomonty(uint64_t *d, const uint64_t *a)
448 {
449 	/*
450 	 * R2 = 2^520 mod p.
451 	 * If R = 2^260 mod p, then R2 = R^2 mod p; and the Montgomery
452 	 * multiplication of a by R2 is: a*R2/R = a*R mod p, i.e. the
453 	 * conversion to Montgomery representation.
454 	 */
455 	static const uint64_t R2[] = {
456 		0x0000000000300, 0xFFFFFFFF00000, 0xFFFFEFFFFFFFB,
457 		0xFDFFFFFFFFFFF, 0x0000004FFFFFF
458 	};
459 
460 	f256_montymul(d, a, R2);
461 }
462 
463 /*
464  * Convert from Montgomery representation.
465  */
466 static void
467 f256_frommonty(uint64_t *d, const uint64_t *a)
468 {
469 	/*
470 	 * Montgomery multiplication by 1 is division by 2^260 modulo p.
471 	 */
472 	static const uint64_t one[] = { 1, 0, 0, 0, 0 };
473 
474 	f256_montymul(d, a, one);
475 }
476 
477 /*
478  * Inversion in the field. If the source value is 0 modulo p, then this
479  * returns 0 or p. This function uses Montgomery representation.
480  */
481 static void
482 f256_invert(uint64_t *d, const uint64_t *a)
483 {
484 	/*
485 	 * We compute a^(p-2) mod p. The exponent pattern (from high to
486 	 * low) is:
487 	 *  - 32 bits of value 1
488 	 *  - 31 bits of value 0
489 	 *  - 1 bit of value 1
490 	 *  - 96 bits of value 0
491 	 *  - 94 bits of value 1
492 	 *  - 1 bit of value 0
493 	 *  - 1 bit of value 1
494 	 * To speed up the square-and-multiply algorithm, we precompute
495 	 * a^(2^31-1).
496 	 */
497 
498 	uint64_t r[5], t[5];
499 	int i;
500 
501 	memcpy(t, a, sizeof t);
502 	for (i = 0; i < 30; i ++) {
503 		f256_montysquare(t, t);
504 		f256_montymul(t, t, a);
505 	}
506 
507 	memcpy(r, t, sizeof t);
508 	for (i = 224; i >= 0; i --) {
509 		f256_montysquare(r, r);
510 		switch (i) {
511 		case 0:
512 		case 2:
513 		case 192:
514 		case 224:
515 			f256_montymul(r, r, a);
516 			break;
517 		case 3:
518 		case 34:
519 		case 65:
520 			f256_montymul(r, r, t);
521 			break;
522 		}
523 	}
524 	memcpy(d, r, sizeof r);
525 }
526 
527 /*
528  * Finalize reduction.
529  * Input value should be partially reduced.
530  * On output, limbs a[0] to a[3] fit on 52 bits each, limb a[4] fits
531  * on 48 bits, and the integer is less than p.
532  */
533 static inline void
534 f256_final_reduce(uint64_t *a)
535 {
536 	uint64_t r[5], t[5], w, cc;
537 	int i;
538 
539 	/*
540 	 * Propagate carries to ensure that limbs 0 to 3 fit on 52 bits.
541 	 */
542 	cc = 0;
543 	for (i = 0; i < 5; i ++) {
544 		w = a[i] + cc;
545 		r[i] = w & MASK52;
546 		cc = w >> 52;
547 	}
548 
549 	/*
550 	 * We compute t = r + (2^256 - p) = r + 2^224 - 2^192 - 2^96 + 1.
551 	 * If t < 2^256, then r < p, and we return r. Otherwise, we
552 	 * want to return r - p = t - 2^256.
553 	 */
554 
555 	/*
556 	 * Add 2^224 + 1, and propagate carries to ensure that limbs
557 	 * t[0] to t[3] fit in 52 bits each.
558 	 */
559 	w = r[0] + 1;
560 	t[0] = w & MASK52;
561 	cc = w >> 52;
562 	w = r[1] + cc;
563 	t[1] = w & MASK52;
564 	cc = w >> 52;
565 	w = r[2] + cc;
566 	t[2] = w & MASK52;
567 	cc = w >> 52;
568 	w = r[3] + cc;
569 	t[3] = w & MASK52;
570 	cc = w >> 52;
571 	t[4] = r[4] + cc + BIT(16);
572 
573 	/*
574 	 * Subtract 2^192 + 2^96. Since we just added 2^224 + 1, the
575 	 * result cannot be negative.
576 	 */
577 	w = t[1] - BIT(44);
578 	t[1] = w & MASK52;
579 	cc = w >> 63;
580 	w = t[2] - cc;
581 	t[2] = w & MASK52;
582 	cc = w >> 63;
583 	w = t[3] - BIT(36);
584 	t[3] = w & MASK52;
585 	cc = w >> 63;
586 	t[4] -= cc;
587 
588 	/*
589 	 * If the top limb t[4] fits on 48 bits, then r[] is already
590 	 * in the proper range. Otherwise, t[] is the value to return
591 	 * (truncated to 256 bits).
592 	 */
593 	cc = -(t[4] >> 48);
594 	t[4] &= MASK48;
595 	for (i = 0; i < 5; i ++) {
596 		a[i] = r[i] ^ (cc & (r[i] ^ t[i]));
597 	}
598 }
599 
600 /*
601  * Points in affine and Jacobian coordinates.
602  *
603  *  - In affine coordinates, the point-at-infinity cannot be encoded.
604  *  - Jacobian coordinates (X,Y,Z) correspond to affine (X/Z^2,Y/Z^3);
605  *    if Z = 0 then this is the point-at-infinity.
606  */
607 typedef struct {
608 	uint64_t x[5];
609 	uint64_t y[5];
610 } p256_affine;
611 
612 typedef struct {
613 	uint64_t x[5];
614 	uint64_t y[5];
615 	uint64_t z[5];
616 } p256_jacobian;
617 
618 /*
619  * Decode a field element (unsigned big endian notation).
620  */
621 static void
622 f256_decode(uint64_t *a, const unsigned char *buf)
623 {
624 	uint64_t w0, w1, w2, w3;
625 
626 	w3 = br_dec64be(buf +  0);
627 	w2 = br_dec64be(buf +  8);
628 	w1 = br_dec64be(buf + 16);
629 	w0 = br_dec64be(buf + 24);
630 	a[0] = w0 & MASK52;
631 	a[1] = ((w0 >> 52) | (w1 << 12)) & MASK52;
632 	a[2] = ((w1 >> 40) | (w2 << 24)) & MASK52;
633 	a[3] = ((w2 >> 28) | (w3 << 36)) & MASK52;
634 	a[4] = w3 >> 16;
635 }
636 
637 /*
638  * Encode a field element (unsigned big endian notation). The field
639  * element MUST be fully reduced.
640  */
641 static void
642 f256_encode(unsigned char *buf, const uint64_t *a)
643 {
644 	uint64_t w0, w1, w2, w3;
645 
646 	w0 = a[0] | (a[1] << 52);
647 	w1 = (a[1] >> 12) | (a[2] << 40);
648 	w2 = (a[2] >> 24) | (a[3] << 28);
649 	w3 = (a[3] >> 36) | (a[4] << 16);
650 	br_enc64be(buf +  0, w3);
651 	br_enc64be(buf +  8, w2);
652 	br_enc64be(buf + 16, w1);
653 	br_enc64be(buf + 24, w0);
654 }
655 
656 /*
657  * Decode a point. The returned point is in Jacobian coordinates, but
658  * with z = 1. If the encoding is invalid, or encodes a point which is
659  * not on the curve, or encodes the point at infinity, then this function
660  * returns 0. Otherwise, 1 is returned.
661  *
662  * The buffer is assumed to have length exactly 65 bytes.
663  */
664 static uint32_t
665 point_decode(p256_jacobian *P, const unsigned char *buf)
666 {
667 	uint64_t x[5], y[5], t[5], x3[5], tt;
668 	uint32_t r;
669 
670 	/*
671 	 * Header byte shall be 0x04.
672 	 */
673 	r = EQ(buf[0], 0x04);
674 
675 	/*
676 	 * Decode X and Y coordinates, and convert them into
677 	 * Montgomery representation.
678 	 */
679 	f256_decode(x, buf +  1);
680 	f256_decode(y, buf + 33);
681 	f256_tomonty(x, x);
682 	f256_tomonty(y, y);
683 
684 	/*
685 	 * Verify y^2 = x^3 + A*x + B. In curve P-256, A = -3.
686 	 * Note that the Montgomery representation of 0 is 0. We must
687 	 * take care to apply the final reduction to make sure we have
688 	 * 0 and not p.
689 	 */
690 	f256_montysquare(t, y);
691 	f256_montysquare(x3, x);
692 	f256_montymul(x3, x3, x);
693 	f256_sub(t, t, x3);
694 	f256_add(t, t, x);
695 	f256_add(t, t, x);
696 	f256_add(t, t, x);
697 	f256_sub(t, t, P256_B_MONTY);
698 	f256_final_reduce(t);
699 	tt = t[0] | t[1] | t[2] | t[3] | t[4];
700 	r &= EQ((uint32_t)(tt | (tt >> 32)), 0);
701 
702 	/*
703 	 * Return the point in Jacobian coordinates (and Montgomery
704 	 * representation).
705 	 */
706 	memcpy(P->x, x, sizeof x);
707 	memcpy(P->y, y, sizeof y);
708 	memcpy(P->z, F256_R, sizeof F256_R);
709 	return r;
710 }
711 
712 /*
713  * Final conversion for a point:
714  *  - The point is converted back to affine coordinates.
715  *  - Final reduction is performed.
716  *  - The point is encoded into the provided buffer.
717  *
718  * If the point is the point-at-infinity, all operations are performed,
719  * but the buffer contents are indeterminate, and 0 is returned. Otherwise,
720  * the encoded point is written in the buffer, and 1 is returned.
721  */
722 static uint32_t
723 point_encode(unsigned char *buf, const p256_jacobian *P)
724 {
725 	uint64_t t1[5], t2[5], z;
726 
727 	/* Set t1 = 1/z^2 and t2 = 1/z^3. */
728 	f256_invert(t2, P->z);
729 	f256_montysquare(t1, t2);
730 	f256_montymul(t2, t2, t1);
731 
732 	/* Compute affine coordinates x (in t1) and y (in t2). */
733 	f256_montymul(t1, P->x, t1);
734 	f256_montymul(t2, P->y, t2);
735 
736 	/* Convert back from Montgomery representation, and finalize
737 	   reductions. */
738 	f256_frommonty(t1, t1);
739 	f256_frommonty(t2, t2);
740 	f256_final_reduce(t1);
741 	f256_final_reduce(t2);
742 
743 	/* Encode. */
744 	buf[0] = 0x04;
745 	f256_encode(buf +  1, t1);
746 	f256_encode(buf + 33, t2);
747 
748 	/* Return success if and only if P->z != 0. */
749 	z = P->z[0] | P->z[1] | P->z[2] | P->z[3] | P->z[4];
750 	return NEQ((uint32_t)(z | z >> 32), 0);
751 }
752 
753 /*
754  * Point doubling in Jacobian coordinates: point P is doubled.
755  * Note: if the source point is the point-at-infinity, then the result is
756  * still the point-at-infinity, which is correct. Moreover, if the three
757  * coordinates were zero, then they still are zero in the returned value.
758  */
759 static void
760 p256_double(p256_jacobian *P)
761 {
762 	/*
763 	 * Doubling formulas are:
764 	 *
765 	 *   s = 4*x*y^2
766 	 *   m = 3*(x + z^2)*(x - z^2)
767 	 *   x' = m^2 - 2*s
768 	 *   y' = m*(s - x') - 8*y^4
769 	 *   z' = 2*y*z
770 	 *
771 	 * These formulas work for all points, including points of order 2
772 	 * and points at infinity:
773 	 *   - If y = 0 then z' = 0. But there is no such point in P-256
774 	 *     anyway.
775 	 *   - If z = 0 then z' = 0.
776 	 */
777 	uint64_t t1[5], t2[5], t3[5], t4[5];
778 
779 	/*
780 	 * Compute z^2 in t1.
781 	 */
782 	f256_montysquare(t1, P->z);
783 
784 	/*
785 	 * Compute x-z^2 in t2 and x+z^2 in t1.
786 	 */
787 	f256_add(t2, P->x, t1);
788 	f256_sub(t1, P->x, t1);
789 
790 	/*
791 	 * Compute 3*(x+z^2)*(x-z^2) in t1.
792 	 */
793 	f256_montymul(t3, t1, t2);
794 	f256_add(t1, t3, t3);
795 	f256_add(t1, t3, t1);
796 
797 	/*
798 	 * Compute 4*x*y^2 (in t2) and 2*y^2 (in t3).
799 	 */
800 	f256_montysquare(t3, P->y);
801 	f256_add(t3, t3, t3);
802 	f256_montymul(t2, P->x, t3);
803 	f256_add(t2, t2, t2);
804 
805 	/*
806 	 * Compute x' = m^2 - 2*s.
807 	 */
808 	f256_montysquare(P->x, t1);
809 	f256_sub(P->x, P->x, t2);
810 	f256_sub(P->x, P->x, t2);
811 
812 	/*
813 	 * Compute z' = 2*y*z.
814 	 */
815 	f256_montymul(t4, P->y, P->z);
816 	f256_add(P->z, t4, t4);
817 	f256_partial_reduce(P->z);
818 
819 	/*
820 	 * Compute y' = m*(s - x') - 8*y^4. Note that we already have
821 	 * 2*y^2 in t3.
822 	 */
823 	f256_sub(t2, t2, P->x);
824 	f256_montymul(P->y, t1, t2);
825 	f256_montysquare(t4, t3);
826 	f256_add(t4, t4, t4);
827 	f256_sub(P->y, P->y, t4);
828 }
829 
830 /*
831  * Point addition (Jacobian coordinates): P1 is replaced with P1+P2.
832  * This function computes the wrong result in the following cases:
833  *
834  *   - If P1 == 0 but P2 != 0
835  *   - If P1 != 0 but P2 == 0
836  *   - If P1 == P2
837  *
838  * In all three cases, P1 is set to the point at infinity.
839  *
840  * Returned value is 0 if one of the following occurs:
841  *
842  *   - P1 and P2 have the same Y coordinate.
843  *   - P1 == 0 and P2 == 0.
844  *   - The Y coordinate of one of the points is 0 and the other point is
845  *     the point at infinity.
846  *
847  * The third case cannot actually happen with valid points, since a point
848  * with Y == 0 is a point of order 2, and there is no point of order 2 on
849  * curve P-256.
850  *
851  * Therefore, assuming that P1 != 0 and P2 != 0 on input, then the caller
852  * can apply the following:
853  *
854  *   - If the result is not the point at infinity, then it is correct.
855  *   - Otherwise, if the returned value is 1, then this is a case of
856  *     P1+P2 == 0, so the result is indeed the point at infinity.
857  *   - Otherwise, P1 == P2, so a "double" operation should have been
858  *     performed.
859  *
860  * Note that you can get a returned value of 0 with a correct result,
861  * e.g. if P1 and P2 have the same Y coordinate, but distinct X coordinates.
862  */
863 static uint32_t
864 p256_add(p256_jacobian *P1, const p256_jacobian *P2)
865 {
866 	/*
867 	 * Addtions formulas are:
868 	 *
869 	 *   u1 = x1 * z2^2
870 	 *   u2 = x2 * z1^2
871 	 *   s1 = y1 * z2^3
872 	 *   s2 = y2 * z1^3
873 	 *   h = u2 - u1
874 	 *   r = s2 - s1
875 	 *   x3 = r^2 - h^3 - 2 * u1 * h^2
876 	 *   y3 = r * (u1 * h^2 - x3) - s1 * h^3
877 	 *   z3 = h * z1 * z2
878 	 */
879 	uint64_t t1[5], t2[5], t3[5], t4[5], t5[5], t6[5], t7[5], tt;
880 	uint32_t ret;
881 
882 	/*
883 	 * Compute u1 = x1*z2^2 (in t1) and s1 = y1*z2^3 (in t3).
884 	 */
885 	f256_montysquare(t3, P2->z);
886 	f256_montymul(t1, P1->x, t3);
887 	f256_montymul(t4, P2->z, t3);
888 	f256_montymul(t3, P1->y, t4);
889 
890 	/*
891 	 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4).
892 	 */
893 	f256_montysquare(t4, P1->z);
894 	f256_montymul(t2, P2->x, t4);
895 	f256_montymul(t5, P1->z, t4);
896 	f256_montymul(t4, P2->y, t5);
897 
898 	/*
899 	 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4).
900 	 * We need to test whether r is zero, so we will do some extra
901 	 * reduce.
902 	 */
903 	f256_sub(t2, t2, t1);
904 	f256_sub(t4, t4, t3);
905 	f256_final_reduce(t4);
906 	tt = t4[0] | t4[1] | t4[2] | t4[3] | t4[4];
907 	ret = (uint32_t)(tt | (tt >> 32));
908 	ret = (ret | -ret) >> 31;
909 
910 	/*
911 	 * Compute u1*h^2 (in t6) and h^3 (in t5);
912 	 */
913 	f256_montysquare(t7, t2);
914 	f256_montymul(t6, t1, t7);
915 	f256_montymul(t5, t7, t2);
916 
917 	/*
918 	 * Compute x3 = r^2 - h^3 - 2*u1*h^2.
919 	 */
920 	f256_montysquare(P1->x, t4);
921 	f256_sub(P1->x, P1->x, t5);
922 	f256_sub(P1->x, P1->x, t6);
923 	f256_sub(P1->x, P1->x, t6);
924 
925 	/*
926 	 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3.
927 	 */
928 	f256_sub(t6, t6, P1->x);
929 	f256_montymul(P1->y, t4, t6);
930 	f256_montymul(t1, t5, t3);
931 	f256_sub(P1->y, P1->y, t1);
932 
933 	/*
934 	 * Compute z3 = h*z1*z2.
935 	 */
936 	f256_montymul(t1, P1->z, P2->z);
937 	f256_montymul(P1->z, t1, t2);
938 
939 	return ret;
940 }
941 
942 /*
943  * Point addition (mixed coordinates): P1 is replaced with P1+P2.
944  * This is a specialised function for the case when P2 is a non-zero point
945  * in affine coordinates.
946  *
947  * This function computes the wrong result in the following cases:
948  *
949  *   - If P1 == 0
950  *   - If P1 == P2
951  *
952  * In both cases, P1 is set to the point at infinity.
953  *
954  * Returned value is 0 if one of the following occurs:
955  *
956  *   - P1 and P2 have the same Y (affine) coordinate.
957  *   - The Y coordinate of P2 is 0 and P1 is the point at infinity.
958  *
959  * The second case cannot actually happen with valid points, since a point
960  * with Y == 0 is a point of order 2, and there is no point of order 2 on
961  * curve P-256.
962  *
963  * Therefore, assuming that P1 != 0 on input, then the caller
964  * can apply the following:
965  *
966  *   - If the result is not the point at infinity, then it is correct.
967  *   - Otherwise, if the returned value is 1, then this is a case of
968  *     P1+P2 == 0, so the result is indeed the point at infinity.
969  *   - Otherwise, P1 == P2, so a "double" operation should have been
970  *     performed.
971  *
972  * Again, a value of 0 may be returned in some cases where the addition
973  * result is correct.
974  */
975 static uint32_t
976 p256_add_mixed(p256_jacobian *P1, const p256_affine *P2)
977 {
978 	/*
979 	 * Addtions formulas are:
980 	 *
981 	 *   u1 = x1
982 	 *   u2 = x2 * z1^2
983 	 *   s1 = y1
984 	 *   s2 = y2 * z1^3
985 	 *   h = u2 - u1
986 	 *   r = s2 - s1
987 	 *   x3 = r^2 - h^3 - 2 * u1 * h^2
988 	 *   y3 = r * (u1 * h^2 - x3) - s1 * h^3
989 	 *   z3 = h * z1
990 	 */
991 	uint64_t t1[5], t2[5], t3[5], t4[5], t5[5], t6[5], t7[5], tt;
992 	uint32_t ret;
993 
994 	/*
995 	 * Compute u1 = x1 (in t1) and s1 = y1 (in t3).
996 	 */
997 	memcpy(t1, P1->x, sizeof t1);
998 	memcpy(t3, P1->y, sizeof t3);
999 
1000 	/*
1001 	 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4).
1002 	 */
1003 	f256_montysquare(t4, P1->z);
1004 	f256_montymul(t2, P2->x, t4);
1005 	f256_montymul(t5, P1->z, t4);
1006 	f256_montymul(t4, P2->y, t5);
1007 
1008 	/*
1009 	 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4).
1010 	 * We need to test whether r is zero, so we will do some extra
1011 	 * reduce.
1012 	 */
1013 	f256_sub(t2, t2, t1);
1014 	f256_sub(t4, t4, t3);
1015 	f256_final_reduce(t4);
1016 	tt = t4[0] | t4[1] | t4[2] | t4[3] | t4[4];
1017 	ret = (uint32_t)(tt | (tt >> 32));
1018 	ret = (ret | -ret) >> 31;
1019 
1020 	/*
1021 	 * Compute u1*h^2 (in t6) and h^3 (in t5);
1022 	 */
1023 	f256_montysquare(t7, t2);
1024 	f256_montymul(t6, t1, t7);
1025 	f256_montymul(t5, t7, t2);
1026 
1027 	/*
1028 	 * Compute x3 = r^2 - h^3 - 2*u1*h^2.
1029 	 */
1030 	f256_montysquare(P1->x, t4);
1031 	f256_sub(P1->x, P1->x, t5);
1032 	f256_sub(P1->x, P1->x, t6);
1033 	f256_sub(P1->x, P1->x, t6);
1034 
1035 	/*
1036 	 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3.
1037 	 */
1038 	f256_sub(t6, t6, P1->x);
1039 	f256_montymul(P1->y, t4, t6);
1040 	f256_montymul(t1, t5, t3);
1041 	f256_sub(P1->y, P1->y, t1);
1042 
1043 	/*
1044 	 * Compute z3 = h*z1*z2.
1045 	 */
1046 	f256_montymul(P1->z, P1->z, t2);
1047 
1048 	return ret;
1049 }
1050 
1051 #if 0
1052 /* unused */
1053 /*
1054  * Point addition (mixed coordinates, complete): P1 is replaced with P1+P2.
1055  * This is a specialised function for the case when P2 is a non-zero point
1056  * in affine coordinates.
1057  *
1058  * This function returns the correct result in all cases.
1059  */
1060 static uint32_t
1061 p256_add_complete_mixed(p256_jacobian *P1, const p256_affine *P2)
1062 {
1063 	/*
1064 	 * Addtions formulas, in the general case, are:
1065 	 *
1066 	 *   u1 = x1
1067 	 *   u2 = x2 * z1^2
1068 	 *   s1 = y1
1069 	 *   s2 = y2 * z1^3
1070 	 *   h = u2 - u1
1071 	 *   r = s2 - s1
1072 	 *   x3 = r^2 - h^3 - 2 * u1 * h^2
1073 	 *   y3 = r * (u1 * h^2 - x3) - s1 * h^3
1074 	 *   z3 = h * z1
1075 	 *
1076 	 * These formulas mishandle the two following cases:
1077 	 *
1078 	 *  - If P1 is the point-at-infinity (z1 = 0), then z3 is
1079 	 *    incorrectly set to 0.
1080 	 *
1081 	 *  - If P1 = P2, then u1 = u2 and s1 = s2, and x3, y3 and z3
1082 	 *    are all set to 0.
1083 	 *
1084 	 * However, if P1 + P2 = 0, then u1 = u2 but s1 != s2, and then
1085 	 * we correctly get z3 = 0 (the point-at-infinity).
1086 	 *
1087 	 * To fix the case P1 = 0, we perform at the end a copy of P2
1088 	 * over P1, conditional to z1 = 0.
1089 	 *
1090 	 * For P1 = P2: in that case, both h and r are set to 0, and
1091 	 * we get x3, y3 and z3 equal to 0. We can test for that
1092 	 * occurrence to make a mask which will be all-one if P1 = P2,
1093 	 * or all-zero otherwise; then we can compute the double of P2
1094 	 * and add it, combined with the mask, to (x3,y3,z3).
1095 	 *
1096 	 * Using the doubling formulas in p256_double() on (x2,y2),
1097 	 * simplifying since P2 is affine (i.e. z2 = 1, implicitly),
1098 	 * we get:
1099 	 *   s = 4*x2*y2^2
1100 	 *   m = 3*(x2 + 1)*(x2 - 1)
1101 	 *   x' = m^2 - 2*s
1102 	 *   y' = m*(s - x') - 8*y2^4
1103 	 *   z' = 2*y2
1104 	 * which requires only 6 multiplications. Added to the 11
1105 	 * multiplications of the normal mixed addition in Jacobian
1106 	 * coordinates, we get a cost of 17 multiplications in total.
1107 	 */
1108 	uint64_t t1[5], t2[5], t3[5], t4[5], t5[5], t6[5], t7[5], tt, zz;
1109 	int i;
1110 
1111 	/*
1112 	 * Set zz to -1 if P1 is the point at infinity, 0 otherwise.
1113 	 */
1114 	zz = P1->z[0] | P1->z[1] | P1->z[2] | P1->z[3] | P1->z[4];
1115 	zz = ((zz | -zz) >> 63) - (uint64_t)1;
1116 
1117 	/*
1118 	 * Compute u1 = x1 (in t1) and s1 = y1 (in t3).
1119 	 */
1120 	memcpy(t1, P1->x, sizeof t1);
1121 	memcpy(t3, P1->y, sizeof t3);
1122 
1123 	/*
1124 	 * Compute u2 = x2*z1^2 (in t2) and s2 = y2*z1^3 (in t4).
1125 	 */
1126 	f256_montysquare(t4, P1->z);
1127 	f256_montymul(t2, P2->x, t4);
1128 	f256_montymul(t5, P1->z, t4);
1129 	f256_montymul(t4, P2->y, t5);
1130 
1131 	/*
1132 	 * Compute h = h2 - u1 (in t2) and r = s2 - s1 (in t4).
1133 	 * reduce.
1134 	 */
1135 	f256_sub(t2, t2, t1);
1136 	f256_sub(t4, t4, t3);
1137 
1138 	/*
1139 	 * If both h = 0 and r = 0, then P1 = P2, and we want to set
1140 	 * the mask tt to -1; otherwise, the mask will be 0.
1141 	 */
1142 	f256_final_reduce(t2);
1143 	f256_final_reduce(t4);
1144 	tt = t2[0] | t2[1] | t2[2] | t2[3] | t2[4]
1145 		| t4[0] | t4[1] | t4[2] | t4[3] | t4[4];
1146 	tt = ((tt | -tt) >> 63) - (uint64_t)1;
1147 
1148 	/*
1149 	 * Compute u1*h^2 (in t6) and h^3 (in t5);
1150 	 */
1151 	f256_montysquare(t7, t2);
1152 	f256_montymul(t6, t1, t7);
1153 	f256_montymul(t5, t7, t2);
1154 
1155 	/*
1156 	 * Compute x3 = r^2 - h^3 - 2*u1*h^2.
1157 	 */
1158 	f256_montysquare(P1->x, t4);
1159 	f256_sub(P1->x, P1->x, t5);
1160 	f256_sub(P1->x, P1->x, t6);
1161 	f256_sub(P1->x, P1->x, t6);
1162 
1163 	/*
1164 	 * Compute y3 = r*(u1*h^2 - x3) - s1*h^3.
1165 	 */
1166 	f256_sub(t6, t6, P1->x);
1167 	f256_montymul(P1->y, t4, t6);
1168 	f256_montymul(t1, t5, t3);
1169 	f256_sub(P1->y, P1->y, t1);
1170 
1171 	/*
1172 	 * Compute z3 = h*z1.
1173 	 */
1174 	f256_montymul(P1->z, P1->z, t2);
1175 
1176 	/*
1177 	 * The "double" result, in case P1 = P2.
1178 	 */
1179 
1180 	/*
1181 	 * Compute z' = 2*y2 (in t1).
1182 	 */
1183 	f256_add(t1, P2->y, P2->y);
1184 	f256_partial_reduce(t1);
1185 
1186 	/*
1187 	 * Compute 2*(y2^2) (in t2) and s = 4*x2*(y2^2) (in t3).
1188 	 */
1189 	f256_montysquare(t2, P2->y);
1190 	f256_add(t2, t2, t2);
1191 	f256_add(t3, t2, t2);
1192 	f256_montymul(t3, P2->x, t3);
1193 
1194 	/*
1195 	 * Compute m = 3*(x2^2 - 1) (in t4).
1196 	 */
1197 	f256_montysquare(t4, P2->x);
1198 	f256_sub(t4, t4, F256_R);
1199 	f256_add(t5, t4, t4);
1200 	f256_add(t4, t4, t5);
1201 
1202 	/*
1203 	 * Compute x' = m^2 - 2*s (in t5).
1204 	 */
1205 	f256_montysquare(t5, t4);
1206 	f256_sub(t5, t3);
1207 	f256_sub(t5, t3);
1208 
1209 	/*
1210 	 * Compute y' = m*(s - x') - 8*y2^4 (in t6).
1211 	 */
1212 	f256_sub(t6, t3, t5);
1213 	f256_montymul(t6, t6, t4);
1214 	f256_montysquare(t7, t2);
1215 	f256_sub(t6, t6, t7);
1216 	f256_sub(t6, t6, t7);
1217 
1218 	/*
1219 	 * We now have the alternate (doubling) coordinates in (t5,t6,t1).
1220 	 * We combine them with (x3,y3,z3).
1221 	 */
1222 	for (i = 0; i < 5; i ++) {
1223 		P1->x[i] |= tt & t5[i];
1224 		P1->y[i] |= tt & t6[i];
1225 		P1->z[i] |= tt & t1[i];
1226 	}
1227 
1228 	/*
1229 	 * If P1 = 0, then we get z3 = 0 (which is invalid); if z1 is 0,
1230 	 * then we want to replace the result with a copy of P2. The
1231 	 * test on z1 was done at the start, in the zz mask.
1232 	 */
1233 	for (i = 0; i < 5; i ++) {
1234 		P1->x[i] ^= zz & (P1->x[i] ^ P2->x[i]);
1235 		P1->y[i] ^= zz & (P1->y[i] ^ P2->y[i]);
1236 		P1->z[i] ^= zz & (P1->z[i] ^ F256_R[i]);
1237 	}
1238 }
1239 #endif
1240 
1241 /*
1242  * Inner function for computing a point multiplication. A window is
1243  * provided, with points 1*P to 15*P in affine coordinates.
1244  *
1245  * Assumptions:
1246  *  - All provided points are valid points on the curve.
1247  *  - Multiplier is non-zero, and smaller than the curve order.
1248  *  - Everything is in Montgomery representation.
1249  */
1250 static void
1251 point_mul_inner(p256_jacobian *R, const p256_affine *W,
1252 	const unsigned char *k, size_t klen)
1253 {
1254 	p256_jacobian Q;
1255 	uint32_t qz;
1256 
1257 	memset(&Q, 0, sizeof Q);
1258 	qz = 1;
1259 	while (klen -- > 0) {
1260 		int i;
1261 		unsigned bk;
1262 
1263 		bk = *k ++;
1264 		for (i = 0; i < 2; i ++) {
1265 			uint32_t bits;
1266 			uint32_t bnz;
1267 			p256_affine T;
1268 			p256_jacobian U;
1269 			uint32_t n;
1270 			int j;
1271 			uint64_t m;
1272 
1273 			p256_double(&Q);
1274 			p256_double(&Q);
1275 			p256_double(&Q);
1276 			p256_double(&Q);
1277 			bits = (bk >> 4) & 0x0F;
1278 			bnz = NEQ(bits, 0);
1279 
1280 			/*
1281 			 * Lookup point in window. If the bits are 0,
1282 			 * we get something invalid, which is not a
1283 			 * problem because we will use it only if the
1284 			 * bits are non-zero.
1285 			 */
1286 			memset(&T, 0, sizeof T);
1287 			for (n = 0; n < 15; n ++) {
1288 				m = -(uint64_t)EQ(bits, n + 1);
1289 				T.x[0] |= m & W[n].x[0];
1290 				T.x[1] |= m & W[n].x[1];
1291 				T.x[2] |= m & W[n].x[2];
1292 				T.x[3] |= m & W[n].x[3];
1293 				T.x[4] |= m & W[n].x[4];
1294 				T.y[0] |= m & W[n].y[0];
1295 				T.y[1] |= m & W[n].y[1];
1296 				T.y[2] |= m & W[n].y[2];
1297 				T.y[3] |= m & W[n].y[3];
1298 				T.y[4] |= m & W[n].y[4];
1299 			}
1300 
1301 			U = Q;
1302 			p256_add_mixed(&U, &T);
1303 
1304 			/*
1305 			 * If qz is still 1, then Q was all-zeros, and this
1306 			 * is conserved through p256_double().
1307 			 */
1308 			m = -(uint64_t)(bnz & qz);
1309 			for (j = 0; j < 5; j ++) {
1310 				Q.x[j] ^= m & (Q.x[j] ^ T.x[j]);
1311 				Q.y[j] ^= m & (Q.y[j] ^ T.y[j]);
1312 				Q.z[j] ^= m & (Q.z[j] ^ F256_R[j]);
1313 			}
1314 			CCOPY(bnz & ~qz, &Q, &U, sizeof Q);
1315 			qz &= ~bnz;
1316 			bk <<= 4;
1317 		}
1318 	}
1319 	*R = Q;
1320 }
1321 
1322 /*
1323  * Convert a window from Jacobian to affine coordinates. A single
1324  * field inversion is used. This function works for windows up to
1325  * 32 elements.
1326  *
1327  * The destination array (aff[]) and the source array (jac[]) may
1328  * overlap, provided that the start of aff[] is not after the start of
1329  * jac[]. Even if the arrays do _not_ overlap, the source array is
1330  * modified.
1331  */
1332 static void
1333 window_to_affine(p256_affine *aff, p256_jacobian *jac, int num)
1334 {
1335 	/*
1336 	 * Convert the window points to affine coordinates. We use the
1337 	 * following trick to mutualize the inversion computation: if
1338 	 * we have z1, z2, z3, and z4, and want to invert all of them,
1339 	 * we compute u = 1/(z1*z2*z3*z4), and then we have:
1340 	 *   1/z1 = u*z2*z3*z4
1341 	 *   1/z2 = u*z1*z3*z4
1342 	 *   1/z3 = u*z1*z2*z4
1343 	 *   1/z4 = u*z1*z2*z3
1344 	 *
1345 	 * The partial products are computed recursively:
1346 	 *
1347 	 *  - on input (z_1,z_2), return (z_2,z_1) and z_1*z_2
1348 	 *  - on input (z_1,z_2,... z_n):
1349 	 *       recurse on (z_1,z_2,... z_(n/2)) -> r1 and m1
1350 	 *       recurse on (z_(n/2+1),z_(n/2+2)... z_n) -> r2 and m2
1351 	 *       multiply elements of r1 by m2 -> s1
1352 	 *       multiply elements of r2 by m1 -> s2
1353 	 *       return r1||r2 and m1*m2
1354 	 *
1355 	 * In the example below, we suppose that we have 14 elements.
1356 	 * Let z1, z2,... zE be the 14 values to invert (index noted in
1357 	 * hexadecimal, starting at 1).
1358 	 *
1359 	 *  - Depth 1:
1360 	 *      swap(z1, z2); z12 = z1*z2
1361 	 *      swap(z3, z4); z34 = z3*z4
1362 	 *      swap(z5, z6); z56 = z5*z6
1363 	 *      swap(z7, z8); z78 = z7*z8
1364 	 *      swap(z9, zA); z9A = z9*zA
1365 	 *      swap(zB, zC); zBC = zB*zC
1366 	 *      swap(zD, zE); zDE = zD*zE
1367 	 *
1368 	 *  - Depth 2:
1369 	 *      z1 <- z1*z34, z2 <- z2*z34, z3 <- z3*z12, z4 <- z4*z12
1370 	 *      z1234 = z12*z34
1371 	 *      z5 <- z5*z78, z6 <- z6*z78, z7 <- z7*z56, z8 <- z8*z56
1372 	 *      z5678 = z56*z78
1373 	 *      z9 <- z9*zBC, zA <- zA*zBC, zB <- zB*z9A, zC <- zC*z9A
1374 	 *      z9ABC = z9A*zBC
1375 	 *
1376 	 *  - Depth 3:
1377 	 *      z1 <- z1*z5678, z2 <- z2*z5678, z3 <- z3*z5678, z4 <- z4*z5678
1378 	 *      z5 <- z5*z1234, z6 <- z6*z1234, z7 <- z7*z1234, z8 <- z8*z1234
1379 	 *      z12345678 = z1234*z5678
1380 	 *      z9 <- z9*zDE, zA <- zA*zDE, zB <- zB*zDE, zC <- zC*zDE
1381 	 *      zD <- zD*z9ABC, zE*z9ABC
1382 	 *      z9ABCDE = z9ABC*zDE
1383 	 *
1384 	 *  - Depth 4:
1385 	 *      multiply z1..z8 by z9ABCDE
1386 	 *      multiply z9..zE by z12345678
1387 	 *      final z = z12345678*z9ABCDE
1388 	 */
1389 
1390 	uint64_t z[16][5];
1391 	int i, k, s;
1392 #define zt   (z[15])
1393 #define zu   (z[14])
1394 #define zv   (z[13])
1395 
1396 	/*
1397 	 * First recursion step (pairwise swapping and multiplication).
1398 	 * If there is an odd number of elements, then we "invent" an
1399 	 * extra one with coordinate Z = 1 (in Montgomery representation).
1400 	 */
1401 	for (i = 0; (i + 1) < num; i += 2) {
1402 		memcpy(zt, jac[i].z, sizeof zt);
1403 		memcpy(jac[i].z, jac[i + 1].z, sizeof zt);
1404 		memcpy(jac[i + 1].z, zt, sizeof zt);
1405 		f256_montymul(z[i >> 1], jac[i].z, jac[i + 1].z);
1406 	}
1407 	if ((num & 1) != 0) {
1408 		memcpy(z[num >> 1], jac[num - 1].z, sizeof zt);
1409 		memcpy(jac[num - 1].z, F256_R, sizeof F256_R);
1410 	}
1411 
1412 	/*
1413 	 * Perform further recursion steps. At the entry of each step,
1414 	 * the process has been done for groups of 's' points. The
1415 	 * integer k is the log2 of s.
1416 	 */
1417 	for (k = 1, s = 2; s < num; k ++, s <<= 1) {
1418 		int n;
1419 
1420 		for (i = 0; i < num; i ++) {
1421 			f256_montymul(jac[i].z, jac[i].z, z[(i >> k) ^ 1]);
1422 		}
1423 		n = (num + s - 1) >> k;
1424 		for (i = 0; i < (n >> 1); i ++) {
1425 			f256_montymul(z[i], z[i << 1], z[(i << 1) + 1]);
1426 		}
1427 		if ((n & 1) != 0) {
1428 			memmove(z[n >> 1], z[n], sizeof zt);
1429 		}
1430 	}
1431 
1432 	/*
1433 	 * Invert the final result, and convert all points.
1434 	 */
1435 	f256_invert(zt, z[0]);
1436 	for (i = 0; i < num; i ++) {
1437 		f256_montymul(zv, jac[i].z, zt);
1438 		f256_montysquare(zu, zv);
1439 		f256_montymul(zv, zv, zu);
1440 		f256_montymul(aff[i].x, jac[i].x, zu);
1441 		f256_montymul(aff[i].y, jac[i].y, zv);
1442 	}
1443 }
1444 
1445 /*
1446  * Multiply the provided point by an integer.
1447  * Assumptions:
1448  *  - Source point is a valid curve point.
1449  *  - Source point is not the point-at-infinity.
1450  *  - Integer is not 0, and is lower than the curve order.
1451  * If these conditions are not met, then the result is indeterminate
1452  * (but the process is still constant-time).
1453  */
1454 static void
1455 p256_mul(p256_jacobian *P, const unsigned char *k, size_t klen)
1456 {
1457 	union {
1458 		p256_affine aff[15];
1459 		p256_jacobian jac[15];
1460 	} window;
1461 	int i;
1462 
1463 	/*
1464 	 * Compute window, in Jacobian coordinates.
1465 	 */
1466 	window.jac[0] = *P;
1467 	for (i = 2; i < 16; i ++) {
1468 		window.jac[i - 1] = window.jac[(i >> 1) - 1];
1469 		if ((i & 1) == 0) {
1470 			p256_double(&window.jac[i - 1]);
1471 		} else {
1472 			p256_add(&window.jac[i - 1], &window.jac[i >> 1]);
1473 		}
1474 	}
1475 
1476 	/*
1477 	 * Convert the window points to affine coordinates. Point
1478 	 * window[0] is the source point, already in affine coordinates.
1479 	 */
1480 	window_to_affine(window.aff, window.jac, 15);
1481 
1482 	/*
1483 	 * Perform point multiplication.
1484 	 */
1485 	point_mul_inner(P, window.aff, k, klen);
1486 }
1487 
1488 /*
1489  * Precomputed window for the conventional generator: P256_Gwin[n]
1490  * contains (n+1)*G (affine coordinates, in Montgomery representation).
1491  */
1492 static const p256_affine P256_Gwin[] = {
1493 	{
1494 		{ 0x30D418A9143C1, 0xC4FEDB60179E7, 0x62251075BA95F,
1495 		  0x5C669FB732B77, 0x08905F76B5375 },
1496 		{ 0x5357CE95560A8, 0x43A19E45CDDF2, 0x21F3258B4AB8E,
1497 		  0xD8552E88688DD, 0x0571FF18A5885 }
1498 	},
1499 	{
1500 		{ 0x46D410DDD64DF, 0x0B433827D8500, 0x1490D9AA6AE3C,
1501 		  0xA3A832205038D, 0x06BB32E52DCF3 },
1502 		{ 0x48D361BEE1A57, 0xB7B236FF82F36, 0x042DBE152CD7C,
1503 		  0xA3AA9A8FB0E92, 0x08C577517A5B8 }
1504 	},
1505 	{
1506 		{ 0x3F904EEBC1272, 0x9E87D81FBFFAC, 0xCBBC98B027F84,
1507 		  0x47E46AD77DD87, 0x06936A3FD6FF7 },
1508 		{ 0x5C1FC983A7EBD, 0xC3861FE1AB04C, 0x2EE98E583E47A,
1509 		  0xC06A88208311A, 0x05F06A2AB587C }
1510 	},
1511 	{
1512 		{ 0xB50D46918DCC5, 0xD7623C17374B0, 0x100AF24650A6E,
1513 		  0x76ABCDAACACE8, 0x077362F591B01 },
1514 		{ 0xF24CE4CBABA68, 0x17AD6F4472D96, 0xDDD22E1762847,
1515 		  0x862EB6C36DEE5, 0x04B14C39CC5AB }
1516 	},
1517 	{
1518 		{ 0x8AAEC45C61F5C, 0x9D4B9537DBE1B, 0x76C20C90EC649,
1519 		  0x3C7D41CB5AAD0, 0x0907960649052 },
1520 		{ 0x9B4AE7BA4F107, 0xF75EB882BEB30, 0x7A1F6873C568E,
1521 		  0x915C540A9877E, 0x03A076BB9DD1E }
1522 	},
1523 	{
1524 		{ 0x47373E77664A1, 0xF246CEE3E4039, 0x17A3AD55AE744,
1525 		  0x673C50A961A5B, 0x03074B5964213 },
1526 		{ 0x6220D377E44BA, 0x30DFF14B593D3, 0x639F11299C2B5,
1527 		  0x75F5424D44CEF, 0x04C9916DEA07F }
1528 	},
1529 	{
1530 		{ 0x354EA0173B4F1, 0x3C23C00F70746, 0x23BB082BD2021,
1531 		  0xE03E43EAAB50C, 0x03BA5119D3123 },
1532 		{ 0xD0303F5B9D4DE, 0x17DA67BDD2847, 0xC941956742F2F,
1533 		  0x8670F933BDC77, 0x0AEDD9164E240 }
1534 	},
1535 	{
1536 		{ 0x4CD19499A78FB, 0x4BF9B345527F1, 0x2CFC6B462AB5C,
1537 		  0x30CDF90F02AF0, 0x0763891F62652 },
1538 		{ 0xA3A9532D49775, 0xD7F9EBA15F59D, 0x60BBF021E3327,
1539 		  0xF75C23C7B84BE, 0x06EC12F2C706D }
1540 	},
1541 	{
1542 		{ 0x6E8F264E20E8E, 0xC79A7A84175C9, 0xC8EB00ABE6BFE,
1543 		  0x16A4CC09C0444, 0x005B3081D0C4E },
1544 		{ 0x777AA45F33140, 0xDCE5D45E31EB7, 0xB12F1A56AF7BE,
1545 		  0xF9B2B6E019A88, 0x086659CDFD835 }
1546 	},
1547 	{
1548 		{ 0xDBD19DC21EC8C, 0x94FCF81392C18, 0x250B4998F9868,
1549 		  0x28EB37D2CD648, 0x0C61C947E4B34 },
1550 		{ 0x407880DD9E767, 0x0C83FBE080C2B, 0x9BE5D2C43A899,
1551 		  0xAB4EF7D2D6577, 0x08719A555B3B4 }
1552 	},
1553 	{
1554 		{ 0x260A6245E4043, 0x53E7FDFE0EA7D, 0xAC1AB59DE4079,
1555 		  0x072EFF3A4158D, 0x0E7090F1949C9 },
1556 		{ 0x85612B944E886, 0xE857F61C81A76, 0xAD643D250F939,
1557 		  0x88DAC0DAA891E, 0x089300244125B }
1558 	},
1559 	{
1560 		{ 0x1AA7D26977684, 0x58A345A3304B7, 0x37385EABDEDEF,
1561 		  0x155E409D29DEE, 0x0EE1DF780B83E },
1562 		{ 0x12D91CBB5B437, 0x65A8956370CAC, 0xDE6D66170ED2F,
1563 		  0xAC9B8228CFA8A, 0x0FF57C95C3238 }
1564 	},
1565 	{
1566 		{ 0x25634B2ED7097, 0x9156FD30DCCC4, 0x9E98110E35676,
1567 		  0x7594CBCD43F55, 0x038477ACC395B },
1568 		{ 0x2B90C00EE17FF, 0xF842ED2E33575, 0x1F5BC16874838,
1569 		  0x7968CD06422BD, 0x0BC0876AB9E7B }
1570 	},
1571 	{
1572 		{ 0xA35BB0CF664AF, 0x68F9707E3A242, 0x832660126E48F,
1573 		  0x72D2717BF54C6, 0x0AAE7333ED12C },
1574 		{ 0x2DB7995D586B1, 0xE732237C227B5, 0x65E7DBBE29569,
1575 		  0xBBBD8E4193E2A, 0x052706DC3EAA1 }
1576 	},
1577 	{
1578 		{ 0xD8B7BC60055BE, 0xD76E27E4B72BC, 0x81937003CC23E,
1579 		  0xA090E337424E4, 0x02AA0E43EAD3D },
1580 		{ 0x524F6383C45D2, 0x422A41B2540B8, 0x8A4797D766355,
1581 		  0xDF444EFA6DE77, 0x0042170A9079A }
1582 	},
1583 };
1584 
1585 /*
1586  * Multiply the conventional generator of the curve by the provided
1587  * integer. Return is written in *P.
1588  *
1589  * Assumptions:
1590  *  - Integer is not 0, and is lower than the curve order.
1591  * If this conditions is not met, then the result is indeterminate
1592  * (but the process is still constant-time).
1593  */
1594 static void
1595 p256_mulgen(p256_jacobian *P, const unsigned char *k, size_t klen)
1596 {
1597 	point_mul_inner(P, P256_Gwin, k, klen);
1598 }
1599 
1600 /*
1601  * Return 1 if all of the following hold:
1602  *  - klen <= 32
1603  *  - k != 0
1604  *  - k is lower than the curve order
1605  * Otherwise, return 0.
1606  *
1607  * Constant-time behaviour: only klen may be observable.
1608  */
1609 static uint32_t
1610 check_scalar(const unsigned char *k, size_t klen)
1611 {
1612 	uint32_t z;
1613 	int32_t c;
1614 	size_t u;
1615 
1616 	if (klen > 32) {
1617 		return 0;
1618 	}
1619 	z = 0;
1620 	for (u = 0; u < klen; u ++) {
1621 		z |= k[u];
1622 	}
1623 	if (klen == 32) {
1624 		c = 0;
1625 		for (u = 0; u < klen; u ++) {
1626 			c |= -(int32_t)EQ0(c) & CMP(k[u], P256_N[u]);
1627 		}
1628 	} else {
1629 		c = -1;
1630 	}
1631 	return NEQ(z, 0) & LT0(c);
1632 }
1633 
1634 static uint32_t
1635 api_mul(unsigned char *G, size_t Glen,
1636 	const unsigned char *k, size_t klen, int curve)
1637 {
1638 	uint32_t r;
1639 	p256_jacobian P;
1640 
1641 	(void)curve;
1642 	if (Glen != 65) {
1643 		return 0;
1644 	}
1645 	r = check_scalar(k, klen);
1646 	r &= point_decode(&P, G);
1647 	p256_mul(&P, k, klen);
1648 	r &= point_encode(G, &P);
1649 	return r;
1650 }
1651 
1652 static size_t
1653 api_mulgen(unsigned char *R,
1654 	const unsigned char *k, size_t klen, int curve)
1655 {
1656 	p256_jacobian P;
1657 
1658 	(void)curve;
1659 	p256_mulgen(&P, k, klen);
1660 	point_encode(R, &P);
1661 	return 65;
1662 }
1663 
1664 static uint32_t
1665 api_muladd(unsigned char *A, const unsigned char *B, size_t len,
1666 	const unsigned char *x, size_t xlen,
1667 	const unsigned char *y, size_t ylen, int curve)
1668 {
1669 	/*
1670 	 * We might want to use Shamir's trick here: make a composite
1671 	 * window of u*P+v*Q points, to merge the two doubling-ladders
1672 	 * into one. This, however, has some complications:
1673 	 *
1674 	 *  - During the computation, we may hit the point-at-infinity.
1675 	 *    Thus, we would need p256_add_complete_mixed() (complete
1676 	 *    formulas for point addition), with a higher cost (17 muls
1677 	 *    instead of 11).
1678 	 *
1679 	 *  - A 4-bit window would be too large, since it would involve
1680 	 *    16*16-1 = 255 points. For the same window size as in the
1681 	 *    p256_mul() case, we would need to reduce the window size
1682 	 *    to 2 bits, and thus perform twice as many non-doubling
1683 	 *    point additions.
1684 	 *
1685 	 *  - The window may itself contain the point-at-infinity, and
1686 	 *    thus cannot be in all generality be made of affine points.
1687 	 *    Instead, we would need to make it a window of points in
1688 	 *    Jacobian coordinates. Even p256_add_complete_mixed() would
1689 	 *    be inappropriate.
1690 	 *
1691 	 * For these reasons, the code below performs two separate
1692 	 * point multiplications, then computes the final point addition
1693 	 * (which is both a "normal" addition, and a doubling, to handle
1694 	 * all cases).
1695 	 */
1696 
1697 	p256_jacobian P, Q;
1698 	uint32_t r, t, s;
1699 	uint64_t z;
1700 
1701 	(void)curve;
1702 	if (len != 65) {
1703 		return 0;
1704 	}
1705 	r = point_decode(&P, A);
1706 	p256_mul(&P, x, xlen);
1707 	if (B == NULL) {
1708 		p256_mulgen(&Q, y, ylen);
1709 	} else {
1710 		r &= point_decode(&Q, B);
1711 		p256_mul(&Q, y, ylen);
1712 	}
1713 
1714 	/*
1715 	 * The final addition may fail in case both points are equal.
1716 	 */
1717 	t = p256_add(&P, &Q);
1718 	f256_final_reduce(P.z);
1719 	z = P.z[0] | P.z[1] | P.z[2] | P.z[3] | P.z[4];
1720 	s = EQ((uint32_t)(z | (z >> 32)), 0);
1721 	p256_double(&Q);
1722 
1723 	/*
1724 	 * If s is 1 then either P+Q = 0 (t = 1) or P = Q (t = 0). So we
1725 	 * have the following:
1726 	 *
1727 	 *   s = 0, t = 0   return P (normal addition)
1728 	 *   s = 0, t = 1   return P (normal addition)
1729 	 *   s = 1, t = 0   return Q (a 'double' case)
1730 	 *   s = 1, t = 1   report an error (P+Q = 0)
1731 	 */
1732 	CCOPY(s & ~t, &P, &Q, sizeof Q);
1733 	point_encode(A, &P);
1734 	r &= ~(s & t);
1735 	return r;
1736 }
1737 
1738 /* see bearssl_ec.h */
1739 const br_ec_impl br_ec_p256_m62 = {
1740 	(uint32_t)0x00800000,
1741 	&api_generator,
1742 	&api_order,
1743 	&api_xoff,
1744 	&api_mul,
1745 	&api_mulgen,
1746 	&api_muladd
1747 };
1748 
1749 /* see bearssl_ec.h */
1750 const br_ec_impl *
1751 br_ec_p256_m62_get(void)
1752 {
1753 	return &br_ec_p256_m62;
1754 }
1755 
1756 #else
1757 
1758 /* see bearssl_ec.h */
1759 const br_ec_impl *
1760 br_ec_p256_m62_get(void)
1761 {
1762 	return 0;
1763 }
1764 
1765 #endif
1766