2
0

rsaz_exp_x2.c 23 KB

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