xref: /freebsd/contrib/bearssl/src/ec/ec_c25519_m31.c (revision cfd6422a5217410fbd66f7a7a8a64d9d85e61229)
1 /*
2  * Copyright (c) 2017 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 /* obsolete
28 #include <stdio.h>
29 #include <stdlib.h>
30 static void
31 print_int(const char *name, const uint32_t *x)
32 {
33 	size_t u;
34 	unsigned char tmp[40];
35 
36 	printf("%s = ", name);
37 	for (u = 0; u < 9; u ++) {
38 		if (x[u] > 0x3FFFFFFF) {
39 			printf("INVALID:");
40 			for (u = 0; u < 9; u ++) {
41 				printf(" %08X", x[u]);
42 			}
43 			printf("\n");
44 			return;
45 		}
46 	}
47 	memset(tmp, 0, sizeof tmp);
48 	for (u = 0; u < 9; u ++) {
49 		uint64_t w;
50 		int j, k;
51 
52 		w = x[u];
53 		j = 30 * (int)u;
54 		k = j & 7;
55 		if (k != 0) {
56 			w <<= k;
57 			j -= k;
58 		}
59 		k = j >> 3;
60 		for (j = 0; j < 8; j ++) {
61 			tmp[39 - k - j] |= (unsigned char)w;
62 			w >>= 8;
63 		}
64 	}
65 	for (u = 8; u < 40; u ++) {
66 		printf("%02X", tmp[u]);
67 	}
68 	printf("\n");
69 }
70 */
71 
72 /*
73  * If BR_NO_ARITH_SHIFT is undefined, or defined to 0, then we _assume_
74  * that right-shifting a signed negative integer copies the sign bit
75  * (arithmetic right-shift). This is "implementation-defined behaviour",
76  * i.e. it is not undefined, but it may differ between compilers. Each
77  * compiler is supposed to document its behaviour in that respect. GCC
78  * explicitly defines that an arithmetic right shift is used. We expect
79  * all other compilers to do the same, because underlying CPU offer an
80  * arithmetic right shift opcode that could not be used otherwise.
81  */
82 #if BR_NO_ARITH_SHIFT
83 #define ARSH(x, n)   (((uint32_t)(x) >> (n)) \
84                     | ((-((uint32_t)(x) >> 31)) << (32 - (n))))
85 #else
86 #define ARSH(x, n)   ((*(int32_t *)&(x)) >> (n))
87 #endif
88 
89 /*
90  * Convert an integer from unsigned little-endian encoding to a sequence of
91  * 30-bit words in little-endian order. The final "partial" word is
92  * returned.
93  */
94 static uint32_t
95 le8_to_le30(uint32_t *dst, const unsigned char *src, size_t len)
96 {
97 	uint32_t acc;
98 	int acc_len;
99 
100 	acc = 0;
101 	acc_len = 0;
102 	while (len -- > 0) {
103 		uint32_t b;
104 
105 		b = *src ++;
106 		if (acc_len < 22) {
107 			acc |= b << acc_len;
108 			acc_len += 8;
109 		} else {
110 			*dst ++ = (acc | (b << acc_len)) & 0x3FFFFFFF;
111 			acc = b >> (30 - acc_len);
112 			acc_len -= 22;
113 		}
114 	}
115 	return acc;
116 }
117 
118 /*
119  * Convert an integer (30-bit words, little-endian) to unsigned
120  * little-endian encoding. The total encoding length is provided; all
121  * the destination bytes will be filled.
122  */
123 static void
124 le30_to_le8(unsigned char *dst, size_t len, const uint32_t *src)
125 {
126 	uint32_t acc;
127 	int acc_len;
128 
129 	acc = 0;
130 	acc_len = 0;
131 	while (len -- > 0) {
132 		if (acc_len < 8) {
133 			uint32_t w;
134 
135 			w = *src ++;
136 			*dst ++ = (unsigned char)(acc | (w << acc_len));
137 			acc = w >> (8 - acc_len);
138 			acc_len += 22;
139 		} else {
140 			*dst ++ = (unsigned char)acc;
141 			acc >>= 8;
142 			acc_len -= 8;
143 		}
144 	}
145 }
146 
147 /*
148  * Multiply two integers. Source integers are represented as arrays of
149  * nine 30-bit words, for values up to 2^270-1. Result is encoded over
150  * 18 words of 30 bits each.
151  */
152 static void
153 mul9(uint32_t *d, const uint32_t *a, const uint32_t *b)
154 {
155 	/*
156 	 * Maximum intermediate result is no more than
157 	 * 10376293531797946367, which fits in 64 bits. Reason:
158 	 *
159 	 *   10376293531797946367 = 9 * (2^30-1)^2 + 9663676406
160 	 *   10376293531797946367 < 9663676407 * 2^30
161 	 *
162 	 * Thus, adding together 9 products of 30-bit integers, with
163 	 * a carry of at most 9663676406, yields an integer that fits
164 	 * on 64 bits and generates a carry of at most 9663676406.
165 	 */
166 	uint64_t t[17];
167 	uint64_t cc;
168 	int i;
169 
170 	t[ 0] = MUL31(a[0], b[0]);
171 	t[ 1] = MUL31(a[0], b[1])
172 		+ MUL31(a[1], b[0]);
173 	t[ 2] = MUL31(a[0], b[2])
174 		+ MUL31(a[1], b[1])
175 		+ MUL31(a[2], b[0]);
176 	t[ 3] = MUL31(a[0], b[3])
177 		+ MUL31(a[1], b[2])
178 		+ MUL31(a[2], b[1])
179 		+ MUL31(a[3], b[0]);
180 	t[ 4] = MUL31(a[0], b[4])
181 		+ MUL31(a[1], b[3])
182 		+ MUL31(a[2], b[2])
183 		+ MUL31(a[3], b[1])
184 		+ MUL31(a[4], b[0]);
185 	t[ 5] = MUL31(a[0], b[5])
186 		+ MUL31(a[1], b[4])
187 		+ MUL31(a[2], b[3])
188 		+ MUL31(a[3], b[2])
189 		+ MUL31(a[4], b[1])
190 		+ MUL31(a[5], b[0]);
191 	t[ 6] = MUL31(a[0], b[6])
192 		+ MUL31(a[1], b[5])
193 		+ MUL31(a[2], b[4])
194 		+ MUL31(a[3], b[3])
195 		+ MUL31(a[4], b[2])
196 		+ MUL31(a[5], b[1])
197 		+ MUL31(a[6], b[0]);
198 	t[ 7] = MUL31(a[0], b[7])
199 		+ MUL31(a[1], b[6])
200 		+ MUL31(a[2], b[5])
201 		+ MUL31(a[3], b[4])
202 		+ MUL31(a[4], b[3])
203 		+ MUL31(a[5], b[2])
204 		+ MUL31(a[6], b[1])
205 		+ MUL31(a[7], b[0]);
206 	t[ 8] = MUL31(a[0], b[8])
207 		+ MUL31(a[1], b[7])
208 		+ MUL31(a[2], b[6])
209 		+ MUL31(a[3], b[5])
210 		+ MUL31(a[4], b[4])
211 		+ MUL31(a[5], b[3])
212 		+ MUL31(a[6], b[2])
213 		+ MUL31(a[7], b[1])
214 		+ MUL31(a[8], b[0]);
215 	t[ 9] = MUL31(a[1], b[8])
216 		+ MUL31(a[2], b[7])
217 		+ MUL31(a[3], b[6])
218 		+ MUL31(a[4], b[5])
219 		+ MUL31(a[5], b[4])
220 		+ MUL31(a[6], b[3])
221 		+ MUL31(a[7], b[2])
222 		+ MUL31(a[8], b[1]);
223 	t[10] = MUL31(a[2], b[8])
224 		+ MUL31(a[3], b[7])
225 		+ MUL31(a[4], b[6])
226 		+ MUL31(a[5], b[5])
227 		+ MUL31(a[6], b[4])
228 		+ MUL31(a[7], b[3])
229 		+ MUL31(a[8], b[2]);
230 	t[11] = MUL31(a[3], b[8])
231 		+ MUL31(a[4], b[7])
232 		+ MUL31(a[5], b[6])
233 		+ MUL31(a[6], b[5])
234 		+ MUL31(a[7], b[4])
235 		+ MUL31(a[8], b[3]);
236 	t[12] = MUL31(a[4], b[8])
237 		+ MUL31(a[5], b[7])
238 		+ MUL31(a[6], b[6])
239 		+ MUL31(a[7], b[5])
240 		+ MUL31(a[8], b[4]);
241 	t[13] = MUL31(a[5], b[8])
242 		+ MUL31(a[6], b[7])
243 		+ MUL31(a[7], b[6])
244 		+ MUL31(a[8], b[5]);
245 	t[14] = MUL31(a[6], b[8])
246 		+ MUL31(a[7], b[7])
247 		+ MUL31(a[8], b[6]);
248 	t[15] = MUL31(a[7], b[8])
249 		+ MUL31(a[8], b[7]);
250 	t[16] = MUL31(a[8], b[8]);
251 
252 	/*
253 	 * Propagate carries.
254 	 */
255 	cc = 0;
256 	for (i = 0; i < 17; i ++) {
257 		uint64_t w;
258 
259 		w = t[i] + cc;
260 		d[i] = (uint32_t)w & 0x3FFFFFFF;
261 		cc = w >> 30;
262 	}
263 	d[17] = (uint32_t)cc;
264 }
265 
266 /*
267  * Square a 270-bit integer, represented as an array of nine 30-bit words.
268  * Result uses 18 words of 30 bits each.
269  */
270 static void
271 square9(uint32_t *d, const uint32_t *a)
272 {
273 	uint64_t t[17];
274 	uint64_t cc;
275 	int i;
276 
277 	t[ 0] = MUL31(a[0], a[0]);
278 	t[ 1] = ((MUL31(a[0], a[1])) << 1);
279 	t[ 2] = MUL31(a[1], a[1])
280 		+ ((MUL31(a[0], a[2])) << 1);
281 	t[ 3] = ((MUL31(a[0], a[3])
282 		+ MUL31(a[1], a[2])) << 1);
283 	t[ 4] = MUL31(a[2], a[2])
284 		+ ((MUL31(a[0], a[4])
285 		+ MUL31(a[1], a[3])) << 1);
286 	t[ 5] = ((MUL31(a[0], a[5])
287 		+ MUL31(a[1], a[4])
288 		+ MUL31(a[2], a[3])) << 1);
289 	t[ 6] = MUL31(a[3], a[3])
290 		+ ((MUL31(a[0], a[6])
291 		+ MUL31(a[1], a[5])
292 		+ MUL31(a[2], a[4])) << 1);
293 	t[ 7] = ((MUL31(a[0], a[7])
294 		+ MUL31(a[1], a[6])
295 		+ MUL31(a[2], a[5])
296 		+ MUL31(a[3], a[4])) << 1);
297 	t[ 8] = MUL31(a[4], a[4])
298 		+ ((MUL31(a[0], a[8])
299 		+ MUL31(a[1], a[7])
300 		+ MUL31(a[2], a[6])
301 		+ MUL31(a[3], a[5])) << 1);
302 	t[ 9] = ((MUL31(a[1], a[8])
303 		+ MUL31(a[2], a[7])
304 		+ MUL31(a[3], a[6])
305 		+ MUL31(a[4], a[5])) << 1);
306 	t[10] = MUL31(a[5], a[5])
307 		+ ((MUL31(a[2], a[8])
308 		+ MUL31(a[3], a[7])
309 		+ MUL31(a[4], a[6])) << 1);
310 	t[11] = ((MUL31(a[3], a[8])
311 		+ MUL31(a[4], a[7])
312 		+ MUL31(a[5], a[6])) << 1);
313 	t[12] = MUL31(a[6], a[6])
314 		+ ((MUL31(a[4], a[8])
315 		+ MUL31(a[5], a[7])) << 1);
316 	t[13] = ((MUL31(a[5], a[8])
317 		+ MUL31(a[6], a[7])) << 1);
318 	t[14] = MUL31(a[7], a[7])
319 		+ ((MUL31(a[6], a[8])) << 1);
320 	t[15] = ((MUL31(a[7], a[8])) << 1);
321 	t[16] = MUL31(a[8], a[8]);
322 
323 	/*
324 	 * Propagate carries.
325 	 */
326 	cc = 0;
327 	for (i = 0; i < 17; i ++) {
328 		uint64_t w;
329 
330 		w = t[i] + cc;
331 		d[i] = (uint32_t)w & 0x3FFFFFFF;
332 		cc = w >> 30;
333 	}
334 	d[17] = (uint32_t)cc;
335 }
336 
337 /*
338  * Perform a "final reduction" in field F255 (field for Curve25519)
339  * The source value must be less than twice the modulus. If the value
340  * is not lower than the modulus, then the modulus is subtracted and
341  * this function returns 1; otherwise, it leaves it untouched and it
342  * returns 0.
343  */
344 static uint32_t
345 reduce_final_f255(uint32_t *d)
346 {
347 	uint32_t t[9];
348 	uint32_t cc;
349 	int i;
350 
351 	memcpy(t, d, sizeof t);
352 	cc = 19;
353 	for (i = 0; i < 9; i ++) {
354 		uint32_t w;
355 
356 		w = t[i] + cc;
357 		cc = w >> 30;
358 		t[i] = w & 0x3FFFFFFF;
359 	}
360 	cc = t[8] >> 15;
361 	t[8] &= 0x7FFF;
362 	CCOPY(cc, d, t, sizeof t);
363 	return cc;
364 }
365 
366 /*
367  * Perform a multiplication of two integers modulo 2^255-19.
368  * Operands are arrays of 9 words, each containing 30 bits of data, in
369  * little-endian order. Input value may be up to 2^256-1; on output, value
370  * fits on 256 bits and is lower than twice the modulus.
371  */
372 static void
373 f255_mul(uint32_t *d, const uint32_t *a, const uint32_t *b)
374 {
375 	uint32_t t[18], cc;
376 	int i;
377 
378 	/*
379 	 * Compute raw multiplication. All result words fit in 30 bits
380 	 * each; upper word (t[17]) must fit on 2 bits, since the product
381 	 * of two 256-bit integers must fit on 512 bits.
382 	 */
383 	mul9(t, a, b);
384 
385 	/*
386 	 * Modular reduction: each high word is added where necessary.
387 	 * Since the modulus is 2^255-19 and word 9 corresponds to
388 	 * offset 9*30 = 270, word 9+k must be added to word k with
389 	 * a factor of 19*2^15 = 622592. The extra bits in word 8 are also
390 	 * added that way.
391 	 *
392 	 * Keeping the carry on 32 bits helps with 32-bit architectures,
393 	 * and does not noticeably impact performance on 64-bit systems.
394 	 */
395 	cc = MUL15(t[8] >> 15, 19);  /* at most 19*(2^15-1) = 622573 */
396 	t[8] &= 0x7FFF;
397 	for (i = 0; i < 9; i ++) {
398 		uint64_t w;
399 
400 		w = (uint64_t)t[i] + (uint64_t)cc + MUL31(t[i + 9], 622592);
401 		t[i] = (uint32_t)w & 0x3FFFFFFF;
402 		cc = (uint32_t)(w >> 30);  /* at most 622592 */
403 	}
404 
405 	/*
406 	 * Original product was up to (2^256-1)^2, i.e. a 512-bit integer.
407 	 * This was split into two parts (upper of 257 bits, lower of 255
408 	 * bits), and the upper was added to the lower with a factor 19,
409 	 * which means that the intermediate value is less than 77*2^255
410 	 * (19*2^257 + 2^255). Therefore, the extra bits "t[8] >> 15" are
411 	 * less than 77, and the initial carry cc is at most 76*19 = 1444.
412 	 */
413 	cc = MUL15(t[8] >> 15, 19);
414 	t[8] &= 0x7FFF;
415 	for (i = 0; i < 9; i ++) {
416 		uint32_t z;
417 
418 		z = t[i] + cc;
419 		d[i] = z & 0x3FFFFFFF;
420 		cc = z >> 30;
421 	}
422 
423 	/*
424 	 * Final result is at most 2^255 + 1443. In particular, the last
425 	 * carry is necessarily 0, since t[8] was truncated to 15 bits.
426 	 */
427 }
428 
429 /*
430  * Perform a squaring of an integer modulo 2^255-19.
431  * Operands are arrays of 9 words, each containing 30 bits of data, in
432  * little-endian order. Input value may be up to 2^256-1; on output, value
433  * fits on 256 bits and is lower than twice the modulus.
434  */
435 static void
436 f255_square(uint32_t *d, const uint32_t *a)
437 {
438 	uint32_t t[18], cc;
439 	int i;
440 
441 	/*
442 	 * Compute raw squaring. All result words fit in 30 bits
443 	 * each; upper word (t[17]) must fit on 2 bits, since the square
444 	 * of a 256-bit integers must fit on 512 bits.
445 	 */
446 	square9(t, a);
447 
448 	/*
449 	 * Modular reduction: each high word is added where necessary.
450 	 * See f255_mul() for details on the reduction and carry limits.
451 	 */
452 	cc = MUL15(t[8] >> 15, 19);
453 	t[8] &= 0x7FFF;
454 	for (i = 0; i < 9; i ++) {
455 		uint64_t w;
456 
457 		w = (uint64_t)t[i] + (uint64_t)cc + MUL31(t[i + 9], 622592);
458 		t[i] = (uint32_t)w & 0x3FFFFFFF;
459 		cc = (uint32_t)(w >> 30);
460 	}
461 	cc = MUL15(t[8] >> 15, 19);
462 	t[8] &= 0x7FFF;
463 	for (i = 0; i < 9; i ++) {
464 		uint32_t z;
465 
466 		z = t[i] + cc;
467 		d[i] = z & 0x3FFFFFFF;
468 		cc = z >> 30;
469 	}
470 }
471 
472 /*
473  * Add two values in F255. Partial reduction is performed (down to less
474  * than twice the modulus).
475  */
476 static void
477 f255_add(uint32_t *d, const uint32_t *a, const uint32_t *b)
478 {
479 	/*
480 	 * Since operand words fit on 30 bits, we can use 32-bit
481 	 * variables throughout.
482 	 */
483 	int i;
484 	uint32_t cc, w;
485 
486 	cc = 0;
487 	for (i = 0; i < 9; i ++) {
488 		w = a[i] + b[i] + cc;
489 		d[i] = w & 0x3FFFFFFF;
490 		cc = w >> 30;
491 	}
492 	cc = MUL15(w >> 15, 19);
493 	d[8] &= 0x7FFF;
494 	for (i = 0; i < 9; i ++) {
495 		w = d[i] + cc;
496 		d[i] = w & 0x3FFFFFFF;
497 		cc = w >> 30;
498 	}
499 }
500 
501 /*
502  * Subtract one value from another in F255. Partial reduction is
503  * performed (down to less than twice the modulus).
504  */
505 static void
506 f255_sub(uint32_t *d, const uint32_t *a, const uint32_t *b)
507 {
508 	/*
509 	 * We actually compute a - b + 2*p, so that the final value is
510 	 * necessarily positive.
511 	 */
512 	int i;
513 	uint32_t cc, w;
514 
515 	cc = (uint32_t)-38;
516 	for (i = 0; i < 9; i ++) {
517 		w = a[i] - b[i] + cc;
518 		d[i] = w & 0x3FFFFFFF;
519 		cc = ARSH(w, 30);
520 	}
521 	cc = MUL15((w + 0x10000) >> 15, 19);
522 	d[8] &= 0x7FFF;
523 	for (i = 0; i < 9; i ++) {
524 		w = d[i] + cc;
525 		d[i] = w & 0x3FFFFFFF;
526 		cc = w >> 30;
527 	}
528 }
529 
530 /*
531  * Multiply an integer by the 'A24' constant (121665). Partial reduction
532  * is performed (down to less than twice the modulus).
533  */
534 static void
535 f255_mul_a24(uint32_t *d, const uint32_t *a)
536 {
537 	int i;
538 	uint64_t w;
539 	uint32_t cc;
540 
541 	/*
542 	 * a[] is over 256 bits, thus a[8] has length at most 16 bits.
543 	 * We single out the processing of the last word: intermediate
544 	 * value w is up to 121665*2^16, yielding a carry for the next
545 	 * loop of at most 19*(121665*2^16/2^15) = 4623289.
546 	 */
547 	cc = 0;
548 	for (i = 0; i < 8; i ++) {
549 		w = MUL31(a[i], 121665) + (uint64_t)cc;
550 		d[i] = (uint32_t)w & 0x3FFFFFFF;
551 		cc = (uint32_t)(w >> 30);
552 	}
553 	w = MUL31(a[8], 121665) + (uint64_t)cc;
554 	d[8] = (uint32_t)w & 0x7FFF;
555 	cc = MUL15((uint32_t)(w >> 15), 19);
556 
557 	for (i = 0; i < 9; i ++) {
558 		uint32_t z;
559 
560 		z = d[i] + cc;
561 		d[i] = z & 0x3FFFFFFF;
562 		cc = z >> 30;
563 	}
564 }
565 
566 static const unsigned char GEN[] = {
567 	0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
568 	0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
569 	0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
570 	0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
571 };
572 
573 static const unsigned char ORDER[] = {
574 	0x7F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
575 	0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
576 	0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
577 	0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
578 };
579 
580 static const unsigned char *
581 api_generator(int curve, size_t *len)
582 {
583 	(void)curve;
584 	*len = 32;
585 	return GEN;
586 }
587 
588 static const unsigned char *
589 api_order(int curve, size_t *len)
590 {
591 	(void)curve;
592 	*len = 32;
593 	return ORDER;
594 }
595 
596 static size_t
597 api_xoff(int curve, size_t *len)
598 {
599 	(void)curve;
600 	*len = 32;
601 	return 0;
602 }
603 
604 static void
605 cswap(uint32_t *a, uint32_t *b, uint32_t ctl)
606 {
607 	int i;
608 
609 	ctl = -ctl;
610 	for (i = 0; i < 9; i ++) {
611 		uint32_t aw, bw, tw;
612 
613 		aw = a[i];
614 		bw = b[i];
615 		tw = ctl & (aw ^ bw);
616 		a[i] = aw ^ tw;
617 		b[i] = bw ^ tw;
618 	}
619 }
620 
621 static uint32_t
622 api_mul(unsigned char *G, size_t Glen,
623 	const unsigned char *kb, size_t kblen, int curve)
624 {
625 	uint32_t x1[9], x2[9], x3[9], z2[9], z3[9];
626 	uint32_t a[9], aa[9], b[9], bb[9];
627 	uint32_t c[9], d[9], e[9], da[9], cb[9];
628 	unsigned char k[32];
629 	uint32_t swap;
630 	int i;
631 
632 	(void)curve;
633 
634 	/*
635 	 * Points are encoded over exactly 32 bytes. Multipliers must fit
636 	 * in 32 bytes as well.
637 	 * RFC 7748 mandates that the high bit of the last point byte must
638 	 * be ignored/cleared.
639 	 */
640 	if (Glen != 32 || kblen > 32) {
641 		return 0;
642 	}
643 	G[31] &= 0x7F;
644 
645 	/*
646 	 * Initialise variables x1, x2, z2, x3 and z3. We set all of them
647 	 * into Montgomery representation.
648 	 */
649 	x1[8] = le8_to_le30(x1, G, 32);
650 	memcpy(x3, x1, sizeof x1);
651 	memset(z2, 0, sizeof z2);
652 	memset(x2, 0, sizeof x2);
653 	x2[0] = 1;
654 	memset(z3, 0, sizeof z3);
655 	z3[0] = 1;
656 
657 	memset(k, 0, (sizeof k) - kblen);
658 	memcpy(k + (sizeof k) - kblen, kb, kblen);
659 	k[31] &= 0xF8;
660 	k[0] &= 0x7F;
661 	k[0] |= 0x40;
662 
663 	/* obsolete
664 	print_int("x1", x1);
665 	*/
666 
667 	swap = 0;
668 	for (i = 254; i >= 0; i --) {
669 		uint32_t kt;
670 
671 		kt = (k[31 - (i >> 3)] >> (i & 7)) & 1;
672 		swap ^= kt;
673 		cswap(x2, x3, swap);
674 		cswap(z2, z3, swap);
675 		swap = kt;
676 
677 		/* obsolete
678 		print_int("x2", x2);
679 		print_int("z2", z2);
680 		print_int("x3", x3);
681 		print_int("z3", z3);
682 		*/
683 
684 		f255_add(a, x2, z2);
685 		f255_square(aa, a);
686 		f255_sub(b, x2, z2);
687 		f255_square(bb, b);
688 		f255_sub(e, aa, bb);
689 		f255_add(c, x3, z3);
690 		f255_sub(d, x3, z3);
691 		f255_mul(da, d, a);
692 		f255_mul(cb, c, b);
693 
694 		/* obsolete
695 		print_int("a ", a);
696 		print_int("aa", aa);
697 		print_int("b ", b);
698 		print_int("bb", bb);
699 		print_int("e ", e);
700 		print_int("c ", c);
701 		print_int("d ", d);
702 		print_int("da", da);
703 		print_int("cb", cb);
704 		*/
705 
706 		f255_add(x3, da, cb);
707 		f255_square(x3, x3);
708 		f255_sub(z3, da, cb);
709 		f255_square(z3, z3);
710 		f255_mul(z3, z3, x1);
711 		f255_mul(x2, aa, bb);
712 		f255_mul_a24(z2, e);
713 		f255_add(z2, z2, aa);
714 		f255_mul(z2, e, z2);
715 
716 		/* obsolete
717 		print_int("x2", x2);
718 		print_int("z2", z2);
719 		print_int("x3", x3);
720 		print_int("z3", z3);
721 		*/
722 	}
723 	cswap(x2, x3, swap);
724 	cswap(z2, z3, swap);
725 
726 	/*
727 	 * Inverse z2 with a modular exponentiation. This is a simple
728 	 * square-and-multiply algorithm; we mutualise most non-squarings
729 	 * since the exponent contains almost only ones.
730 	 */
731 	memcpy(a, z2, sizeof z2);
732 	for (i = 0; i < 15; i ++) {
733 		f255_square(a, a);
734 		f255_mul(a, a, z2);
735 	}
736 	memcpy(b, a, sizeof a);
737 	for (i = 0; i < 14; i ++) {
738 		int j;
739 
740 		for (j = 0; j < 16; j ++) {
741 			f255_square(b, b);
742 		}
743 		f255_mul(b, b, a);
744 	}
745 	for (i = 14; i >= 0; i --) {
746 		f255_square(b, b);
747 		if ((0xFFEB >> i) & 1) {
748 			f255_mul(b, z2, b);
749 		}
750 	}
751 	f255_mul(x2, x2, b);
752 	reduce_final_f255(x2);
753 	le30_to_le8(G, 32, x2);
754 	return 1;
755 }
756 
757 static size_t
758 api_mulgen(unsigned char *R,
759 	const unsigned char *x, size_t xlen, int curve)
760 {
761 	const unsigned char *G;
762 	size_t Glen;
763 
764 	G = api_generator(curve, &Glen);
765 	memcpy(R, G, Glen);
766 	api_mul(R, Glen, x, xlen, curve);
767 	return Glen;
768 }
769 
770 static uint32_t
771 api_muladd(unsigned char *A, const unsigned char *B, size_t len,
772 	const unsigned char *x, size_t xlen,
773 	const unsigned char *y, size_t ylen, int curve)
774 {
775 	/*
776 	 * We don't implement this method, since it is used for ECDSA
777 	 * only, and there is no ECDSA over Curve25519 (which instead
778 	 * uses EdDSA).
779 	 */
780 	(void)A;
781 	(void)B;
782 	(void)len;
783 	(void)x;
784 	(void)xlen;
785 	(void)y;
786 	(void)ylen;
787 	(void)curve;
788 	return 0;
789 }
790 
791 /* see bearssl_ec.h */
792 const br_ec_impl br_ec_c25519_m31 = {
793 	(uint32_t)0x20000000,
794 	&api_generator,
795 	&api_order,
796 	&api_xoff,
797 	&api_mul,
798 	&api_mulgen,
799 	&api_muladd
800 };
801