2 * Copyright 2020-2023 The OpenSSL Project Authors. All Rights Reserved.
3 * Copyright (c) 2020-2021, Intel Corporation. All Rights Reserved.
5 * Licensed under the Apache License 2.0 (the "License"). You may not use
6 * this file except in compliance with the License. You can obtain a copy
7 * in the file LICENSE in the source distribution or at
8 * https://www.openssl.org/source/license.html
11 * Originally written by Sergey Kirillov and Andrey Matyukov.
12 * Special thanks to Ilya Albrekht for his valuable hints.
17 #include <openssl/opensslconf.h>
18 #include <openssl/crypto.h>
22 NON_EMPTY_TRANSLATION_UNIT
27 # define ALIGN_OF(ptr, boundary) \
28 ((unsigned char *)(ptr) + (boundary - (((size_t)(ptr)) & (boundary - 1))))
31 # define DIGIT_SIZE (52)
33 # define DIGIT_MASK ((uint64_t)0xFFFFFFFFFFFFF)
35 # define BITS2WORD8_SIZE(x) (((x) + 7) >> 3)
36 # define BITS2WORD64_SIZE(x) (((x) + 63) >> 6)
38 /* Number of registers required to hold |digits_num| amount of qword digits */
39 # define NUMBER_OF_REGISTERS(digits_num, register_size) \
40 (((digits_num) * 64 + (register_size) - 1) / (register_size))
42 static ossl_inline
uint64_t get_digit(const uint8_t *in
, int in_len
);
43 static ossl_inline
void put_digit(uint8_t *out
, int out_len
, uint64_t digit
);
44 static void to_words52(BN_ULONG
*out
, int out_len
, const BN_ULONG
*in
,
46 static void from_words52(BN_ULONG
*bn_out
, int out_bitsize
, const BN_ULONG
*in
);
47 static ossl_inline
void set_bit(BN_ULONG
*a
, int idx
);
49 /* Number of |digit_size|-bit digits in |bitsize|-bit value */
50 static ossl_inline
int number_of_digits(int bitsize
, int digit_size
)
52 return (bitsize
+ digit_size
- 1) / digit_size
;
56 * For details of the methods declared below please refer to
57 * crypto/bn/asm/rsaz-avx512.pl
60 * amm = Almost Montgomery Multiplication
61 * ams = Almost Montgomery Squaring
62 * 52xZZ - data represented as array of ZZ digits in 52-bit radix
63 * _x1_/_x2_ - 1 or 2 independent inputs/outputs
64 * _ifma256 - uses 256-bit wide IFMA ISA (AVX512_IFMA256)
67 void ossl_rsaz_amm52x20_x1_ifma256(BN_ULONG
*res
, const BN_ULONG
*a
,
68 const BN_ULONG
*b
, const BN_ULONG
*m
,
70 void ossl_rsaz_amm52x20_x2_ifma256(BN_ULONG
*out
, const BN_ULONG
*a
,
71 const BN_ULONG
*b
, const BN_ULONG
*m
,
72 const BN_ULONG k0
[2]);
73 void ossl_extract_multiplier_2x20_win5(BN_ULONG
*red_Y
,
74 const BN_ULONG
*red_table
,
75 int red_table_idx1
, int red_table_idx2
);
77 void ossl_rsaz_amm52x30_x1_ifma256(BN_ULONG
*res
, const BN_ULONG
*a
,
78 const BN_ULONG
*b
, const BN_ULONG
*m
,
80 void ossl_rsaz_amm52x30_x2_ifma256(BN_ULONG
*out
, const BN_ULONG
*a
,
81 const BN_ULONG
*b
, const BN_ULONG
*m
,
82 const BN_ULONG k0
[2]);
83 void ossl_extract_multiplier_2x30_win5(BN_ULONG
*red_Y
,
84 const BN_ULONG
*red_table
,
85 int red_table_idx1
, int red_table_idx2
);
87 void ossl_rsaz_amm52x40_x1_ifma256(BN_ULONG
*res
, const BN_ULONG
*a
,
88 const BN_ULONG
*b
, const BN_ULONG
*m
,
90 void ossl_rsaz_amm52x40_x2_ifma256(BN_ULONG
*out
, const BN_ULONG
*a
,
91 const BN_ULONG
*b
, const BN_ULONG
*m
,
92 const BN_ULONG k0
[2]);
93 void ossl_extract_multiplier_2x40_win5(BN_ULONG
*red_Y
,
94 const BN_ULONG
*red_table
,
95 int red_table_idx1
, int red_table_idx2
);
97 static int RSAZ_mod_exp_x2_ifma256(BN_ULONG
*res
, const BN_ULONG
*base
,
98 const BN_ULONG
*exp
[2], const BN_ULONG
*m
,
99 const BN_ULONG
*rr
, const BN_ULONG k0
[2],
100 int modulus_bitsize
);
103 * Dual Montgomery modular exponentiation using prime moduli of the
104 * same bit size, optimized with AVX512 ISA.
106 * Input and output parameters for each exponentiation are independent and
107 * denoted here by index |i|, i = 1..2.
109 * Input and output are all in regular 2^64 radix.
111 * Each moduli shall be |factor_size| bit size.
118 * [out] res|i| - result of modular exponentiation: array of qword values
119 * in regular (2^64) radix. Size of array shall be enough
120 * to hold |factor_size| bits.
121 * [in] base|i| - base
122 * [in] exp|i| - exponent
124 * [in] rr|i| - Montgomery parameter RR = R^2 mod m|i|
125 * [in] k0_|i| - Montgomery parameter k0 = -1/m|i| mod 2^64
126 * [in] factor_size - moduli bit size
128 * \return 0 in case of failure,
129 * 1 in case of success.
131 int ossl_rsaz_mod_exp_avx512_x2(BN_ULONG
*res1
,
132 const BN_ULONG
*base1
,
133 const BN_ULONG
*exp1
,
138 const BN_ULONG
*base2
,
139 const BN_ULONG
*exp2
,
145 typedef void (*AMM
)(BN_ULONG
*res
, const BN_ULONG
*a
,
146 const BN_ULONG
*b
, const BN_ULONG
*m
, BN_ULONG k0
);
150 * Number of word-size (BN_ULONG) digits to store exponent in redundant
153 int exp_digits
= number_of_digits(factor_size
+ 2, DIGIT_SIZE
);
154 int coeff_pow
= 4 * (DIGIT_SIZE
* exp_digits
- factor_size
);
156 /* Number of YMM registers required to store exponent's digits */
157 int ymm_regs_num
= NUMBER_OF_REGISTERS(exp_digits
, 256 /* ymm bit size */);
158 /* Capacity of the register set (in qwords) to store exponent */
159 int regs_capacity
= ymm_regs_num
* 4;
161 BN_ULONG
*base1_red
, *m1_red
, *rr1_red
;
162 BN_ULONG
*base2_red
, *m2_red
, *rr2_red
;
164 BN_ULONG
*storage
= NULL
;
165 BN_ULONG
*storage_aligned
= NULL
;
166 int storage_len_bytes
= 7 * regs_capacity
* sizeof(BN_ULONG
)
167 + 64 /* alignment */;
169 const BN_ULONG
*exp
[2] = {0};
170 BN_ULONG k0
[2] = {0};
171 /* AMM = Almost Montgomery Multiplication */
174 switch (factor_size
) {
176 amm
= ossl_rsaz_amm52x20_x1_ifma256
;
179 amm
= ossl_rsaz_amm52x30_x1_ifma256
;
182 amm
= ossl_rsaz_amm52x40_x1_ifma256
;
188 storage
= (BN_ULONG
*)OPENSSL_malloc(storage_len_bytes
);
191 storage_aligned
= (BN_ULONG
*)ALIGN_OF(storage
, 64);
193 /* Memory layout for red(undant) representations */
194 base1_red
= storage_aligned
;
195 base2_red
= storage_aligned
+ 1 * regs_capacity
;
196 m1_red
= storage_aligned
+ 2 * regs_capacity
;
197 m2_red
= storage_aligned
+ 3 * regs_capacity
;
198 rr1_red
= storage_aligned
+ 4 * regs_capacity
;
199 rr2_red
= storage_aligned
+ 5 * regs_capacity
;
200 coeff_red
= storage_aligned
+ 6 * regs_capacity
;
202 /* Convert base_i, m_i, rr_i, from regular to 52-bit radix */
203 to_words52(base1_red
, regs_capacity
, base1
, factor_size
);
204 to_words52(base2_red
, regs_capacity
, base2
, factor_size
);
205 to_words52(m1_red
, regs_capacity
, m1
, factor_size
);
206 to_words52(m2_red
, regs_capacity
, m2
, factor_size
);
207 to_words52(rr1_red
, regs_capacity
, rr1
, factor_size
);
208 to_words52(rr2_red
, regs_capacity
, rr2
, factor_size
);
211 * Compute target domain Montgomery converters RR' for each modulus
212 * based on precomputed original domain's RR.
214 * RR -> RR' transformation steps:
216 * (2) t = AMM(RR,RR) = RR^2 / R' mod m
217 * (3) RR' = AMM(t, coeff) = RR^2 * 2^k / R'^2 mod m
219 * k = 4 * (52 * digits52 - modlen)
220 * R = 2^(64 * ceil(modlen/64)) mod m
222 * R' = 2^(52 * ceil(modlen/52)) mod m
224 * EX/ modlen = 1024: k = 64, RR = 2^2048 mod m, RR' = 2^2080 mod m
226 memset(coeff_red
, 0, exp_digits
* sizeof(BN_ULONG
));
227 /* (1) in reduced domain representation */
228 set_bit(coeff_red
, 64 * (int)(coeff_pow
/ 52) + coeff_pow
% 52);
230 amm(rr1_red
, rr1_red
, rr1_red
, m1_red
, k0_1
); /* (2) for m1 */
231 amm(rr1_red
, rr1_red
, coeff_red
, m1_red
, k0_1
); /* (3) for m1 */
233 amm(rr2_red
, rr2_red
, rr2_red
, m2_red
, k0_2
); /* (2) for m2 */
234 amm(rr2_red
, rr2_red
, coeff_red
, m2_red
, k0_2
); /* (3) for m2 */
242 /* Dual (2-exps in parallel) exponentiation */
243 ret
= RSAZ_mod_exp_x2_ifma256(rr1_red
, base1_red
, exp
, m1_red
, rr1_red
,
248 /* Convert rr_i back to regular radix */
249 from_words52(res1
, factor_size
, rr1_red
);
250 from_words52(res2
, factor_size
, rr2_red
);
252 /* bn_reduce_once_in_place expects number of BN_ULONG, not bit size */
253 factor_size
/= sizeof(BN_ULONG
) * 8;
255 bn_reduce_once_in_place(res1
, /*carry=*/0, m1
, storage
, factor_size
);
256 bn_reduce_once_in_place(res2
, /*carry=*/0, m2
, storage
, factor_size
);
259 if (storage
!= NULL
) {
260 OPENSSL_cleanse(storage
, storage_len_bytes
);
261 OPENSSL_free(storage
);
267 * Dual {1024,1536,2048}-bit w-ary modular exponentiation using prime moduli of
268 * the same bit size using Almost Montgomery Multiplication, optimized with
269 * AVX512_IFMA256 ISA.
271 * The parameter w (window size) = 5.
273 * [out] res - result of modular exponentiation: 2x{20,30,40} qword
274 * values in 2^52 radix.
275 * [in] base - base (2x{20,30,40} qword values in 2^52 radix)
276 * [in] exp - array of 2 pointers to {16,24,32} qword values in 2^64 radix.
277 * Exponent is not converted to redundant representation.
278 * [in] m - moduli (2x{20,30,40} qword values in 2^52 radix)
279 * [in] rr - Montgomery parameter for 2 moduli:
280 * RR(1024) = 2^2080 mod m.
281 * RR(1536) = 2^3120 mod m.
282 * RR(2048) = 2^4160 mod m.
283 * (2x{20,30,40} qword values in 2^52 radix)
284 * [in] k0 - Montgomery parameter for 2 moduli: k0 = -1/m mod 2^64
288 int RSAZ_mod_exp_x2_ifma256(BN_ULONG
*out
,
289 const BN_ULONG
*base
,
290 const BN_ULONG
*exp
[2],
293 const BN_ULONG k0
[2],
296 typedef void (*DAMM
)(BN_ULONG
*res
, const BN_ULONG
*a
,
297 const BN_ULONG
*b
, const BN_ULONG
*m
,
298 const BN_ULONG k0
[2]);
299 typedef void (*DEXTRACT
)(BN_ULONG
*res
, const BN_ULONG
*red_table
,
300 int red_table_idx
, int tbl_idx
);
305 /* Exponent window size */
306 int exp_win_size
= 5;
307 int exp_win_mask
= (1U << exp_win_size
) - 1;
310 * Number of digits (64-bit words) in redundant representation to handle
316 BN_ULONG
*storage
= NULL
;
317 BN_ULONG
*storage_aligned
= NULL
;
318 int storage_len_bytes
= 0;
320 /* Red(undant) result Y and multiplier X */
321 BN_ULONG
*red_Y
= NULL
; /* [2][red_digits] */
322 BN_ULONG
*red_X
= NULL
; /* [2][red_digits] */
323 /* Pre-computed table of base powers */
324 BN_ULONG
*red_table
= NULL
; /* [1U << exp_win_size][2][red_digits] */
325 /* Expanded exponent */
326 BN_ULONG
*expz
= NULL
; /* [2][exp_digits + 1] */
330 /* Extractor from red_table */
331 DEXTRACT extract
= NULL
;
334 * Squaring is done using multiplication now. That can be a subject of
335 * optimization in future.
337 # define DAMS(r,a,m,k0) damm((r),(a),(a),(m),(k0))
339 switch (modulus_bitsize
) {
343 damm
= ossl_rsaz_amm52x20_x2_ifma256
;
344 extract
= ossl_extract_multiplier_2x20_win5
;
347 /* Extended with 2 digits padding to avoid mask ops in high YMM register */
350 damm
= ossl_rsaz_amm52x30_x2_ifma256
;
351 extract
= ossl_extract_multiplier_2x30_win5
;
356 damm
= ossl_rsaz_amm52x40_x2_ifma256
;
357 extract
= ossl_extract_multiplier_2x40_win5
;
363 storage_len_bytes
= (2 * red_digits
/* red_Y */
364 + 2 * red_digits
/* red_X */
365 + 2 * red_digits
* (1U << exp_win_size
) /* red_table */
366 + 2 * (exp_digits
+ 1)) /* expz */
368 + 64; /* alignment */
370 storage
= (BN_ULONG
*)OPENSSL_zalloc(storage_len_bytes
);
373 storage_aligned
= (BN_ULONG
*)ALIGN_OF(storage
, 64);
375 red_Y
= storage_aligned
;
376 red_X
= red_Y
+ 2 * red_digits
;
377 red_table
= red_X
+ 2 * red_digits
;
378 expz
= red_table
+ 2 * red_digits
* (1U << exp_win_size
);
381 * Compute table of powers base^i, i = 0, ..., (2^EXP_WIN_SIZE) - 1
382 * table[0] = mont(x^0) = mont(1)
383 * table[1] = mont(x^1) = mont(x)
385 red_X
[0 * red_digits
] = 1;
386 red_X
[1 * red_digits
] = 1;
387 damm(&red_table
[0 * 2 * red_digits
], (const BN_ULONG
*)red_X
, rr
, m
, k0
);
388 damm(&red_table
[1 * 2 * red_digits
], base
, rr
, m
, k0
);
390 for (idx
= 1; idx
< (int)((1U << exp_win_size
) / 2); idx
++) {
391 DAMS(&red_table
[(2 * idx
+ 0) * 2 * red_digits
],
392 &red_table
[(1 * idx
) * 2 * red_digits
], m
, k0
);
393 damm(&red_table
[(2 * idx
+ 1) * 2 * red_digits
],
394 &red_table
[(2 * idx
) * 2 * red_digits
],
395 &red_table
[1 * 2 * red_digits
], m
, k0
);
398 /* Copy and expand exponents */
399 memcpy(&expz
[0 * (exp_digits
+ 1)], exp
[0], exp_digits
* sizeof(BN_ULONG
));
400 expz
[1 * (exp_digits
+ 1) - 1] = 0;
401 memcpy(&expz
[1 * (exp_digits
+ 1)], exp
[1], exp_digits
* sizeof(BN_ULONG
));
402 expz
[2 * (exp_digits
+ 1) - 1] = 0;
406 const int rem
= modulus_bitsize
% exp_win_size
;
407 const BN_ULONG table_idx_mask
= exp_win_mask
;
409 int exp_bit_no
= modulus_bitsize
- rem
;
410 int exp_chunk_no
= exp_bit_no
/ 64;
411 int exp_chunk_shift
= exp_bit_no
% 64;
413 BN_ULONG red_table_idx_0
, red_table_idx_1
;
417 * exp_bit_no = modulus_bitsize - exp_win_size
418 * However, this isn't possible because rem is { 1024, 1536, 2048 } % 5
419 * which is { 4, 1, 3 } respectively.
421 * If this assertion ever fails the fix above is easy.
423 OPENSSL_assert(rem
!= 0);
425 /* Process 1-st exp window - just init result */
426 red_table_idx_0
= expz
[exp_chunk_no
+ 0 * (exp_digits
+ 1)];
427 red_table_idx_1
= expz
[exp_chunk_no
+ 1 * (exp_digits
+ 1)];
430 * The function operates with fixed moduli sizes divisible by 64,
431 * thus table index here is always in supported range [0, EXP_WIN_SIZE).
433 red_table_idx_0
>>= exp_chunk_shift
;
434 red_table_idx_1
>>= exp_chunk_shift
;
436 extract(&red_Y
[0 * red_digits
], (const BN_ULONG
*)red_table
, (int)red_table_idx_0
, (int)red_table_idx_1
);
438 /* Process other exp windows */
439 for (exp_bit_no
-= exp_win_size
; exp_bit_no
>= 0; exp_bit_no
-= exp_win_size
) {
440 /* Extract pre-computed multiplier from the table */
444 exp_chunk_no
= exp_bit_no
/ 64;
445 exp_chunk_shift
= exp_bit_no
% 64;
447 red_table_idx_0
= expz
[exp_chunk_no
+ 0 * (exp_digits
+ 1)];
448 T
= expz
[exp_chunk_no
+ 1 + 0 * (exp_digits
+ 1)];
450 red_table_idx_0
>>= exp_chunk_shift
;
452 * Get additional bits from then next quadword
453 * when 64-bit boundaries are crossed.
455 if (exp_chunk_shift
> 64 - exp_win_size
) {
456 T
<<= (64 - exp_chunk_shift
);
457 red_table_idx_0
^= T
;
459 red_table_idx_0
&= table_idx_mask
;
462 red_table_idx_1
= expz
[exp_chunk_no
+ 1 * (exp_digits
+ 1)];
463 T
= expz
[exp_chunk_no
+ 1 + 1 * (exp_digits
+ 1)];
465 red_table_idx_1
>>= exp_chunk_shift
;
467 * Get additional bits from then next quadword
468 * when 64-bit boundaries are crossed.
470 if (exp_chunk_shift
> 64 - exp_win_size
) {
471 T
<<= (64 - exp_chunk_shift
);
472 red_table_idx_1
^= T
;
474 red_table_idx_1
&= table_idx_mask
;
477 extract(&red_X
[0 * red_digits
], (const BN_ULONG
*)red_table
, (int)red_table_idx_0
, (int)red_table_idx_1
);
480 /* Series of squaring */
481 DAMS((BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_Y
, m
, k0
);
482 DAMS((BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_Y
, m
, k0
);
483 DAMS((BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_Y
, m
, k0
);
484 DAMS((BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_Y
, m
, k0
);
485 DAMS((BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_Y
, m
, k0
);
487 damm((BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_X
, m
, k0
);
493 * NB: After the last AMM of exponentiation in Montgomery domain, the result
494 * may be (modulus_bitsize + 1), but the conversion out of Montgomery domain
495 * performs an AMM(x,1) which guarantees that the final result is less than
496 * |m|, so no conditional subtraction is needed here. See [1] for details.
498 * [1] Gueron, S. Efficient software implementations of modular exponentiation.
499 * DOI: 10.1007/s13389-012-0031-5
502 /* Convert result back in regular 2^52 domain */
503 memset(red_X
, 0, 2 * red_digits
* sizeof(BN_ULONG
));
504 red_X
[0 * red_digits
] = 1;
505 red_X
[1 * red_digits
] = 1;
506 damm(out
, (const BN_ULONG
*)red_Y
, (const BN_ULONG
*)red_X
, m
, k0
);
511 if (storage
!= NULL
) {
512 /* Clear whole storage */
513 OPENSSL_cleanse(storage
, storage_len_bytes
);
514 OPENSSL_free(storage
);
521 static ossl_inline
uint64_t get_digit(const uint8_t *in
, int in_len
)
528 for (; in_len
> 0; in_len
--) {
530 digit
+= (uint64_t)(in
[in_len
- 1]);
536 * Convert array of words in regular (base=2^64) representation to array of
537 * words in redundant (base=2^52) one.
539 static void to_words52(BN_ULONG
*out
, int out_len
,
540 const BN_ULONG
*in
, int in_bitsize
)
542 uint8_t *in_str
= NULL
;
546 /* Check destination buffer capacity */
547 assert(out_len
>= number_of_digits(in_bitsize
, DIGIT_SIZE
));
549 in_str
= (uint8_t *)in
;
551 for (; in_bitsize
>= (2 * DIGIT_SIZE
); in_bitsize
-= (2 * DIGIT_SIZE
), out
+= 2) {
554 memcpy(&digit
, in_str
, sizeof(digit
));
555 out
[0] = digit
& DIGIT_MASK
;
557 memcpy(&digit
, in_str
, sizeof(digit
));
558 out
[1] = (digit
>> 4) & DIGIT_MASK
;
563 if (in_bitsize
> DIGIT_SIZE
) {
564 uint64_t digit
= get_digit(in_str
, 7);
566 out
[0] = digit
& DIGIT_MASK
;
568 in_bitsize
-= DIGIT_SIZE
;
569 digit
= get_digit(in_str
, BITS2WORD8_SIZE(in_bitsize
));
573 } else if (in_bitsize
> 0) {
574 out
[0] = get_digit(in_str
, BITS2WORD8_SIZE(in_bitsize
));
579 while (out_len
> 0) {
586 static ossl_inline
void put_digit(uint8_t *out
, int out_len
, uint64_t digit
)
589 assert(out_len
<= 8);
591 for (; out_len
> 0; out_len
--) {
592 *out
++ = (uint8_t)(digit
& 0xFF);
598 * Convert array of words in redundant (base=2^52) representation to array of
599 * words in regular (base=2^64) one.
601 static void from_words52(BN_ULONG
*out
, int out_bitsize
, const BN_ULONG
*in
)
604 int out_len
= BITS2WORD64_SIZE(out_bitsize
);
609 for (i
= 0; i
< out_len
; i
++)
613 uint8_t *out_str
= (uint8_t *)out
;
615 for (; out_bitsize
>= (2 * DIGIT_SIZE
);
616 out_bitsize
-= (2 * DIGIT_SIZE
), in
+= 2) {
620 memcpy(out_str
, &digit
, sizeof(digit
));
622 digit
= digit
>> 48 | in
[1] << 4;
623 memcpy(out_str
, &digit
, sizeof(digit
));
627 if (out_bitsize
> DIGIT_SIZE
) {
628 put_digit(out_str
, 7, in
[0]);
630 out_bitsize
-= DIGIT_SIZE
;
631 put_digit(out_str
, BITS2WORD8_SIZE(out_bitsize
),
632 (in
[1] << 4 | in
[0] >> 48));
633 } else if (out_bitsize
) {
634 put_digit(out_str
, BITS2WORD8_SIZE(out_bitsize
), in
[0]);
640 * Set bit at index |idx| in the words array |a|.
641 * It does not do any boundaries checks, make sure the index is valid before
642 * calling the function.
644 static ossl_inline
void set_bit(BN_ULONG
*a
, int idx
)
653 a
[i
] |= (((BN_ULONG
)1) << j
);