xref: /openssl/crypto/bn/rsaz_exp_x2.c (revision 7ed6de99)
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     ((unsigned char *)(ptr) + (boundary - (((size_t)(ptr)) & (boundary - 1))))
29 
30 /* Internal radix */
31 # define DIGIT_SIZE (52)
32 /* 52-bit mask */
33 # define DIGIT_MASK ((uint64_t)0xFFFFFFFFFFFFF)
34 
35 # define BITS2WORD8_SIZE(x)  (((x) + 7) >> 3)
36 # 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     (((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 {
52     return (bitsize + digit_size - 1) / digit_size;
53 }
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 {
145     typedef void (*AMM)(BN_ULONG *res, const BN_ULONG *a,
146                         const BN_ULONG *b, const BN_ULONG *m, BN_ULONG k0);
147     int ret = 0;
148 
149     /*
150      * Number of word-size (BN_ULONG) digits to store exponent in redundant
151      * representation.
152      */
153     int exp_digits = number_of_digits(factor_size + 2, DIGIT_SIZE);
154     int coeff_pow = 4 * (DIGIT_SIZE * exp_digits - factor_size);
155 
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;
160 
161     BN_ULONG *base1_red, *m1_red, *rr1_red;
162     BN_ULONG *base2_red, *m2_red, *rr2_red;
163     BN_ULONG *coeff_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 */;
168 
169     const BN_ULONG *exp[2] = {0};
170     BN_ULONG k0[2] = {0};
171     /* AMM = Almost Montgomery Multiplication */
172     AMM amm = NULL;
173 
174     switch (factor_size) {
175     case 1024:
176         amm = ossl_rsaz_amm52x20_x1_ifma256;
177         break;
178     case 1536:
179         amm = ossl_rsaz_amm52x30_x1_ifma256;
180         break;
181     case 2048:
182         amm = ossl_rsaz_amm52x40_x1_ifma256;
183         break;
184     default:
185         goto err;
186     }
187 
188     storage = (BN_ULONG *)OPENSSL_malloc(storage_len_bytes);
189     if (storage == NULL)
190         goto err;
191     storage_aligned = (BN_ULONG *)ALIGN_OF(storage, 64);
192 
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;
201 
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);
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     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);
229 
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 */
232 
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 */
235 
236     exp[0] = exp1;
237     exp[1] = exp2;
238 
239     k0[0] = k0_1;
240     k0[1] = k0_2;
241 
242     /* Dual (2-exps in parallel) exponentiation */
243     ret = RSAZ_mod_exp_x2_ifma256(rr1_red, base1_red, exp, m1_red, rr1_red,
244                                   k0, factor_size);
245     if (!ret)
246         goto err;
247 
248     /* Convert rr_i back to regular radix */
249     from_words52(res1, factor_size, rr1_red);
250     from_words52(res2, factor_size, rr2_red);
251 
252     /* bn_reduce_once_in_place expects number of BN_ULONG, not bit size */
253     factor_size /= sizeof(BN_ULONG) * 8;
254 
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);
257 
258 err:
259     if (storage != NULL) {
260         OPENSSL_cleanse(storage, storage_len_bytes);
261         OPENSSL_free(storage);
262     }
263     return ret;
264 }
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 {
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);
301 
302     int ret = 0;
303     int idx;
304 
305     /* Exponent window size */
306     int exp_win_size = 5;
307     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     int red_digits = 0;
314     int exp_digits = 0;
315 
316     BN_ULONG *storage = NULL;
317     BN_ULONG *storage_aligned = NULL;
318     int storage_len_bytes = 0;
319 
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] */
327 
328     /* Dual AMM */
329     DAMM damm = NULL;
330     /* Extractor from red_table */
331     DEXTRACT extract = NULL;
332 
333 /*
334  * Squaring is done using multiplication now. That can be a subject of
335  * optimization in future.
336  */
337 # define DAMS(r,a,m,k0) damm((r),(a),(a),(m),(k0))
338 
339     switch (modulus_bitsize) {
340     case 1024:
341         red_digits = 20;
342         exp_digits = 16;
343         damm = ossl_rsaz_amm52x20_x2_ifma256;
344         extract = ossl_extract_multiplier_2x20_win5;
345         break;
346     case 1536:
347         /* Extended with 2 digits padding to avoid mask ops in high YMM register */
348         red_digits = 30 + 2;
349         exp_digits = 24;
350         damm = ossl_rsaz_amm52x30_x2_ifma256;
351         extract = ossl_extract_multiplier_2x30_win5;
352         break;
353     case 2048:
354         red_digits = 40;
355         exp_digits = 32;
356         damm = ossl_rsaz_amm52x40_x2_ifma256;
357         extract = ossl_extract_multiplier_2x40_win5;
358         break;
359     default:
360         goto err;
361     }
362 
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      */
367                        * sizeof(BN_ULONG)
368                        + 64;                                    /* alignment */
369 
370     storage = (BN_ULONG *)OPENSSL_zalloc(storage_len_bytes);
371     if (storage == NULL)
372         goto err;
373     storage_aligned = (BN_ULONG *)ALIGN_OF(storage, 64);
374 
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);
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     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);
389 
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);
396     }
397 
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;
403 
404     /* Exponentiation */
405     {
406         const int rem = modulus_bitsize % exp_win_size;
407         const BN_ULONG table_idx_mask = exp_win_mask;
408 
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;
412 
413         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         OPENSSL_assert(rem != 0);
424 
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)];
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         red_table_idx_0 >>= exp_chunk_shift;
434         red_table_idx_1 >>= exp_chunk_shift;
435 
436         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         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             {
442                 BN_ULONG T;
443 
444                 exp_chunk_no = exp_bit_no / 64;
445                 exp_chunk_shift = exp_bit_no % 64;
446                 {
447                     red_table_idx_0 = expz[exp_chunk_no + 0 * (exp_digits + 1)];
448                     T = expz[exp_chunk_no + 1 + 0 * (exp_digits + 1)];
449 
450                     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                     if (exp_chunk_shift > 64 - exp_win_size) {
456                         T <<= (64 - exp_chunk_shift);
457                         red_table_idx_0 ^= T;
458                     }
459                     red_table_idx_0 &= table_idx_mask;
460                 }
461                 {
462                     red_table_idx_1 = expz[exp_chunk_no + 1 * (exp_digits + 1)];
463                     T = expz[exp_chunk_no + 1 + 1 * (exp_digits + 1)];
464 
465                     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                     if (exp_chunk_shift > 64 - exp_win_size) {
471                         T <<= (64 - exp_chunk_shift);
472                         red_table_idx_1 ^= T;
473                     }
474                     red_table_idx_1 &= table_idx_mask;
475                 }
476 
477                 extract(&red_X[0 * red_digits], (const BN_ULONG*)red_table, (int)red_table_idx_0, (int)red_table_idx_1);
478             }
479 
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);
486 
487             damm((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0);
488         }
489     }
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     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);
507 
508     ret = 1;
509 
510 err:
511     if (storage != NULL) {
512         /* Clear whole storage */
513         OPENSSL_cleanse(storage, storage_len_bytes);
514         OPENSSL_free(storage);
515     }
516 
517 #undef DAMS
518     return ret;
519 }
520 
521 static ossl_inline uint64_t get_digit(const uint8_t *in, int in_len)
522 {
523     uint64_t digit = 0;
524 
525     assert(in != NULL);
526     assert(in_len <= 8);
527 
528     for (; in_len > 0; in_len--) {
529         digit <<= 8;
530         digit += (uint64_t)(in[in_len - 1]);
531     }
532     return digit;
533 }
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 {
542     uint8_t *in_str = NULL;
543 
544     assert(out != NULL);
545     assert(in != NULL);
546     /* Check destination buffer capacity */
547     assert(out_len >= number_of_digits(in_bitsize, DIGIT_SIZE));
548 
549     in_str = (uint8_t *)in;
550 
551     for (; in_bitsize >= (2 * DIGIT_SIZE); in_bitsize -= (2 * DIGIT_SIZE), out += 2) {
552         uint64_t digit;
553 
554         memcpy(&digit, in_str, sizeof(digit));
555         out[0] = digit & DIGIT_MASK;
556         in_str += 6;
557         memcpy(&digit, in_str, sizeof(digit));
558         out[1] = (digit >> 4) & DIGIT_MASK;
559         in_str += 7;
560         out_len -= 2;
561     }
562 
563     if (in_bitsize > DIGIT_SIZE) {
564         uint64_t digit = get_digit(in_str, 7);
565 
566         out[0] = digit & DIGIT_MASK;
567         in_str += 6;
568         in_bitsize -= DIGIT_SIZE;
569         digit = get_digit(in_str, BITS2WORD8_SIZE(in_bitsize));
570         out[1] = digit >> 4;
571         out += 2;
572         out_len -= 2;
573     } else if (in_bitsize > 0) {
574         out[0] = get_digit(in_str, BITS2WORD8_SIZE(in_bitsize));
575         out++;
576         out_len--;
577     }
578 
579     memset(out, 0, out_len * sizeof(BN_ULONG));
580 }
581 
582 static ossl_inline void put_digit(uint8_t *out, int out_len, uint64_t digit)
583 {
584     assert(out != NULL);
585     assert(out_len <= 8);
586 
587     for (; out_len > 0; out_len--) {
588         *out++ = (uint8_t)(digit & 0xFF);
589         digit >>= 8;
590     }
591 }
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 {
599     int i;
600     int out_len = BITS2WORD64_SIZE(out_bitsize);
601 
602     assert(out != NULL);
603     assert(in != NULL);
604 
605     for (i = 0; i < out_len; i++)
606         out[i] = 0;
607 
608     {
609         uint8_t *out_str = (uint8_t *)out;
610 
611         for (; out_bitsize >= (2 * DIGIT_SIZE);
612                out_bitsize -= (2 * DIGIT_SIZE), in += 2) {
613             uint64_t digit;
614 
615             digit = in[0];
616             memcpy(out_str, &digit, sizeof(digit));
617             out_str += 6;
618             digit = digit >> 48 | in[1] << 4;
619             memcpy(out_str, &digit, sizeof(digit));
620             out_str += 7;
621         }
622 
623         if (out_bitsize > DIGIT_SIZE) {
624             put_digit(out_str, 7, in[0]);
625             out_str += 6;
626             out_bitsize -= DIGIT_SIZE;
627             put_digit(out_str, BITS2WORD8_SIZE(out_bitsize),
628                         (in[1] << 4 | in[0] >> 48));
629         } else if (out_bitsize) {
630             put_digit(out_str, BITS2WORD8_SIZE(out_bitsize), in[0]);
631         }
632     }
633 }
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 {
642     assert(a != NULL);
643 
644     {
645         int i, j;
646 
647         i = idx / BN_BITS2;
648         j = idx % BN_BITS2;
649         a[i] |= (((BN_ULONG)1) << j);
650     }
651 }
652 
653 #endif
654