Coverage Report

Created: 2024-11-21 06:47

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