1/* 2 * Copyright 2020-2023 The OpenSSL Project Authors. All Rights Reserved. 3 * Copyright (c) 2020, Intel Corporation. All Rights Reserved. 4 * 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 9 * 10 * 11 * Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov 12 * Intel Corporation 13 * 14 */ 15 16#include <openssl/opensslconf.h> 17#include <openssl/crypto.h> 18#include "rsaz_exp.h" 19 20#ifndef RSAZ_ENABLED 21NON_EMPTY_TRANSLATION_UNIT 22#else 23# include <assert.h> 24# include <string.h> 25 26# if defined(__GNUC__) 27# define ALIGN64 __attribute__((aligned(64))) 28# elif defined(_MSC_VER) 29# define ALIGN64 __declspec(align(64)) 30# else 31# define ALIGN64 32# endif 33 34# define ALIGN_OF(ptr, boundary) \ 35 ((unsigned char *)(ptr) + (boundary - (((size_t)(ptr)) & (boundary - 1)))) 36 37/* Internal radix */ 38# define DIGIT_SIZE (52) 39/* 52-bit mask */ 40# define DIGIT_MASK ((uint64_t)0xFFFFFFFFFFFFF) 41 42# define BITS2WORD8_SIZE(x) (((x) + 7) >> 3) 43# define BITS2WORD64_SIZE(x) (((x) + 63) >> 6) 44 45static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len); 46static ossl_inline void put_digit52(uint8_t *out, int out_len, uint64_t digit); 47static void to_words52(BN_ULONG *out, int out_len, const BN_ULONG *in, 48 int in_bitsize); 49static void from_words52(BN_ULONG *bn_out, int out_bitsize, const BN_ULONG *in); 50static ossl_inline void set_bit(BN_ULONG *a, int idx); 51 52/* Number of |digit_size|-bit digits in |bitsize|-bit value */ 53static ossl_inline int number_of_digits(int bitsize, int digit_size) 54{ 55 return (bitsize + digit_size - 1) / digit_size; 56} 57 58typedef void (*AMM52)(BN_ULONG *res, const BN_ULONG *base, 59 const BN_ULONG *exp, const BN_ULONG *m, BN_ULONG k0); 60typedef void (*EXP52_x2)(BN_ULONG *res, const BN_ULONG *base, 61 const BN_ULONG *exp[2], const BN_ULONG *m, 62 const BN_ULONG *rr, const BN_ULONG k0[2]); 63 64/* 65 * For details of the methods declared below please refer to 66 * crypto/bn/asm/rsaz-avx512.pl 67 * 68 * Naming notes: 69 * amm = Almost Montgomery Multiplication 70 * ams = Almost Montgomery Squaring 71 * 52x20 - data represented as array of 20 digits in 52-bit radix 72 * _x1_/_x2_ - 1 or 2 independent inputs/outputs 73 * _256 suffix - uses 256-bit (AVX512VL) registers 74 */ 75 76/*AMM = Almost Montgomery Multiplication. */ 77void ossl_rsaz_amm52x20_x1_256(BN_ULONG *res, const BN_ULONG *base, 78 const BN_ULONG *exp, const BN_ULONG *m, 79 BN_ULONG k0); 80static void RSAZ_exp52x20_x2_256(BN_ULONG *res, const BN_ULONG *base, 81 const BN_ULONG *exp[2], const BN_ULONG *m, 82 const BN_ULONG *rr, const BN_ULONG k0[2]); 83void ossl_rsaz_amm52x20_x2_256(BN_ULONG *out, const BN_ULONG *a, 84 const BN_ULONG *b, const BN_ULONG *m, 85 const BN_ULONG k0[2]); 86void ossl_extract_multiplier_2x20_win5(BN_ULONG *red_Y, 87 const BN_ULONG *red_table, 88 int red_table_idx, int tbl_idx); 89 90/* 91 * Dual Montgomery modular exponentiation using prime moduli of the 92 * same bit size, optimized with AVX512 ISA. 93 * 94 * Input and output parameters for each exponentiation are independent and 95 * denoted here by index |i|, i = 1..2. 96 * 97 * Input and output are all in regular 2^64 radix. 98 * 99 * Each moduli shall be |factor_size| bit size. 100 * 101 * NOTE: currently only 2x1024 case is supported. 102 * 103 * [out] res|i| - result of modular exponentiation: array of qword values 104 * in regular (2^64) radix. Size of array shall be enough 105 * to hold |factor_size| bits. 106 * [in] base|i| - base 107 * [in] exp|i| - exponent 108 * [in] m|i| - moduli 109 * [in] rr|i| - Montgomery parameter RR = R^2 mod m|i| 110 * [in] k0_|i| - Montgomery parameter k0 = -1/m|i| mod 2^64 111 * [in] factor_size - moduli bit size 112 * 113 * \return 0 in case of failure, 114 * 1 in case of success. 115 */ 116int ossl_rsaz_mod_exp_avx512_x2(BN_ULONG *res1, 117 const BN_ULONG *base1, 118 const BN_ULONG *exp1, 119 const BN_ULONG *m1, 120 const BN_ULONG *rr1, 121 BN_ULONG k0_1, 122 BN_ULONG *res2, 123 const BN_ULONG *base2, 124 const BN_ULONG *exp2, 125 const BN_ULONG *m2, 126 const BN_ULONG *rr2, 127 BN_ULONG k0_2, 128 int factor_size) 129{ 130 int ret = 0; 131 132 /* 133 * Number of word-size (BN_ULONG) digits to store exponent in redundant 134 * representation. 135 */ 136 int exp_digits = number_of_digits(factor_size + 2, DIGIT_SIZE); 137 int coeff_pow = 4 * (DIGIT_SIZE * exp_digits - factor_size); 138 BN_ULONG *base1_red, *m1_red, *rr1_red; 139 BN_ULONG *base2_red, *m2_red, *rr2_red; 140 BN_ULONG *coeff_red; 141 BN_ULONG *storage = NULL; 142 BN_ULONG *storage_aligned = NULL; 143 BN_ULONG storage_len_bytes = 7 * exp_digits * sizeof(BN_ULONG); 144 145 /* AMM = Almost Montgomery Multiplication */ 146 AMM52 amm = NULL; 147 /* Dual (2-exps in parallel) exponentiation */ 148 EXP52_x2 exp_x2 = NULL; 149 150 const BN_ULONG *exp[2] = {0}; 151 BN_ULONG k0[2] = {0}; 152 153 /* Only 1024-bit factor size is supported now */ 154 switch (factor_size) { 155 case 1024: 156 amm = ossl_rsaz_amm52x20_x1_256; 157 exp_x2 = RSAZ_exp52x20_x2_256; 158 break; 159 default: 160 goto err; 161 } 162 163 storage = (BN_ULONG *)OPENSSL_malloc(storage_len_bytes + 64); 164 if (storage == NULL) 165 goto err; 166 storage_aligned = (BN_ULONG *)ALIGN_OF(storage, 64); 167 168 /* Memory layout for red(undant) representations */ 169 base1_red = storage_aligned; 170 base2_red = storage_aligned + 1 * exp_digits; 171 m1_red = storage_aligned + 2 * exp_digits; 172 m2_red = storage_aligned + 3 * exp_digits; 173 rr1_red = storage_aligned + 4 * exp_digits; 174 rr2_red = storage_aligned + 5 * exp_digits; 175 coeff_red = storage_aligned + 6 * exp_digits; 176 177 /* Convert base_i, m_i, rr_i, from regular to 52-bit radix */ 178 to_words52(base1_red, exp_digits, base1, factor_size); 179 to_words52(base2_red, exp_digits, base2, factor_size); 180 to_words52(m1_red, exp_digits, m1, factor_size); 181 to_words52(m2_red, exp_digits, m2, factor_size); 182 to_words52(rr1_red, exp_digits, rr1, factor_size); 183 to_words52(rr2_red, exp_digits, rr2, factor_size); 184 185 /* 186 * Compute target domain Montgomery converters RR' for each modulus 187 * based on precomputed original domain's RR. 188 * 189 * RR -> RR' transformation steps: 190 * (1) coeff = 2^k 191 * (2) t = AMM(RR,RR) = RR^2 / R' mod m 192 * (3) RR' = AMM(t, coeff) = RR^2 * 2^k / R'^2 mod m 193 * where 194 * k = 4 * (52 * digits52 - modlen) 195 * R = 2^(64 * ceil(modlen/64)) mod m 196 * RR = R^2 mod M 197 * R' = 2^(52 * ceil(modlen/52)) mod m 198 * 199 * modlen = 1024: k = 64, RR = 2^2048 mod m, RR' = 2^2080 mod m 200 */ 201 memset(coeff_red, 0, exp_digits * sizeof(BN_ULONG)); 202 /* (1) in reduced domain representation */ 203 set_bit(coeff_red, 64 * (int)(coeff_pow / 52) + coeff_pow % 52); 204 205 amm(rr1_red, rr1_red, rr1_red, m1_red, k0_1); /* (2) for m1 */ 206 amm(rr1_red, rr1_red, coeff_red, m1_red, k0_1); /* (3) for m1 */ 207 208 amm(rr2_red, rr2_red, rr2_red, m2_red, k0_2); /* (2) for m2 */ 209 amm(rr2_red, rr2_red, coeff_red, m2_red, k0_2); /* (3) for m2 */ 210 211 exp[0] = exp1; 212 exp[1] = exp2; 213 214 k0[0] = k0_1; 215 k0[1] = k0_2; 216 217 exp_x2(rr1_red, base1_red, exp, m1_red, rr1_red, k0); 218 219 /* Convert rr_i back to regular radix */ 220 from_words52(res1, factor_size, rr1_red); 221 from_words52(res2, factor_size, rr2_red); 222 223 /* bn_reduce_once_in_place expects number of BN_ULONG, not bit size */ 224 factor_size /= sizeof(BN_ULONG) * 8; 225 226 bn_reduce_once_in_place(res1, /*carry=*/0, m1, storage, factor_size); 227 bn_reduce_once_in_place(res2, /*carry=*/0, m2, storage, factor_size); 228 229 ret = 1; 230err: 231 if (storage != NULL) { 232 OPENSSL_cleanse(storage, storage_len_bytes); 233 OPENSSL_free(storage); 234 } 235 return ret; 236} 237 238/* 239 * Dual 1024-bit w-ary modular exponentiation using prime moduli of the same 240 * bit size using Almost Montgomery Multiplication, optimized with AVX512_IFMA 241 * ISA. 242 * 243 * The parameter w (window size) = 5. 244 * 245 * [out] res - result of modular exponentiation: 2x20 qword 246 * values in 2^52 radix. 247 * [in] base - base (2x20 qword values in 2^52 radix) 248 * [in] exp - array of 2 pointers to 16 qword values in 2^64 radix. 249 * Exponent is not converted to redundant representation. 250 * [in] m - moduli (2x20 qword values in 2^52 radix) 251 * [in] rr - Montgomery parameter for 2 moduli: RR = 2^2080 mod m. 252 * (2x20 qword values in 2^52 radix) 253 * [in] k0 - Montgomery parameter for 2 moduli: k0 = -1/m mod 2^64 254 * 255 * \return (void). 256 */ 257static void RSAZ_exp52x20_x2_256(BN_ULONG *out, /* [2][20] */ 258 const BN_ULONG *base, /* [2][20] */ 259 const BN_ULONG *exp[2], /* 2x16 */ 260 const BN_ULONG *m, /* [2][20] */ 261 const BN_ULONG *rr, /* [2][20] */ 262 const BN_ULONG k0[2]) 263{ 264# define BITSIZE_MODULUS (1024) 265# define EXP_WIN_SIZE (5) 266# define EXP_WIN_MASK ((1U << EXP_WIN_SIZE) - 1) 267/* 268 * Number of digits (64-bit words) in redundant representation to handle 269 * modulus bits 270 */ 271# define RED_DIGITS (20) 272# define EXP_DIGITS (16) 273# define DAMM ossl_rsaz_amm52x20_x2_256 274/* 275 * Squaring is done using multiplication now. That can be a subject of 276 * optimization in future. 277 */ 278# define DAMS(r,a,m,k0) \ 279 ossl_rsaz_amm52x20_x2_256((r),(a),(a),(m),(k0)) 280 281 /* Allocate stack for red(undant) result Y and multiplier X */ 282 ALIGN64 BN_ULONG red_Y[2][RED_DIGITS]; 283 ALIGN64 BN_ULONG red_X[2][RED_DIGITS]; 284 285 /* Allocate expanded exponent */ 286 ALIGN64 BN_ULONG expz[2][EXP_DIGITS + 1]; 287 288 /* Pre-computed table of base powers */ 289 ALIGN64 BN_ULONG red_table[1U << EXP_WIN_SIZE][2][RED_DIGITS]; 290 291 int idx; 292 293 memset(red_Y, 0, sizeof(red_Y)); 294 memset(red_table, 0, sizeof(red_table)); 295 memset(red_X, 0, sizeof(red_X)); 296 297 /* 298 * Compute table of powers base^i, i = 0, ..., (2^EXP_WIN_SIZE) - 1 299 * table[0] = mont(x^0) = mont(1) 300 * table[1] = mont(x^1) = mont(x) 301 */ 302 red_X[0][0] = 1; 303 red_X[1][0] = 1; 304 DAMM(red_table[0][0], (const BN_ULONG*)red_X, rr, m, k0); 305 DAMM(red_table[1][0], base, rr, m, k0); 306 307 for (idx = 1; idx < (int)((1U << EXP_WIN_SIZE) / 2); idx++) { 308 DAMS(red_table[2 * idx + 0][0], red_table[1 * idx][0], m, k0); 309 DAMM(red_table[2 * idx + 1][0], red_table[2 * idx][0], red_table[1][0], m, k0); 310 } 311 312 /* Copy and expand exponents */ 313 memcpy(expz[0], exp[0], EXP_DIGITS * sizeof(BN_ULONG)); 314 expz[0][EXP_DIGITS] = 0; 315 memcpy(expz[1], exp[1], EXP_DIGITS * sizeof(BN_ULONG)); 316 expz[1][EXP_DIGITS] = 0; 317 318 /* Exponentiation */ 319 { 320 const int rem = BITSIZE_MODULUS % EXP_WIN_SIZE; 321 BN_ULONG table_idx_mask = EXP_WIN_MASK; 322 323 int exp_bit_no = BITSIZE_MODULUS - rem; 324 int exp_chunk_no = exp_bit_no / 64; 325 int exp_chunk_shift = exp_bit_no % 64; 326 327 BN_ULONG red_table_idx_0, red_table_idx_1; 328 329 /* 330 * If rem == 0, then 331 * exp_bit_no = modulus_bitsize - exp_win_size 332 * However, this isn't possible because rem is { 1024, 1536, 2048 } % 5 333 * which is { 4, 1, 3 } respectively. 334 * 335 * If this assertion ever fails the fix above is easy. 336 */ 337 OPENSSL_assert(rem != 0); 338 339 /* Process 1-st exp window - just init result */ 340 red_table_idx_0 = expz[0][exp_chunk_no]; 341 red_table_idx_1 = expz[1][exp_chunk_no]; 342 /* 343 * The function operates with fixed moduli sizes divisible by 64, 344 * thus table index here is always in supported range [0, EXP_WIN_SIZE). 345 */ 346 red_table_idx_0 >>= exp_chunk_shift; 347 red_table_idx_1 >>= exp_chunk_shift; 348 349 ossl_extract_multiplier_2x20_win5(red_Y[0], (const BN_ULONG*)red_table, 350 (int)red_table_idx_0, 0); 351 ossl_extract_multiplier_2x20_win5(red_Y[1], (const BN_ULONG*)red_table, 352 (int)red_table_idx_1, 1); 353 354 /* Process other exp windows */ 355 for (exp_bit_no -= EXP_WIN_SIZE; exp_bit_no >= 0; exp_bit_no -= EXP_WIN_SIZE) { 356 /* Extract pre-computed multiplier from the table */ 357 { 358 BN_ULONG T; 359 360 exp_chunk_no = exp_bit_no / 64; 361 exp_chunk_shift = exp_bit_no % 64; 362 { 363 red_table_idx_0 = expz[0][exp_chunk_no]; 364 T = expz[0][exp_chunk_no + 1]; 365 366 red_table_idx_0 >>= exp_chunk_shift; 367 /* 368 * Get additional bits from then next quadword 369 * when 64-bit boundaries are crossed. 370 */ 371 if (exp_chunk_shift > 64 - EXP_WIN_SIZE) { 372 T <<= (64 - exp_chunk_shift); 373 red_table_idx_0 ^= T; 374 } 375 red_table_idx_0 &= table_idx_mask; 376 377 ossl_extract_multiplier_2x20_win5(red_X[0], 378 (const BN_ULONG*)red_table, 379 (int)red_table_idx_0, 0); 380 } 381 { 382 red_table_idx_1 = expz[1][exp_chunk_no]; 383 T = expz[1][exp_chunk_no + 1]; 384 385 red_table_idx_1 >>= exp_chunk_shift; 386 /* 387 * Get additional bits from then next quadword 388 * when 64-bit boundaries are crossed. 389 */ 390 if (exp_chunk_shift > 64 - EXP_WIN_SIZE) { 391 T <<= (64 - exp_chunk_shift); 392 red_table_idx_1 ^= T; 393 } 394 red_table_idx_1 &= table_idx_mask; 395 396 ossl_extract_multiplier_2x20_win5(red_X[1], 397 (const BN_ULONG*)red_table, 398 (int)red_table_idx_1, 1); 399 } 400 } 401 402 /* Series of squaring */ 403 DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0); 404 DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0); 405 DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0); 406 DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0); 407 DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0); 408 409 DAMM((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0); 410 } 411 } 412 413 /* 414 * 415 * NB: After the last AMM of exponentiation in Montgomery domain, the result 416 * may be 1025-bit, but the conversion out of Montgomery domain performs an 417 * AMM(x,1) which guarantees that the final result is less than |m|, so no 418 * conditional subtraction is needed here. See "Efficient Software 419 * Implementations of Modular Exponentiation" (by Shay Gueron) paper for details. 420 */ 421 422 /* Convert result back in regular 2^52 domain */ 423 memset(red_X, 0, sizeof(red_X)); 424 red_X[0][0] = 1; 425 red_X[1][0] = 1; 426 DAMM(out, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0); 427 428 /* Clear exponents */ 429 OPENSSL_cleanse(expz, sizeof(expz)); 430 OPENSSL_cleanse(red_Y, sizeof(red_Y)); 431 432# undef DAMS 433# undef DAMM 434# undef EXP_DIGITS 435# undef RED_DIGITS 436# undef EXP_WIN_MASK 437# undef EXP_WIN_SIZE 438# undef BITSIZE_MODULUS 439} 440 441static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len) 442{ 443 uint64_t digit = 0; 444 445 assert(in != NULL); 446 447 for (; in_len > 0; in_len--) { 448 digit <<= 8; 449 digit += (uint64_t)(in[in_len - 1]); 450 } 451 return digit; 452} 453 454/* 455 * Convert array of words in regular (base=2^64) representation to array of 456 * words in redundant (base=2^52) one. 457 */ 458static void to_words52(BN_ULONG *out, int out_len, 459 const BN_ULONG *in, int in_bitsize) 460{ 461 uint8_t *in_str = NULL; 462 463 assert(out != NULL); 464 assert(in != NULL); 465 /* Check destination buffer capacity */ 466 assert(out_len >= number_of_digits(in_bitsize, DIGIT_SIZE)); 467 468 in_str = (uint8_t *)in; 469 470 for (; in_bitsize >= (2 * DIGIT_SIZE); in_bitsize -= (2 * DIGIT_SIZE), out += 2) { 471 uint64_t digit; 472 473 memcpy(&digit, in_str, sizeof(digit)); 474 out[0] = digit & DIGIT_MASK; 475 in_str += 6; 476 memcpy(&digit, in_str, sizeof(digit)); 477 out[1] = (digit >> 4) & DIGIT_MASK; 478 in_str += 7; 479 out_len -= 2; 480 } 481 482 if (in_bitsize > DIGIT_SIZE) { 483 uint64_t digit = get_digit52(in_str, 7); 484 485 out[0] = digit & DIGIT_MASK; 486 in_str += 6; 487 in_bitsize -= DIGIT_SIZE; 488 digit = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize)); 489 out[1] = digit >> 4; 490 out += 2; 491 out_len -= 2; 492 } else if (in_bitsize > 0) { 493 out[0] = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize)); 494 out++; 495 out_len--; 496 } 497 498 while (out_len > 0) { 499 *out = 0; 500 out_len--; 501 out++; 502 } 503} 504 505static ossl_inline void put_digit52(uint8_t *pStr, int strLen, uint64_t digit) 506{ 507 assert(pStr != NULL); 508 509 for (; strLen > 0; strLen--) { 510 *pStr++ = (uint8_t)(digit & 0xFF); 511 digit >>= 8; 512 } 513} 514 515/* 516 * Convert array of words in redundant (base=2^52) representation to array of 517 * words in regular (base=2^64) one. 518 */ 519static void from_words52(BN_ULONG *out, int out_bitsize, const BN_ULONG *in) 520{ 521 int i; 522 int out_len = BITS2WORD64_SIZE(out_bitsize); 523 524 assert(out != NULL); 525 assert(in != NULL); 526 527 for (i = 0; i < out_len; i++) 528 out[i] = 0; 529 530 { 531 uint8_t *out_str = (uint8_t *)out; 532 533 for (; out_bitsize >= (2 * DIGIT_SIZE); 534 out_bitsize -= (2 * DIGIT_SIZE), in += 2) { 535 uint64_t digit; 536 537 digit = in[0]; 538 memcpy(out_str, &digit, sizeof(digit)); 539 out_str += 6; 540 digit = digit >> 48 | in[1] << 4; 541 memcpy(out_str, &digit, sizeof(digit)); 542 out_str += 7; 543 } 544 545 if (out_bitsize > DIGIT_SIZE) { 546 put_digit52(out_str, 7, in[0]); 547 out_str += 6; 548 out_bitsize -= DIGIT_SIZE; 549 put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize), 550 (in[1] << 4 | in[0] >> 48)); 551 } else if (out_bitsize) { 552 put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize), in[0]); 553 } 554 } 555} 556 557/* 558 * Set bit at index |idx| in the words array |a|. 559 * It does not do any boundaries checks, make sure the index is valid before 560 * calling the function. 561 */ 562static ossl_inline void set_bit(BN_ULONG *a, int idx) 563{ 564 assert(a != NULL); 565 566 { 567 int i, j; 568 569 i = idx / BN_BITS2; 570 j = idx % BN_BITS2; 571 a[i] |= (((BN_ULONG)1) << j); 572 } 573} 574 575#endif 576