2 * Copyright (c) 2017 Thomas Pornin <pornin@bolet.org>
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:
12 * The above copyright notice and this permission notice shall be
13 * included in all copies or substantial portions of the Software.
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
31 print_int(const char *name, const uint32_t *x)
34 unsigned char tmp[36];
36 printf("%s = ", name);
37 for (u = 0; u < 20; u ++) {
40 for (u = 0; u < 20; u ++) {
41 printf(" %04X", x[u]);
47 memset(tmp, 0, sizeof tmp);
48 for (u = 0; u < 20; u ++) {
60 tmp[35 - k] |= (unsigned char)w;
61 tmp[34 - k] |= (unsigned char)(w >> 8);
62 tmp[33 - k] |= (unsigned char)(w >> 16);
63 tmp[32 - k] |= (unsigned char)(w >> 24);
65 for (u = 4; u < 36; u ++) {
66 printf("%02X", tmp[u]);
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.
83 #define ARSH(x, n) (((uint32_t)(x) >> (n)) \
84 | ((-((uint32_t)(x) >> 31)) << (32 - (n))))
86 #define ARSH(x, n) ((*(int32_t *)&(x)) >> (n))
90 * Convert an integer from unsigned little-endian encoding to a sequence of
91 * 13-bit words in little-endian order. The final "partial" word is
95 le8_to_le13(uint32_t *dst
, const unsigned char *src
, size_t len
)
103 acc
|= (uint32_t)(*src
++) << acc_len
;
106 *dst
++ = acc
& 0x1FFF;
115 * Convert an integer (13-bit words, little-endian) to unsigned
116 * little-endian encoding. The total encoding length is provided; all
117 * the destination bytes will be filled.
120 le13_to_le8(unsigned char *dst
, size_t len
, const uint32_t *src
)
129 acc
|= (*src
++) << acc_len
;
132 *dst
++ = (unsigned char)acc
;
139 * Normalise an array of words to a strict 13 bits per word. Returned
140 * value is the resulting carry. The source (w) and destination (d)
141 * arrays may be identical, but shall not overlap partially.
143 static inline uint32_t
144 norm13(uint32_t *d
, const uint32_t *w
, size_t len
)
150 for (u
= 0; u
< len
; u
++) {
161 * mul20() multiplies two 260-bit integers together. Each word must fit
162 * on 13 bits; source operands use 20 words, destination operand
163 * receives 40 words. All overlaps allowed.
165 * square20() computes the square of a 260-bit integer. Each word must
166 * fit on 13 bits; source operand uses 20 words, destination operand
167 * receives 40 words. All overlaps allowed.
173 mul20(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
176 * Two-level Karatsuba: turns a 20x20 multiplication into
177 * nine 5x5 multiplications. We use 13-bit words but do not
178 * propagate carries immediately, so words may expand:
180 * - First Karatsuba decomposition turns the 20x20 mul on
181 * 13-bit words into three 10x10 muls, two on 13-bit words
182 * and one on 14-bit words.
184 * - Second Karatsuba decomposition further splits these into:
186 * * four 5x5 muls on 13-bit words
187 * * four 5x5 muls on 14-bit words
188 * * one 5x5 mul on 15-bit words
190 * Highest word value is 8191, 16382 or 32764, for 13-bit, 14-bit
191 * or 15-bit words, respectively.
193 uint32_t u
[45], v
[45], w
[90];
197 #define ZADD(dw, d_off, s1w, s1_off, s2w, s2_off) do { \
198 (dw)[5 * (d_off) + 0] = (s1w)[5 * (s1_off) + 0] \
199 + (s2w)[5 * (s2_off) + 0]; \
200 (dw)[5 * (d_off) + 1] = (s1w)[5 * (s1_off) + 1] \
201 + (s2w)[5 * (s2_off) + 1]; \
202 (dw)[5 * (d_off) + 2] = (s1w)[5 * (s1_off) + 2] \
203 + (s2w)[5 * (s2_off) + 2]; \
204 (dw)[5 * (d_off) + 3] = (s1w)[5 * (s1_off) + 3] \
205 + (s2w)[5 * (s2_off) + 3]; \
206 (dw)[5 * (d_off) + 4] = (s1w)[5 * (s1_off) + 4] \
207 + (s2w)[5 * (s2_off) + 4]; \
210 #define ZADDT(dw, d_off, sw, s_off) do { \
211 (dw)[5 * (d_off) + 0] += (sw)[5 * (s_off) + 0]; \
212 (dw)[5 * (d_off) + 1] += (sw)[5 * (s_off) + 1]; \
213 (dw)[5 * (d_off) + 2] += (sw)[5 * (s_off) + 2]; \
214 (dw)[5 * (d_off) + 3] += (sw)[5 * (s_off) + 3]; \
215 (dw)[5 * (d_off) + 4] += (sw)[5 * (s_off) + 4]; \
218 #define ZSUB2F(dw, d_off, s1w, s1_off, s2w, s2_off) do { \
219 (dw)[5 * (d_off) + 0] -= (s1w)[5 * (s1_off) + 0] \
220 + (s2w)[5 * (s2_off) + 0]; \
221 (dw)[5 * (d_off) + 1] -= (s1w)[5 * (s1_off) + 1] \
222 + (s2w)[5 * (s2_off) + 1]; \
223 (dw)[5 * (d_off) + 2] -= (s1w)[5 * (s1_off) + 2] \
224 + (s2w)[5 * (s2_off) + 2]; \
225 (dw)[5 * (d_off) + 3] -= (s1w)[5 * (s1_off) + 3] \
226 + (s2w)[5 * (s2_off) + 3]; \
227 (dw)[5 * (d_off) + 4] -= (s1w)[5 * (s1_off) + 4] \
228 + (s2w)[5 * (s2_off) + 4]; \
231 #define CPR1(w, cprcc) do { \
232 uint32_t cprz = (w) + cprcc; \
233 (w) = cprz & 0x1FFF; \
234 cprcc = cprz >> 13; \
237 #define CPR(dw, d_off) do { \
240 CPR1((dw)[(d_off) + 0], cprcc); \
241 CPR1((dw)[(d_off) + 1], cprcc); \
242 CPR1((dw)[(d_off) + 2], cprcc); \
243 CPR1((dw)[(d_off) + 3], cprcc); \
244 CPR1((dw)[(d_off) + 4], cprcc); \
245 CPR1((dw)[(d_off) + 5], cprcc); \
246 CPR1((dw)[(d_off) + 6], cprcc); \
247 CPR1((dw)[(d_off) + 7], cprcc); \
248 CPR1((dw)[(d_off) + 8], cprcc); \
249 (dw)[(d_off) + 9] = cprcc; \
252 memcpy(u
, a
, 20 * sizeof *a
);
253 ZADD(u
, 4, a
, 0, a
, 1);
254 ZADD(u
, 5, a
, 2, a
, 3);
255 ZADD(u
, 6, a
, 0, a
, 2);
256 ZADD(u
, 7, a
, 1, a
, 3);
257 ZADD(u
, 8, u
, 6, u
, 7);
259 memcpy(v
, b
, 20 * sizeof *b
);
260 ZADD(v
, 4, b
, 0, b
, 1);
261 ZADD(v
, 5, b
, 2, b
, 3);
262 ZADD(v
, 6, b
, 0, b
, 2);
263 ZADD(v
, 7, b
, 1, b
, 3);
264 ZADD(v
, 8, v
, 6, v
, 7);
267 * Do the eight first 8x8 muls. Source words are at most 16382
268 * each, so we can add product results together "as is" in 32-bit
271 for (i
= 0; i
< 40; i
+= 5) {
272 w
[(i
<< 1) + 0] = MUL15(u
[i
+ 0], v
[i
+ 0]);
273 w
[(i
<< 1) + 1] = MUL15(u
[i
+ 0], v
[i
+ 1])
274 + MUL15(u
[i
+ 1], v
[i
+ 0]);
275 w
[(i
<< 1) + 2] = MUL15(u
[i
+ 0], v
[i
+ 2])
276 + MUL15(u
[i
+ 1], v
[i
+ 1])
277 + MUL15(u
[i
+ 2], v
[i
+ 0]);
278 w
[(i
<< 1) + 3] = MUL15(u
[i
+ 0], v
[i
+ 3])
279 + MUL15(u
[i
+ 1], v
[i
+ 2])
280 + MUL15(u
[i
+ 2], v
[i
+ 1])
281 + MUL15(u
[i
+ 3], v
[i
+ 0]);
282 w
[(i
<< 1) + 4] = MUL15(u
[i
+ 0], v
[i
+ 4])
283 + MUL15(u
[i
+ 1], v
[i
+ 3])
284 + MUL15(u
[i
+ 2], v
[i
+ 2])
285 + MUL15(u
[i
+ 3], v
[i
+ 1])
286 + MUL15(u
[i
+ 4], v
[i
+ 0]);
287 w
[(i
<< 1) + 5] = MUL15(u
[i
+ 1], v
[i
+ 4])
288 + MUL15(u
[i
+ 2], v
[i
+ 3])
289 + MUL15(u
[i
+ 3], v
[i
+ 2])
290 + MUL15(u
[i
+ 4], v
[i
+ 1]);
291 w
[(i
<< 1) + 6] = MUL15(u
[i
+ 2], v
[i
+ 4])
292 + MUL15(u
[i
+ 3], v
[i
+ 3])
293 + MUL15(u
[i
+ 4], v
[i
+ 2]);
294 w
[(i
<< 1) + 7] = MUL15(u
[i
+ 3], v
[i
+ 4])
295 + MUL15(u
[i
+ 4], v
[i
+ 3]);
296 w
[(i
<< 1) + 8] = MUL15(u
[i
+ 4], v
[i
+ 4]);
301 * For the 9th multiplication, source words are up to 32764,
302 * so we must do some carry propagation. If we add up to
303 * 4 products and the carry is no more than 524224, then the
304 * result fits in 32 bits, and the next carry will be no more
305 * than 524224 (because 4*(32764^2)+524224 < 8192*524225).
307 * We thus just skip one of the products in the middle word,
308 * then do a carry propagation (this reduces words to 13 bits
309 * each, except possibly the last, which may use up to 17 bits
310 * or so), then add the missing product.
312 w
[80 + 0] = MUL15(u
[40 + 0], v
[40 + 0]);
313 w
[80 + 1] = MUL15(u
[40 + 0], v
[40 + 1])
314 + MUL15(u
[40 + 1], v
[40 + 0]);
315 w
[80 + 2] = MUL15(u
[40 + 0], v
[40 + 2])
316 + MUL15(u
[40 + 1], v
[40 + 1])
317 + MUL15(u
[40 + 2], v
[40 + 0]);
318 w
[80 + 3] = MUL15(u
[40 + 0], v
[40 + 3])
319 + MUL15(u
[40 + 1], v
[40 + 2])
320 + MUL15(u
[40 + 2], v
[40 + 1])
321 + MUL15(u
[40 + 3], v
[40 + 0]);
322 w
[80 + 4] = MUL15(u
[40 + 0], v
[40 + 4])
323 + MUL15(u
[40 + 1], v
[40 + 3])
324 + MUL15(u
[40 + 2], v
[40 + 2])
325 + MUL15(u
[40 + 3], v
[40 + 1]);
326 /* + MUL15(u[40 + 4], v[40 + 0]) */
327 w
[80 + 5] = MUL15(u
[40 + 1], v
[40 + 4])
328 + MUL15(u
[40 + 2], v
[40 + 3])
329 + MUL15(u
[40 + 3], v
[40 + 2])
330 + MUL15(u
[40 + 4], v
[40 + 1]);
331 w
[80 + 6] = MUL15(u
[40 + 2], v
[40 + 4])
332 + MUL15(u
[40 + 3], v
[40 + 3])
333 + MUL15(u
[40 + 4], v
[40 + 2]);
334 w
[80 + 7] = MUL15(u
[40 + 3], v
[40 + 4])
335 + MUL15(u
[40 + 4], v
[40 + 3]);
336 w
[80 + 8] = MUL15(u
[40 + 4], v
[40 + 4]);
340 w
[80 + 4] += MUL15(u
[40 + 4], v
[40 + 0]);
343 * The products on 14-bit words in slots 6 and 7 yield values
344 * up to 5*(16382^2) each, and we need to subtract two such
345 * values from the higher word. We need the subtraction to fit
346 * in a _signed_ 32-bit integer, i.e. 31 bits + a sign bit.
347 * However, 10*(16382^2) does not fit. So we must perform a
348 * bit of reduction here.
357 /* 0..1*0..1 into 0..3 */
358 ZSUB2F(w
, 8, w
, 0, w
, 2);
359 ZSUB2F(w
, 9, w
, 1, w
, 3);
363 /* 2..3*2..3 into 4..7 */
364 ZSUB2F(w
, 10, w
, 4, w
, 6);
365 ZSUB2F(w
, 11, w
, 5, w
, 7);
369 /* (0..1+2..3)*(0..1+2..3) into 12..15 */
370 ZSUB2F(w
, 16, w
, 12, w
, 14);
371 ZSUB2F(w
, 17, w
, 13, w
, 15);
375 /* first-level recomposition */
376 ZSUB2F(w
, 12, w
, 0, w
, 4);
377 ZSUB2F(w
, 13, w
, 1, w
, 5);
378 ZSUB2F(w
, 14, w
, 2, w
, 6);
379 ZSUB2F(w
, 15, w
, 3, w
, 7);
386 * Perform carry propagation to bring all words down to 13 bits.
388 cc
= norm13(d
, w
, 40);
399 square20(uint32_t *d
, const uint32_t *a
)
407 mul20(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
411 t
[ 0] = MUL15(a
[ 0], b
[ 0]);
412 t
[ 1] = MUL15(a
[ 0], b
[ 1])
413 + MUL15(a
[ 1], b
[ 0]);
414 t
[ 2] = MUL15(a
[ 0], b
[ 2])
415 + MUL15(a
[ 1], b
[ 1])
416 + MUL15(a
[ 2], b
[ 0]);
417 t
[ 3] = MUL15(a
[ 0], b
[ 3])
418 + MUL15(a
[ 1], b
[ 2])
419 + MUL15(a
[ 2], b
[ 1])
420 + MUL15(a
[ 3], b
[ 0]);
421 t
[ 4] = MUL15(a
[ 0], b
[ 4])
422 + MUL15(a
[ 1], b
[ 3])
423 + MUL15(a
[ 2], b
[ 2])
424 + MUL15(a
[ 3], b
[ 1])
425 + MUL15(a
[ 4], b
[ 0]);
426 t
[ 5] = MUL15(a
[ 0], b
[ 5])
427 + MUL15(a
[ 1], b
[ 4])
428 + MUL15(a
[ 2], b
[ 3])
429 + MUL15(a
[ 3], b
[ 2])
430 + MUL15(a
[ 4], b
[ 1])
431 + MUL15(a
[ 5], b
[ 0]);
432 t
[ 6] = MUL15(a
[ 0], b
[ 6])
433 + MUL15(a
[ 1], b
[ 5])
434 + MUL15(a
[ 2], b
[ 4])
435 + MUL15(a
[ 3], b
[ 3])
436 + MUL15(a
[ 4], b
[ 2])
437 + MUL15(a
[ 5], b
[ 1])
438 + MUL15(a
[ 6], b
[ 0]);
439 t
[ 7] = MUL15(a
[ 0], b
[ 7])
440 + MUL15(a
[ 1], b
[ 6])
441 + MUL15(a
[ 2], b
[ 5])
442 + MUL15(a
[ 3], b
[ 4])
443 + MUL15(a
[ 4], b
[ 3])
444 + MUL15(a
[ 5], b
[ 2])
445 + MUL15(a
[ 6], b
[ 1])
446 + MUL15(a
[ 7], b
[ 0]);
447 t
[ 8] = MUL15(a
[ 0], b
[ 8])
448 + MUL15(a
[ 1], b
[ 7])
449 + MUL15(a
[ 2], b
[ 6])
450 + MUL15(a
[ 3], b
[ 5])
451 + MUL15(a
[ 4], b
[ 4])
452 + MUL15(a
[ 5], b
[ 3])
453 + MUL15(a
[ 6], b
[ 2])
454 + MUL15(a
[ 7], b
[ 1])
455 + MUL15(a
[ 8], b
[ 0]);
456 t
[ 9] = MUL15(a
[ 0], b
[ 9])
457 + MUL15(a
[ 1], b
[ 8])
458 + MUL15(a
[ 2], b
[ 7])
459 + MUL15(a
[ 3], b
[ 6])
460 + MUL15(a
[ 4], b
[ 5])
461 + MUL15(a
[ 5], b
[ 4])
462 + MUL15(a
[ 6], b
[ 3])
463 + MUL15(a
[ 7], b
[ 2])
464 + MUL15(a
[ 8], b
[ 1])
465 + MUL15(a
[ 9], b
[ 0]);
466 t
[10] = MUL15(a
[ 0], b
[10])
467 + MUL15(a
[ 1], b
[ 9])
468 + MUL15(a
[ 2], b
[ 8])
469 + MUL15(a
[ 3], b
[ 7])
470 + MUL15(a
[ 4], b
[ 6])
471 + MUL15(a
[ 5], b
[ 5])
472 + MUL15(a
[ 6], b
[ 4])
473 + MUL15(a
[ 7], b
[ 3])
474 + MUL15(a
[ 8], b
[ 2])
475 + MUL15(a
[ 9], b
[ 1])
476 + MUL15(a
[10], b
[ 0]);
477 t
[11] = MUL15(a
[ 0], b
[11])
478 + MUL15(a
[ 1], b
[10])
479 + MUL15(a
[ 2], b
[ 9])
480 + MUL15(a
[ 3], b
[ 8])
481 + MUL15(a
[ 4], b
[ 7])
482 + MUL15(a
[ 5], b
[ 6])
483 + MUL15(a
[ 6], b
[ 5])
484 + MUL15(a
[ 7], b
[ 4])
485 + MUL15(a
[ 8], b
[ 3])
486 + MUL15(a
[ 9], b
[ 2])
487 + MUL15(a
[10], b
[ 1])
488 + MUL15(a
[11], b
[ 0]);
489 t
[12] = MUL15(a
[ 0], b
[12])
490 + MUL15(a
[ 1], b
[11])
491 + MUL15(a
[ 2], b
[10])
492 + MUL15(a
[ 3], b
[ 9])
493 + MUL15(a
[ 4], b
[ 8])
494 + MUL15(a
[ 5], b
[ 7])
495 + MUL15(a
[ 6], b
[ 6])
496 + MUL15(a
[ 7], b
[ 5])
497 + MUL15(a
[ 8], b
[ 4])
498 + MUL15(a
[ 9], b
[ 3])
499 + MUL15(a
[10], b
[ 2])
500 + MUL15(a
[11], b
[ 1])
501 + MUL15(a
[12], b
[ 0]);
502 t
[13] = MUL15(a
[ 0], b
[13])
503 + MUL15(a
[ 1], b
[12])
504 + MUL15(a
[ 2], b
[11])
505 + MUL15(a
[ 3], b
[10])
506 + MUL15(a
[ 4], b
[ 9])
507 + MUL15(a
[ 5], b
[ 8])
508 + MUL15(a
[ 6], b
[ 7])
509 + MUL15(a
[ 7], b
[ 6])
510 + MUL15(a
[ 8], b
[ 5])
511 + MUL15(a
[ 9], b
[ 4])
512 + MUL15(a
[10], b
[ 3])
513 + MUL15(a
[11], b
[ 2])
514 + MUL15(a
[12], b
[ 1])
515 + MUL15(a
[13], b
[ 0]);
516 t
[14] = MUL15(a
[ 0], b
[14])
517 + MUL15(a
[ 1], b
[13])
518 + MUL15(a
[ 2], b
[12])
519 + MUL15(a
[ 3], b
[11])
520 + MUL15(a
[ 4], b
[10])
521 + MUL15(a
[ 5], b
[ 9])
522 + MUL15(a
[ 6], b
[ 8])
523 + MUL15(a
[ 7], b
[ 7])
524 + MUL15(a
[ 8], b
[ 6])
525 + MUL15(a
[ 9], b
[ 5])
526 + MUL15(a
[10], b
[ 4])
527 + MUL15(a
[11], b
[ 3])
528 + MUL15(a
[12], b
[ 2])
529 + MUL15(a
[13], b
[ 1])
530 + MUL15(a
[14], b
[ 0]);
531 t
[15] = MUL15(a
[ 0], b
[15])
532 + MUL15(a
[ 1], b
[14])
533 + MUL15(a
[ 2], b
[13])
534 + MUL15(a
[ 3], b
[12])
535 + MUL15(a
[ 4], b
[11])
536 + MUL15(a
[ 5], b
[10])
537 + MUL15(a
[ 6], b
[ 9])
538 + MUL15(a
[ 7], b
[ 8])
539 + MUL15(a
[ 8], b
[ 7])
540 + MUL15(a
[ 9], b
[ 6])
541 + MUL15(a
[10], b
[ 5])
542 + MUL15(a
[11], b
[ 4])
543 + MUL15(a
[12], b
[ 3])
544 + MUL15(a
[13], b
[ 2])
545 + MUL15(a
[14], b
[ 1])
546 + MUL15(a
[15], b
[ 0]);
547 t
[16] = MUL15(a
[ 0], b
[16])
548 + MUL15(a
[ 1], b
[15])
549 + MUL15(a
[ 2], b
[14])
550 + MUL15(a
[ 3], b
[13])
551 + MUL15(a
[ 4], b
[12])
552 + MUL15(a
[ 5], b
[11])
553 + MUL15(a
[ 6], b
[10])
554 + MUL15(a
[ 7], b
[ 9])
555 + MUL15(a
[ 8], b
[ 8])
556 + MUL15(a
[ 9], b
[ 7])
557 + MUL15(a
[10], b
[ 6])
558 + MUL15(a
[11], b
[ 5])
559 + MUL15(a
[12], b
[ 4])
560 + MUL15(a
[13], b
[ 3])
561 + MUL15(a
[14], b
[ 2])
562 + MUL15(a
[15], b
[ 1])
563 + MUL15(a
[16], b
[ 0]);
564 t
[17] = MUL15(a
[ 0], b
[17])
565 + MUL15(a
[ 1], b
[16])
566 + MUL15(a
[ 2], b
[15])
567 + MUL15(a
[ 3], b
[14])
568 + MUL15(a
[ 4], b
[13])
569 + MUL15(a
[ 5], b
[12])
570 + MUL15(a
[ 6], b
[11])
571 + MUL15(a
[ 7], b
[10])
572 + MUL15(a
[ 8], b
[ 9])
573 + MUL15(a
[ 9], b
[ 8])
574 + MUL15(a
[10], b
[ 7])
575 + MUL15(a
[11], b
[ 6])
576 + MUL15(a
[12], b
[ 5])
577 + MUL15(a
[13], b
[ 4])
578 + MUL15(a
[14], b
[ 3])
579 + MUL15(a
[15], b
[ 2])
580 + MUL15(a
[16], b
[ 1])
581 + MUL15(a
[17], b
[ 0]);
582 t
[18] = MUL15(a
[ 0], b
[18])
583 + MUL15(a
[ 1], b
[17])
584 + MUL15(a
[ 2], b
[16])
585 + MUL15(a
[ 3], b
[15])
586 + MUL15(a
[ 4], b
[14])
587 + MUL15(a
[ 5], b
[13])
588 + MUL15(a
[ 6], b
[12])
589 + MUL15(a
[ 7], b
[11])
590 + MUL15(a
[ 8], b
[10])
591 + MUL15(a
[ 9], b
[ 9])
592 + MUL15(a
[10], b
[ 8])
593 + MUL15(a
[11], b
[ 7])
594 + MUL15(a
[12], b
[ 6])
595 + MUL15(a
[13], b
[ 5])
596 + MUL15(a
[14], b
[ 4])
597 + MUL15(a
[15], b
[ 3])
598 + MUL15(a
[16], b
[ 2])
599 + MUL15(a
[17], b
[ 1])
600 + MUL15(a
[18], b
[ 0]);
601 t
[19] = MUL15(a
[ 0], b
[19])
602 + MUL15(a
[ 1], b
[18])
603 + MUL15(a
[ 2], b
[17])
604 + MUL15(a
[ 3], b
[16])
605 + MUL15(a
[ 4], b
[15])
606 + MUL15(a
[ 5], b
[14])
607 + MUL15(a
[ 6], b
[13])
608 + MUL15(a
[ 7], b
[12])
609 + MUL15(a
[ 8], b
[11])
610 + MUL15(a
[ 9], b
[10])
611 + MUL15(a
[10], b
[ 9])
612 + MUL15(a
[11], b
[ 8])
613 + MUL15(a
[12], b
[ 7])
614 + MUL15(a
[13], b
[ 6])
615 + MUL15(a
[14], b
[ 5])
616 + MUL15(a
[15], b
[ 4])
617 + MUL15(a
[16], b
[ 3])
618 + MUL15(a
[17], b
[ 2])
619 + MUL15(a
[18], b
[ 1])
620 + MUL15(a
[19], b
[ 0]);
621 t
[20] = MUL15(a
[ 1], b
[19])
622 + MUL15(a
[ 2], b
[18])
623 + MUL15(a
[ 3], b
[17])
624 + MUL15(a
[ 4], b
[16])
625 + MUL15(a
[ 5], b
[15])
626 + MUL15(a
[ 6], b
[14])
627 + MUL15(a
[ 7], b
[13])
628 + MUL15(a
[ 8], b
[12])
629 + MUL15(a
[ 9], b
[11])
630 + MUL15(a
[10], b
[10])
631 + MUL15(a
[11], b
[ 9])
632 + MUL15(a
[12], b
[ 8])
633 + MUL15(a
[13], b
[ 7])
634 + MUL15(a
[14], b
[ 6])
635 + MUL15(a
[15], b
[ 5])
636 + MUL15(a
[16], b
[ 4])
637 + MUL15(a
[17], b
[ 3])
638 + MUL15(a
[18], b
[ 2])
639 + MUL15(a
[19], b
[ 1]);
640 t
[21] = MUL15(a
[ 2], b
[19])
641 + MUL15(a
[ 3], b
[18])
642 + MUL15(a
[ 4], b
[17])
643 + MUL15(a
[ 5], b
[16])
644 + MUL15(a
[ 6], b
[15])
645 + MUL15(a
[ 7], b
[14])
646 + MUL15(a
[ 8], b
[13])
647 + MUL15(a
[ 9], b
[12])
648 + MUL15(a
[10], b
[11])
649 + MUL15(a
[11], b
[10])
650 + MUL15(a
[12], b
[ 9])
651 + MUL15(a
[13], b
[ 8])
652 + MUL15(a
[14], b
[ 7])
653 + MUL15(a
[15], b
[ 6])
654 + MUL15(a
[16], b
[ 5])
655 + MUL15(a
[17], b
[ 4])
656 + MUL15(a
[18], b
[ 3])
657 + MUL15(a
[19], b
[ 2]);
658 t
[22] = MUL15(a
[ 3], b
[19])
659 + MUL15(a
[ 4], b
[18])
660 + MUL15(a
[ 5], b
[17])
661 + MUL15(a
[ 6], b
[16])
662 + MUL15(a
[ 7], b
[15])
663 + MUL15(a
[ 8], b
[14])
664 + MUL15(a
[ 9], b
[13])
665 + MUL15(a
[10], b
[12])
666 + MUL15(a
[11], b
[11])
667 + MUL15(a
[12], b
[10])
668 + MUL15(a
[13], b
[ 9])
669 + MUL15(a
[14], b
[ 8])
670 + MUL15(a
[15], b
[ 7])
671 + MUL15(a
[16], b
[ 6])
672 + MUL15(a
[17], b
[ 5])
673 + MUL15(a
[18], b
[ 4])
674 + MUL15(a
[19], b
[ 3]);
675 t
[23] = MUL15(a
[ 4], b
[19])
676 + MUL15(a
[ 5], b
[18])
677 + MUL15(a
[ 6], b
[17])
678 + MUL15(a
[ 7], b
[16])
679 + MUL15(a
[ 8], b
[15])
680 + MUL15(a
[ 9], b
[14])
681 + MUL15(a
[10], b
[13])
682 + MUL15(a
[11], b
[12])
683 + MUL15(a
[12], b
[11])
684 + MUL15(a
[13], b
[10])
685 + MUL15(a
[14], b
[ 9])
686 + MUL15(a
[15], b
[ 8])
687 + MUL15(a
[16], b
[ 7])
688 + MUL15(a
[17], b
[ 6])
689 + MUL15(a
[18], b
[ 5])
690 + MUL15(a
[19], b
[ 4]);
691 t
[24] = MUL15(a
[ 5], b
[19])
692 + MUL15(a
[ 6], b
[18])
693 + MUL15(a
[ 7], b
[17])
694 + MUL15(a
[ 8], b
[16])
695 + MUL15(a
[ 9], b
[15])
696 + MUL15(a
[10], b
[14])
697 + MUL15(a
[11], b
[13])
698 + MUL15(a
[12], b
[12])
699 + MUL15(a
[13], b
[11])
700 + MUL15(a
[14], b
[10])
701 + MUL15(a
[15], b
[ 9])
702 + MUL15(a
[16], b
[ 8])
703 + MUL15(a
[17], b
[ 7])
704 + MUL15(a
[18], b
[ 6])
705 + MUL15(a
[19], b
[ 5]);
706 t
[25] = MUL15(a
[ 6], b
[19])
707 + MUL15(a
[ 7], b
[18])
708 + MUL15(a
[ 8], b
[17])
709 + MUL15(a
[ 9], b
[16])
710 + MUL15(a
[10], b
[15])
711 + MUL15(a
[11], b
[14])
712 + MUL15(a
[12], b
[13])
713 + MUL15(a
[13], b
[12])
714 + MUL15(a
[14], b
[11])
715 + MUL15(a
[15], b
[10])
716 + MUL15(a
[16], b
[ 9])
717 + MUL15(a
[17], b
[ 8])
718 + MUL15(a
[18], b
[ 7])
719 + MUL15(a
[19], b
[ 6]);
720 t
[26] = MUL15(a
[ 7], b
[19])
721 + MUL15(a
[ 8], b
[18])
722 + MUL15(a
[ 9], b
[17])
723 + MUL15(a
[10], b
[16])
724 + MUL15(a
[11], b
[15])
725 + MUL15(a
[12], b
[14])
726 + MUL15(a
[13], b
[13])
727 + MUL15(a
[14], b
[12])
728 + MUL15(a
[15], b
[11])
729 + MUL15(a
[16], b
[10])
730 + MUL15(a
[17], b
[ 9])
731 + MUL15(a
[18], b
[ 8])
732 + MUL15(a
[19], b
[ 7]);
733 t
[27] = MUL15(a
[ 8], b
[19])
734 + MUL15(a
[ 9], b
[18])
735 + MUL15(a
[10], b
[17])
736 + MUL15(a
[11], b
[16])
737 + MUL15(a
[12], b
[15])
738 + MUL15(a
[13], b
[14])
739 + MUL15(a
[14], b
[13])
740 + MUL15(a
[15], b
[12])
741 + MUL15(a
[16], b
[11])
742 + MUL15(a
[17], b
[10])
743 + MUL15(a
[18], b
[ 9])
744 + MUL15(a
[19], b
[ 8]);
745 t
[28] = MUL15(a
[ 9], b
[19])
746 + MUL15(a
[10], b
[18])
747 + MUL15(a
[11], b
[17])
748 + MUL15(a
[12], b
[16])
749 + MUL15(a
[13], b
[15])
750 + MUL15(a
[14], b
[14])
751 + MUL15(a
[15], b
[13])
752 + MUL15(a
[16], b
[12])
753 + MUL15(a
[17], b
[11])
754 + MUL15(a
[18], b
[10])
755 + MUL15(a
[19], b
[ 9]);
756 t
[29] = MUL15(a
[10], b
[19])
757 + MUL15(a
[11], b
[18])
758 + MUL15(a
[12], b
[17])
759 + MUL15(a
[13], b
[16])
760 + MUL15(a
[14], b
[15])
761 + MUL15(a
[15], b
[14])
762 + MUL15(a
[16], b
[13])
763 + MUL15(a
[17], b
[12])
764 + MUL15(a
[18], b
[11])
765 + MUL15(a
[19], b
[10]);
766 t
[30] = MUL15(a
[11], b
[19])
767 + MUL15(a
[12], b
[18])
768 + MUL15(a
[13], b
[17])
769 + MUL15(a
[14], b
[16])
770 + MUL15(a
[15], b
[15])
771 + MUL15(a
[16], b
[14])
772 + MUL15(a
[17], b
[13])
773 + MUL15(a
[18], b
[12])
774 + MUL15(a
[19], b
[11]);
775 t
[31] = MUL15(a
[12], b
[19])
776 + MUL15(a
[13], b
[18])
777 + MUL15(a
[14], b
[17])
778 + MUL15(a
[15], b
[16])
779 + MUL15(a
[16], b
[15])
780 + MUL15(a
[17], b
[14])
781 + MUL15(a
[18], b
[13])
782 + MUL15(a
[19], b
[12]);
783 t
[32] = MUL15(a
[13], b
[19])
784 + MUL15(a
[14], b
[18])
785 + MUL15(a
[15], b
[17])
786 + MUL15(a
[16], b
[16])
787 + MUL15(a
[17], b
[15])
788 + MUL15(a
[18], b
[14])
789 + MUL15(a
[19], b
[13]);
790 t
[33] = MUL15(a
[14], b
[19])
791 + MUL15(a
[15], b
[18])
792 + MUL15(a
[16], b
[17])
793 + MUL15(a
[17], b
[16])
794 + MUL15(a
[18], b
[15])
795 + MUL15(a
[19], b
[14]);
796 t
[34] = MUL15(a
[15], b
[19])
797 + MUL15(a
[16], b
[18])
798 + MUL15(a
[17], b
[17])
799 + MUL15(a
[18], b
[16])
800 + MUL15(a
[19], b
[15]);
801 t
[35] = MUL15(a
[16], b
[19])
802 + MUL15(a
[17], b
[18])
803 + MUL15(a
[18], b
[17])
804 + MUL15(a
[19], b
[16]);
805 t
[36] = MUL15(a
[17], b
[19])
806 + MUL15(a
[18], b
[18])
807 + MUL15(a
[19], b
[17]);
808 t
[37] = MUL15(a
[18], b
[19])
809 + MUL15(a
[19], b
[18]);
810 t
[38] = MUL15(a
[19], b
[19]);
812 d
[39] = norm13(d
, t
, 39);
816 square20(uint32_t *d
, const uint32_t *a
)
820 t
[ 0] = MUL15(a
[ 0], a
[ 0]);
821 t
[ 1] = ((MUL15(a
[ 0], a
[ 1])) << 1);
822 t
[ 2] = MUL15(a
[ 1], a
[ 1])
823 + ((MUL15(a
[ 0], a
[ 2])) << 1);
824 t
[ 3] = ((MUL15(a
[ 0], a
[ 3])
825 + MUL15(a
[ 1], a
[ 2])) << 1);
826 t
[ 4] = MUL15(a
[ 2], a
[ 2])
827 + ((MUL15(a
[ 0], a
[ 4])
828 + MUL15(a
[ 1], a
[ 3])) << 1);
829 t
[ 5] = ((MUL15(a
[ 0], a
[ 5])
830 + MUL15(a
[ 1], a
[ 4])
831 + MUL15(a
[ 2], a
[ 3])) << 1);
832 t
[ 6] = MUL15(a
[ 3], a
[ 3])
833 + ((MUL15(a
[ 0], a
[ 6])
834 + MUL15(a
[ 1], a
[ 5])
835 + MUL15(a
[ 2], a
[ 4])) << 1);
836 t
[ 7] = ((MUL15(a
[ 0], a
[ 7])
837 + MUL15(a
[ 1], a
[ 6])
838 + MUL15(a
[ 2], a
[ 5])
839 + MUL15(a
[ 3], a
[ 4])) << 1);
840 t
[ 8] = MUL15(a
[ 4], a
[ 4])
841 + ((MUL15(a
[ 0], a
[ 8])
842 + MUL15(a
[ 1], a
[ 7])
843 + MUL15(a
[ 2], a
[ 6])
844 + MUL15(a
[ 3], a
[ 5])) << 1);
845 t
[ 9] = ((MUL15(a
[ 0], a
[ 9])
846 + MUL15(a
[ 1], a
[ 8])
847 + MUL15(a
[ 2], a
[ 7])
848 + MUL15(a
[ 3], a
[ 6])
849 + MUL15(a
[ 4], a
[ 5])) << 1);
850 t
[10] = MUL15(a
[ 5], a
[ 5])
851 + ((MUL15(a
[ 0], a
[10])
852 + MUL15(a
[ 1], a
[ 9])
853 + MUL15(a
[ 2], a
[ 8])
854 + MUL15(a
[ 3], a
[ 7])
855 + MUL15(a
[ 4], a
[ 6])) << 1);
856 t
[11] = ((MUL15(a
[ 0], a
[11])
857 + MUL15(a
[ 1], a
[10])
858 + MUL15(a
[ 2], a
[ 9])
859 + MUL15(a
[ 3], a
[ 8])
860 + MUL15(a
[ 4], a
[ 7])
861 + MUL15(a
[ 5], a
[ 6])) << 1);
862 t
[12] = MUL15(a
[ 6], a
[ 6])
863 + ((MUL15(a
[ 0], a
[12])
864 + MUL15(a
[ 1], a
[11])
865 + MUL15(a
[ 2], a
[10])
866 + MUL15(a
[ 3], a
[ 9])
867 + MUL15(a
[ 4], a
[ 8])
868 + MUL15(a
[ 5], a
[ 7])) << 1);
869 t
[13] = ((MUL15(a
[ 0], a
[13])
870 + MUL15(a
[ 1], a
[12])
871 + MUL15(a
[ 2], a
[11])
872 + MUL15(a
[ 3], a
[10])
873 + MUL15(a
[ 4], a
[ 9])
874 + MUL15(a
[ 5], a
[ 8])
875 + MUL15(a
[ 6], a
[ 7])) << 1);
876 t
[14] = MUL15(a
[ 7], a
[ 7])
877 + ((MUL15(a
[ 0], a
[14])
878 + MUL15(a
[ 1], a
[13])
879 + MUL15(a
[ 2], a
[12])
880 + MUL15(a
[ 3], a
[11])
881 + MUL15(a
[ 4], a
[10])
882 + MUL15(a
[ 5], a
[ 9])
883 + MUL15(a
[ 6], a
[ 8])) << 1);
884 t
[15] = ((MUL15(a
[ 0], a
[15])
885 + MUL15(a
[ 1], a
[14])
886 + MUL15(a
[ 2], a
[13])
887 + MUL15(a
[ 3], a
[12])
888 + MUL15(a
[ 4], a
[11])
889 + MUL15(a
[ 5], a
[10])
890 + MUL15(a
[ 6], a
[ 9])
891 + MUL15(a
[ 7], a
[ 8])) << 1);
892 t
[16] = MUL15(a
[ 8], a
[ 8])
893 + ((MUL15(a
[ 0], a
[16])
894 + MUL15(a
[ 1], a
[15])
895 + MUL15(a
[ 2], a
[14])
896 + MUL15(a
[ 3], a
[13])
897 + MUL15(a
[ 4], a
[12])
898 + MUL15(a
[ 5], a
[11])
899 + MUL15(a
[ 6], a
[10])
900 + MUL15(a
[ 7], a
[ 9])) << 1);
901 t
[17] = ((MUL15(a
[ 0], a
[17])
902 + MUL15(a
[ 1], a
[16])
903 + MUL15(a
[ 2], a
[15])
904 + MUL15(a
[ 3], a
[14])
905 + MUL15(a
[ 4], a
[13])
906 + MUL15(a
[ 5], a
[12])
907 + MUL15(a
[ 6], a
[11])
908 + MUL15(a
[ 7], a
[10])
909 + MUL15(a
[ 8], a
[ 9])) << 1);
910 t
[18] = MUL15(a
[ 9], a
[ 9])
911 + ((MUL15(a
[ 0], a
[18])
912 + MUL15(a
[ 1], a
[17])
913 + MUL15(a
[ 2], a
[16])
914 + MUL15(a
[ 3], a
[15])
915 + MUL15(a
[ 4], a
[14])
916 + MUL15(a
[ 5], a
[13])
917 + MUL15(a
[ 6], a
[12])
918 + MUL15(a
[ 7], a
[11])
919 + MUL15(a
[ 8], a
[10])) << 1);
920 t
[19] = ((MUL15(a
[ 0], a
[19])
921 + MUL15(a
[ 1], a
[18])
922 + MUL15(a
[ 2], a
[17])
923 + MUL15(a
[ 3], a
[16])
924 + MUL15(a
[ 4], a
[15])
925 + MUL15(a
[ 5], a
[14])
926 + MUL15(a
[ 6], a
[13])
927 + MUL15(a
[ 7], a
[12])
928 + MUL15(a
[ 8], a
[11])
929 + MUL15(a
[ 9], a
[10])) << 1);
930 t
[20] = MUL15(a
[10], a
[10])
931 + ((MUL15(a
[ 1], a
[19])
932 + MUL15(a
[ 2], a
[18])
933 + MUL15(a
[ 3], a
[17])
934 + MUL15(a
[ 4], a
[16])
935 + MUL15(a
[ 5], a
[15])
936 + MUL15(a
[ 6], a
[14])
937 + MUL15(a
[ 7], a
[13])
938 + MUL15(a
[ 8], a
[12])
939 + MUL15(a
[ 9], a
[11])) << 1);
940 t
[21] = ((MUL15(a
[ 2], a
[19])
941 + MUL15(a
[ 3], a
[18])
942 + MUL15(a
[ 4], a
[17])
943 + MUL15(a
[ 5], a
[16])
944 + MUL15(a
[ 6], a
[15])
945 + MUL15(a
[ 7], a
[14])
946 + MUL15(a
[ 8], a
[13])
947 + MUL15(a
[ 9], a
[12])
948 + MUL15(a
[10], a
[11])) << 1);
949 t
[22] = MUL15(a
[11], a
[11])
950 + ((MUL15(a
[ 3], a
[19])
951 + MUL15(a
[ 4], a
[18])
952 + MUL15(a
[ 5], a
[17])
953 + MUL15(a
[ 6], a
[16])
954 + MUL15(a
[ 7], a
[15])
955 + MUL15(a
[ 8], a
[14])
956 + MUL15(a
[ 9], a
[13])
957 + MUL15(a
[10], a
[12])) << 1);
958 t
[23] = ((MUL15(a
[ 4], a
[19])
959 + MUL15(a
[ 5], a
[18])
960 + MUL15(a
[ 6], a
[17])
961 + MUL15(a
[ 7], a
[16])
962 + MUL15(a
[ 8], a
[15])
963 + MUL15(a
[ 9], a
[14])
964 + MUL15(a
[10], a
[13])
965 + MUL15(a
[11], a
[12])) << 1);
966 t
[24] = MUL15(a
[12], a
[12])
967 + ((MUL15(a
[ 5], a
[19])
968 + MUL15(a
[ 6], a
[18])
969 + MUL15(a
[ 7], a
[17])
970 + MUL15(a
[ 8], a
[16])
971 + MUL15(a
[ 9], a
[15])
972 + MUL15(a
[10], a
[14])
973 + MUL15(a
[11], a
[13])) << 1);
974 t
[25] = ((MUL15(a
[ 6], a
[19])
975 + MUL15(a
[ 7], a
[18])
976 + MUL15(a
[ 8], a
[17])
977 + MUL15(a
[ 9], a
[16])
978 + MUL15(a
[10], a
[15])
979 + MUL15(a
[11], a
[14])
980 + MUL15(a
[12], a
[13])) << 1);
981 t
[26] = MUL15(a
[13], a
[13])
982 + ((MUL15(a
[ 7], a
[19])
983 + MUL15(a
[ 8], a
[18])
984 + MUL15(a
[ 9], a
[17])
985 + MUL15(a
[10], a
[16])
986 + MUL15(a
[11], a
[15])
987 + MUL15(a
[12], a
[14])) << 1);
988 t
[27] = ((MUL15(a
[ 8], a
[19])
989 + MUL15(a
[ 9], a
[18])
990 + MUL15(a
[10], a
[17])
991 + MUL15(a
[11], a
[16])
992 + MUL15(a
[12], a
[15])
993 + MUL15(a
[13], a
[14])) << 1);
994 t
[28] = MUL15(a
[14], a
[14])
995 + ((MUL15(a
[ 9], a
[19])
996 + MUL15(a
[10], a
[18])
997 + MUL15(a
[11], a
[17])
998 + MUL15(a
[12], a
[16])
999 + MUL15(a
[13], a
[15])) << 1);
1000 t
[29] = ((MUL15(a
[10], a
[19])
1001 + MUL15(a
[11], a
[18])
1002 + MUL15(a
[12], a
[17])
1003 + MUL15(a
[13], a
[16])
1004 + MUL15(a
[14], a
[15])) << 1);
1005 t
[30] = MUL15(a
[15], a
[15])
1006 + ((MUL15(a
[11], a
[19])
1007 + MUL15(a
[12], a
[18])
1008 + MUL15(a
[13], a
[17])
1009 + MUL15(a
[14], a
[16])) << 1);
1010 t
[31] = ((MUL15(a
[12], a
[19])
1011 + MUL15(a
[13], a
[18])
1012 + MUL15(a
[14], a
[17])
1013 + MUL15(a
[15], a
[16])) << 1);
1014 t
[32] = MUL15(a
[16], a
[16])
1015 + ((MUL15(a
[13], a
[19])
1016 + MUL15(a
[14], a
[18])
1017 + MUL15(a
[15], a
[17])) << 1);
1018 t
[33] = ((MUL15(a
[14], a
[19])
1019 + MUL15(a
[15], a
[18])
1020 + MUL15(a
[16], a
[17])) << 1);
1021 t
[34] = MUL15(a
[17], a
[17])
1022 + ((MUL15(a
[15], a
[19])
1023 + MUL15(a
[16], a
[18])) << 1);
1024 t
[35] = ((MUL15(a
[16], a
[19])
1025 + MUL15(a
[17], a
[18])) << 1);
1026 t
[36] = MUL15(a
[18], a
[18])
1027 + ((MUL15(a
[17], a
[19])) << 1);
1028 t
[37] = ((MUL15(a
[18], a
[19])) << 1);
1029 t
[38] = MUL15(a
[19], a
[19]);
1031 d
[39] = norm13(d
, t
, 39);
1037 * Perform a "final reduction" in field F255 (field for Curve25519)
1038 * The source value must be less than twice the modulus. If the value
1039 * is not lower than the modulus, then the modulus is subtracted and
1040 * this function returns 1; otherwise, it leaves it untouched and it
1044 reduce_final_f255(uint32_t *d
)
1050 memcpy(t
, d
, sizeof t
);
1052 for (i
= 0; i
< 20; i
++) {
1061 CCOPY(cc
, d
, t
, sizeof t
);
1066 f255_mulgen(uint32_t *d
, const uint32_t *a
, const uint32_t *b
, int square
)
1068 uint32_t t
[40], cc
, w
;
1071 * Compute raw multiplication. All result words fit in 13 bits
1072 * each; upper word (t[39]) must fit on 5 bits, since the product
1073 * of two 256-bit integers must fit on 512 bits.
1082 * Modular reduction: each high word is added where necessary.
1083 * Since the modulus is 2^255-19 and word 20 corresponds to
1084 * offset 20*13 = 260, word 20+k must be added to word k with
1085 * a factor of 19*2^5 = 608. The extra bits in word 19 are also
1088 cc
= MUL15(t
[19] >> 8, 19);
1091 #define MM1(x) do { \
1092 w = t[x] + cc + MUL15(t[(x) + 20], 608); \
1093 t[x] = w & 0x1FFF; \
1120 cc
= MUL15(w
>> 8, 19);
1123 #define MM2(x) do { \
1125 d[x] = w & 0x1FFF; \
1154 * Perform a multiplication of two integers modulo 2^255-19.
1155 * Operands are arrays of 20 words, each containing 13 bits of data, in
1156 * little-endian order. Input value may be up to 2^256-1; on output, value
1157 * fits on 256 bits and is lower than twice the modulus.
1159 * f255_mul() is the general multiplication, f255_square() is specialised
1162 #define f255_mul(d, a, b) f255_mulgen(d, a, b, 0)
1163 #define f255_square(d, a) f255_mulgen(d, a, a, 1)
1166 * Add two values in F255. Partial reduction is performed (down to less
1167 * than twice the modulus).
1170 f255_add(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
1176 for (i
= 0; i
< 20; i
++) {
1177 w
= a
[i
] + b
[i
] + cc
;
1181 cc
= MUL15(w
>> 8, 19);
1183 for (i
= 0; i
< 20; i
++) {
1191 * Subtract one value from another in F255. Partial reduction is
1192 * performed (down to less than twice the modulus).
1195 f255_sub(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
1198 * We actually compute a - b + 2*p, so that the final value is
1199 * necessarily positive.
1205 for (i
= 0; i
< 20; i
++) {
1206 w
= a
[i
] - b
[i
] + cc
;
1210 cc
= MUL15((w
+ 0x200) >> 8, 19);
1212 for (i
= 0; i
< 20; i
++) {
1220 * Multiply an integer by the 'A24' constant (121665). Partial reduction
1221 * is performed (down to less than twice the modulus).
1224 f255_mul_a24(uint32_t *d
, const uint32_t *a
)
1230 for (i
= 0; i
< 20; i
++) {
1231 w
= MUL15(a
[i
], 121665) + cc
;
1235 cc
= MUL15(w
>> 8, 19);
1237 for (i
= 0; i
< 20; i
++) {
1244 static const unsigned char GEN
[] = {
1245 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1246 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1247 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1248 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
1251 static const unsigned char ORDER
[] = {
1252 0x7F, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
1253 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
1254 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
1255 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
1258 static const unsigned char *
1259 api_generator(int curve
, size_t *len
)
1266 static const unsigned char *
1267 api_order(int curve
, size_t *len
)
1275 api_xoff(int curve
, size_t *len
)
1283 cswap(uint32_t *a
, uint32_t *b
, uint32_t ctl
)
1288 for (i
= 0; i
< 20; i
++) {
1289 uint32_t aw
, bw
, tw
;
1293 tw
= ctl
& (aw
^ bw
);
1300 api_mul(unsigned char *G
, size_t Glen
,
1301 const unsigned char *kb
, size_t kblen
, int curve
)
1303 uint32_t x1
[20], x2
[20], x3
[20], z2
[20], z3
[20];
1304 uint32_t a
[20], aa
[20], b
[20], bb
[20];
1305 uint32_t c
[20], d
[20], e
[20], da
[20], cb
[20];
1306 unsigned char k
[32];
1313 * Points are encoded over exactly 32 bytes. Multipliers must fit
1314 * in 32 bytes as well.
1315 * RFC 7748 mandates that the high bit of the last point byte must
1316 * be ignored/cleared.
1318 if (Glen
!= 32 || kblen
> 32) {
1324 * Initialise variables x1, x2, z2, x3 and z3. We set all of them
1325 * into Montgomery representation.
1327 x1
[19] = le8_to_le13(x1
, G
, 32);
1328 memcpy(x3
, x1
, sizeof x1
);
1329 memset(z2
, 0, sizeof z2
);
1330 memset(x2
, 0, sizeof x2
);
1332 memset(z3
, 0, sizeof z3
);
1335 memcpy(k
, kb
, kblen
);
1336 memset(k
+ kblen
, 0, (sizeof k
) - kblen
);
1342 print_int("x1", x1);
1346 for (i
= 254; i
>= 0; i
--) {
1349 kt
= (k
[i
>> 3] >> (i
& 7)) & 1;
1351 cswap(x2
, x3
, swap
);
1352 cswap(z2
, z3
, swap
);
1356 print_int("x2", x2);
1357 print_int("z2", z2);
1358 print_int("x3", x3);
1359 print_int("z3", z3);
1362 f255_add(a
, x2
, z2
);
1364 f255_sub(b
, x2
, z2
);
1366 f255_sub(e
, aa
, bb
);
1367 f255_add(c
, x3
, z3
);
1368 f255_sub(d
, x3
, z3
);
1374 print_int("aa", aa);
1376 print_int("bb", bb);
1380 print_int("da", da);
1381 print_int("cb", cb);
1384 f255_add(x3
, da
, cb
);
1385 f255_square(x3
, x3
);
1386 f255_sub(z3
, da
, cb
);
1387 f255_square(z3
, z3
);
1388 f255_mul(z3
, z3
, x1
);
1389 f255_mul(x2
, aa
, bb
);
1390 f255_mul_a24(z2
, e
);
1391 f255_add(z2
, z2
, aa
);
1392 f255_mul(z2
, e
, z2
);
1395 print_int("x2", x2);
1396 print_int("z2", z2);
1397 print_int("x3", x3);
1398 print_int("z3", z3);
1401 cswap(x2
, x3
, swap
);
1402 cswap(z2
, z3
, swap
);
1405 * Inverse z2 with a modular exponentiation. This is a simple
1406 * square-and-multiply algorithm; we mutualise most non-squarings
1407 * since the exponent contains almost only ones.
1409 memcpy(a
, z2
, sizeof z2
);
1410 for (i
= 0; i
< 15; i
++) {
1414 memcpy(b
, a
, sizeof a
);
1415 for (i
= 0; i
< 14; i
++) {
1418 for (j
= 0; j
< 16; j
++) {
1423 for (i
= 14; i
>= 0; i
--) {
1425 if ((0xFFEB >> i
) & 1) {
1429 f255_mul(x2
, x2
, b
);
1430 reduce_final_f255(x2
);
1431 le13_to_le8(G
, 32, x2
);
1436 api_mulgen(unsigned char *R
,
1437 const unsigned char *x
, size_t xlen
, int curve
)
1439 const unsigned char *G
;
1442 G
= api_generator(curve
, &Glen
);
1444 api_mul(R
, Glen
, x
, xlen
, curve
);
1449 api_muladd(unsigned char *A
, const unsigned char *B
, size_t len
,
1450 const unsigned char *x
, size_t xlen
,
1451 const unsigned char *y
, size_t ylen
, int curve
)
1454 * We don't implement this method, since it is used for ECDSA
1455 * only, and there is no ECDSA over Curve25519 (which instead
1469 /* see bearssl_ec.h */
1470 const br_ec_impl br_ec_c25519_m15
= {
1471 (uint32_t)0x20000000,