crypto/bn/rsaz_exp_x2.c

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