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[40];
36 printf("%s = ", name);
37 for (u = 0; u < 9; u ++) {
38 if (x[u] > 0x3FFFFFFF) {
40 for (u = 0; u < 9; u ++) {
41 printf(" %08X", x[u]);
47 memset(tmp, 0, sizeof tmp);
48 for (u = 0; u < 9; u ++) {
60 for (j = 0; j < 8; j ++) {
61 tmp[39 - k - j] |= (unsigned char)w;
65 for (u = 8; u < 40; 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 * 30-bit words in little-endian order. The final "partial" word is
95 le8_to_le30(uint32_t *dst
, const unsigned char *src
, size_t len
)
110 *dst
++ = (acc
| (b
<< acc_len
)) & 0x3FFFFFFF;
111 acc
= b
>> (30 - acc_len
);
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.
124 le30_to_le8(unsigned char *dst
, size_t len
, const uint32_t *src
)
136 *dst
++ = (unsigned char)(acc
| (w
<< acc_len
));
137 acc
= w
>> (8 - acc_len
);
140 *dst
++ = (unsigned char)acc
;
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.
153 mul9(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
156 * Maximum intermediate result is no more than
157 * 10376293531797946367, which fits in 64 bits. Reason:
159 * 10376293531797946367 = 9 * (2^30-1)^2 + 9663676406
160 * 10376293531797946367 < 9663676407 * 2^30
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.
170 t
[ 0] = MUL31(a
[0], b
[0]);
171 t
[ 1] = MUL31(a
[0], b
[1])
173 t
[ 2] = MUL31(a
[0], b
[2])
176 t
[ 3] = MUL31(a
[0], b
[3])
180 t
[ 4] = MUL31(a
[0], b
[4])
185 t
[ 5] = MUL31(a
[0], b
[5])
191 t
[ 6] = MUL31(a
[0], b
[6])
198 t
[ 7] = MUL31(a
[0], b
[7])
206 t
[ 8] = MUL31(a
[0], b
[8])
215 t
[ 9] = MUL31(a
[1], b
[8])
223 t
[10] = MUL31(a
[2], b
[8])
230 t
[11] = MUL31(a
[3], b
[8])
236 t
[12] = MUL31(a
[4], b
[8])
241 t
[13] = MUL31(a
[5], b
[8])
245 t
[14] = MUL31(a
[6], b
[8])
248 t
[15] = MUL31(a
[7], b
[8])
250 t
[16] = MUL31(a
[8], b
[8]);
256 for (i
= 0; i
< 17; i
++) {
260 d
[i
] = (uint32_t)w
& 0x3FFFFFFF;
263 d
[17] = (uint32_t)cc
;
267 * Square a 270-bit integer, represented as an array of nine 30-bit words.
268 * Result uses 18 words of 30 bits each.
271 square9(uint32_t *d
, const uint32_t *a
)
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])
288 + MUL31(a
[2], a
[3])) << 1);
289 t
[ 6] = MUL31(a
[3], a
[3])
290 + ((MUL31(a
[0], a
[6])
292 + MUL31(a
[2], a
[4])) << 1);
293 t
[ 7] = ((MUL31(a
[0], a
[7])
296 + MUL31(a
[3], a
[4])) << 1);
297 t
[ 8] = MUL31(a
[4], a
[4])
298 + ((MUL31(a
[0], a
[8])
301 + MUL31(a
[3], a
[5])) << 1);
302 t
[ 9] = ((MUL31(a
[1], a
[8])
305 + MUL31(a
[4], a
[5])) << 1);
306 t
[10] = MUL31(a
[5], a
[5])
307 + ((MUL31(a
[2], a
[8])
309 + MUL31(a
[4], a
[6])) << 1);
310 t
[11] = ((MUL31(a
[3], a
[8])
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]);
327 for (i
= 0; i
< 17; i
++) {
331 d
[i
] = (uint32_t)w
& 0x3FFFFFFF;
334 d
[17] = (uint32_t)cc
;
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
345 reduce_final_f255(uint32_t *d
)
351 memcpy(t
, d
, sizeof t
);
353 for (i
= 0; i
< 9; i
++) {
358 t
[i
] = w
& 0x3FFFFFFF;
362 CCOPY(cc
, d
, t
, sizeof t
);
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.
373 f255_mul(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
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.
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
392 * Keeping the carry on 32 bits helps with 32-bit architectures,
393 * and does not noticeably impact performance on 64-bit systems.
395 cc
= MUL15(t
[8] >> 15, 19); /* at most 19*(2^15-1) = 622573 */
397 for (i
= 0; i
< 9; i
++) {
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 */
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.
413 cc
= MUL15(t
[8] >> 15, 19);
415 for (i
= 0; i
< 9; i
++) {
419 d
[i
] = z
& 0x3FFFFFFF;
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.
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.
436 f255_square(uint32_t *d
, const uint32_t *a
)
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.
449 * Modular reduction: each high word is added where necessary.
450 * See f255_mul() for details on the reduction and carry limits.
452 cc
= MUL15(t
[8] >> 15, 19);
454 for (i
= 0; i
< 9; i
++) {
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);
461 cc
= MUL15(t
[8] >> 15, 19);
463 for (i
= 0; i
< 9; i
++) {
467 d
[i
] = z
& 0x3FFFFFFF;
473 * Add two values in F255. Partial reduction is performed (down to less
474 * than twice the modulus).
477 f255_add(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
480 * Since operand words fit on 30 bits, we can use 32-bit
481 * variables throughout.
487 for (i
= 0; i
< 9; i
++) {
488 w
= a
[i
] + b
[i
] + cc
;
489 d
[i
] = w
& 0x3FFFFFFF;
492 cc
= MUL15(w
>> 15, 19);
494 for (i
= 0; i
< 9; i
++) {
496 d
[i
] = w
& 0x3FFFFFFF;
502 * Subtract one value from another in F255. Partial reduction is
503 * performed (down to less than twice the modulus).
506 f255_sub(uint32_t *d
, const uint32_t *a
, const uint32_t *b
)
509 * We actually compute a - b + 2*p, so that the final value is
510 * necessarily positive.
516 for (i
= 0; i
< 9; i
++) {
517 w
= a
[i
] - b
[i
] + cc
;
518 d
[i
] = w
& 0x3FFFFFFF;
521 cc
= MUL15((w
+ 0x10000) >> 15, 19);
523 for (i
= 0; i
< 9; i
++) {
525 d
[i
] = w
& 0x3FFFFFFF;
531 * Multiply an integer by the 'A24' constant (121665). Partial reduction
532 * is performed (down to less than twice the modulus).
535 f255_mul_a24(uint32_t *d
, const uint32_t *a
)
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.
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);
553 w
= MUL31(a
[8], 121665) + (uint64_t)cc
;
554 d
[8] = (uint32_t)w
& 0x7FFF;
555 cc
= MUL15((uint32_t)(w
>> 15), 19);
557 for (i
= 0; i
< 9; i
++) {
561 d
[i
] = z
& 0x3FFFFFFF;
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
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
580 static const unsigned char *
581 api_generator(int curve
, size_t *len
)
588 static const unsigned char *
589 api_order(int curve
, size_t *len
)
597 api_xoff(int curve
, size_t *len
)
605 cswap(uint32_t *a
, uint32_t *b
, uint32_t ctl
)
610 for (i
= 0; i
< 9; i
++) {
615 tw
= ctl
& (aw
^ bw
);
622 api_mul(unsigned char *G
, size_t Glen
,
623 const unsigned char *kb
, size_t kblen
, int curve
)
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];
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.
640 if (Glen
!= 32 || kblen
> 32) {
646 * Initialise variables x1, x2, z2, x3 and z3. We set all of them
647 * into Montgomery representation.
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
);
654 memset(z3
, 0, sizeof z3
);
657 memset(k
, 0, (sizeof k
) - kblen
);
658 memcpy(k
+ (sizeof k
) - kblen
, kb
, kblen
);
668 for (i
= 254; i
>= 0; i
--) {
671 kt
= (k
[31 - (i
>> 3)] >> (i
& 7)) & 1;
706 f255_add(x3
, da
, cb
);
708 f255_sub(z3
, da
, cb
);
710 f255_mul(z3
, z3
, x1
);
711 f255_mul(x2
, aa
, bb
);
713 f255_add(z2
, z2
, aa
);
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.
731 memcpy(a
, z2
, sizeof z2
);
732 for (i
= 0; i
< 15; i
++) {
736 memcpy(b
, a
, sizeof a
);
737 for (i
= 0; i
< 14; i
++) {
740 for (j
= 0; j
< 16; j
++) {
745 for (i
= 14; i
>= 0; i
--) {
747 if ((0xFFEB >> i
) & 1) {
752 reduce_final_f255(x2
);
753 le30_to_le8(G
, 32, x2
);
758 api_mulgen(unsigned char *R
,
759 const unsigned char *x
, size_t xlen
, int curve
)
761 const unsigned char *G
;
764 G
= api_generator(curve
, &Glen
);
766 api_mul(R
, Glen
, x
, xlen
, curve
);
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
)
776 * We don't implement this method, since it is used for ECDSA
777 * only, and there is no ECDSA over Curve25519 (which instead
791 /* see bearssl_ec.h */
792 const br_ec_impl br_ec_c25519_m31
= {
793 (uint32_t)0x20000000,