rsaz_exp_x2.c@ 100947

Last change on this file since 100947 was 97372, checked in by vboxsync, 2 years ago
libs: Switch to openssl-3.0.7, bugref:10317
File size: 19.2 KB

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

Note: See TracBrowser for help on using the repository browser.

source: vbox/trunk/src/libs/openssl-3.0.7/crypto/bn/rsaz_exp_x2.c@ 100947

Download in other formats: