crypto/bn/rsaz_exp_x2.c

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