xref: /freebsd/crypto/libecc/src/nn/nn_mod_pow.c (revision f0865ec9906d5a18fa2a3b61381f22ce16e606ad)
1 /*
2  *  Copyright (C) 2021 - This file is part of libecc project
3  *
4  *  Authors:
5  *      Ryad BENADJILA <ryadbenadjila@gmail.com>
6  *      Arnaud EBALARD <arnaud.ebalard@ssi.gouv.fr>
7  *      Jean-Pierre FLORI <jean-pierre.flori@ssi.gouv.fr>
8  *
9  *  Contributors:
10  *      Nicolas VIVET <nicolas.vivet@ssi.gouv.fr>
11  *      Karim KHALFALLAH <karim.khalfallah@ssi.gouv.fr>
12  *
13  *  This software is licensed under a dual BSD and GPL v2 license.
14  *  See LICENSE file at the root folder of the project.
15  */
16 #include <libecc/nn/nn_mul_redc1.h>
17 #include <libecc/nn/nn_div_public.h>
18 #include <libecc/nn/nn_logical.h>
19 #include <libecc/nn/nn_mod_pow.h>
20 #include <libecc/nn/nn_rand.h>
21 #include <libecc/nn/nn.h>
22 
23 /*
24  * NOT constant time with regard to the bitlength of exp.
25  *
26  * Implements a left to right Montgomery Ladder for modular exponentiation.
27  * This is an internal common helper and assumes that all the initialization
28  * and aliasing of inputs/outputs are handled by the callers. Depending on
29  * the inputs, redcification is optionally used.
30  * The base is reduced if necessary.
31  *
32  * Montgomery Ladder is masked using Itoh et al. anti-ADPA
33  * (Address-bit DPA) countermeasure.
34  * See "A Practical Countermeasure against Address-Bit Differential Power Analysis"
35  * by Itoh, Izu and Takenaka for more information.
36 
37  * Returns 0 on success, -1 on error.
38  */
_nn_exp_monty_ladder_ltr(nn_t out,nn_src_t base,nn_src_t exp,nn_src_t mod,nn_src_t r,nn_src_t r_square,word_t mpinv)39 ATTRIBUTE_WARN_UNUSED_RET static int _nn_exp_monty_ladder_ltr(nn_t out, nn_src_t base, nn_src_t exp, nn_src_t mod, nn_src_t r, nn_src_t r_square, word_t mpinv)
40 {
41 	nn T[3];
42 	nn mask;
43 	bitcnt_t explen, oldexplen;
44  	u8 expbit, rbit;
45 	int ret, cmp;
46 	T[0].magic = T[1].magic = T[2].magic = mask.magic = WORD(0);
47 
48 	/* Initialize out */
49 	ret = nn_init(out, 0); EG(ret, err);
50 
51 	ret = nn_init(&T[0], 0); EG(ret, err);
52 	ret = nn_init(&T[1], 0); EG(ret, err);
53 	ret = nn_init(&T[2], 0); EG(ret, err);
54 
55 	/* Generate our Itoh random mask */
56 	ret = nn_get_random_len(&mask, NN_MAX_BYTE_LEN); EG(ret, err);
57 
58 	ret = nn_bitlen(exp, &explen); EG(ret, err);
59 
60 
61 	/* From now on, since we deal with Itoh's countermeasure, we must have at
62 	 * least 2 bits in the exponent. We will deal with the particular cases of 0 and 1
63 	 * bit exponents later.
64 	 */
65 	oldexplen = explen;
66 	explen = (explen < 2) ? 2 : explen;
67 
68 	ret = nn_getbit(&mask, (bitcnt_t)(explen - 1), &rbit); EG(ret, err);
69 
70 	/* Reduce the base if necessary */
71 	ret = nn_cmp(base, mod, &cmp); EG(ret, err);
72 	if(cmp >= 0){
73 		/* Modular reduction */
74 		ret = nn_mod(&T[rbit], base, mod); EG(ret, err);
75 		if(r != NULL){
76 			/* Redcify the base if necessary */
77 			ret = nn_mul_redc1(&T[rbit], &T[rbit], r_square, mod, mpinv); EG(ret, err);
78 		}
79 	}
80 	else{
81 		if(r != NULL){
82 			/* Redcify the base if necessary */
83 			ret = nn_mul_redc1(&T[rbit], base, r_square, mod, mpinv); EG(ret, err);
84 		}
85 		else{
86 			ret = nn_copy(&T[rbit], base); EG(ret, err);
87 		}
88 	}
89 
90 	/* We implement the Montgomery ladder exponentiation with Itoh masking using three
91 	 * registers T[0], T[1] and T[2]. The random mask is in 'mask'.
92 	 */
93 	if(r != NULL){
94 		ret = nn_mul_redc1(&T[1-rbit], &T[rbit], &T[rbit], mod, mpinv); EG(ret, err);
95 	}
96 	else{
97 		ret = nn_mod_mul(&T[1-rbit], &T[rbit], &T[rbit], mod); EG(ret, err);
98 	}
99 
100 	/* Now proceed with the Montgomery Ladder algorithm.
101 	 */
102 	explen = (bitcnt_t)(explen - 1);
103 	while (explen > 0) {
104 		u8 rbit_next;
105 		explen = (bitcnt_t)(explen - 1);
106 
107 		/* rbit is r[i+1], and rbit_next is r[i] */
108 		ret = nn_getbit(&mask, explen, &rbit_next); EG(ret, err);
109 		/* Get the exponent bit */
110 		ret = nn_getbit(exp, explen, &expbit); EG(ret, err);
111 
112 		/* Square */
113 		if(r != NULL){
114 			ret = nn_mul_redc1(&T[2], &T[expbit ^ rbit], &T[expbit ^ rbit], mod, mpinv); EG(ret, err);
115 		}
116 		else{
117 			ret = nn_mod_mul(&T[2], &T[expbit ^ rbit], &T[expbit ^ rbit], mod); EG(ret, err);
118 		}
119 		/* Multiply */
120 		if(r != NULL){
121 			ret = nn_mul_redc1(&T[1], &T[0], &T[1], mod, mpinv); EG(ret, err);
122 		}
123 		else{
124 			ret = nn_mod_mul(&T[1], &T[0], &T[1], mod); EG(ret, err);
125 		}
126 		/* Copy */
127 		ret = nn_copy(&T[0], &T[2 - (expbit ^ rbit_next)]); EG(ret, err);
128 		ret = nn_copy(&T[1], &T[1 + (expbit ^ rbit_next)]); EG(ret, err);
129 		/* Update rbit */
130 		rbit = rbit_next;
131 	}
132 	ret = nn_one(&T[1 - rbit]);
133 	if(r != NULL){
134 		/* Unredcify in out */
135 		ret = nn_mul_redc1(&T[rbit], &T[rbit], &T[1 - rbit], mod, mpinv); EG(ret, err);
136 	}
137 
138 	/* Deal with the particular cases of 0 and 1 bit exponents */
139 	/* Case with 0 bit exponent: T[1 - rbit] contains 1 modulo mod */
140 	ret = nn_mod(&T[1 - rbit], &T[1 - rbit], mod); EG(ret, err);
141 	/* Case with 1 bit exponent */
142 	ret = nn_mod(&T[2], base, mod); EG(ret, err);
143 	/* Proceed with the output */
144 	ret = nn_cnd_swap((oldexplen == 0), out, &T[1 - rbit]);
145 	ret = nn_cnd_swap((oldexplen == 1), out, &T[2]);
146 	ret = nn_cnd_swap(((oldexplen != 0) && (oldexplen != 1)), out, &T[rbit]);
147 
148 err:
149 	nn_uninit(&T[0]);
150 	nn_uninit(&T[1]);
151 	nn_uninit(&T[2]);
152 	nn_uninit(&mask);
153 
154 	return ret;
155 }
156 
157 /*
158  * NOT constant time with regard to the bitlength of exp.
159  *
160  * Reduces the base modulo mod if it is not already reduced,
161  * which is also a small divergence wrt constant time leaking
162  * the information that base <= mod or not: please use with care
163  * in the callers if this information is sensitive.
164  *
165  * Aliasing not supported. Expects caller to check parameters
166  * have been initialized. This is an internal helper.
167  *
168  * Compute (base ** exp) mod (mod) using a Montgomery Ladder algorithm
169  * with Montgomery redcification, hence the Montgomery coefficients as input.
170  * The module "mod" is expected to be odd for redcification to be used.
171  *
172  * Returns 0 on success, -1 on error.
173  */
_nn_mod_pow_redc(nn_t out,nn_src_t base,nn_src_t exp,nn_src_t mod,nn_src_t r,nn_src_t r_square,word_t mpinv)174 ATTRIBUTE_WARN_UNUSED_RET static int _nn_mod_pow_redc(nn_t out, nn_src_t base, nn_src_t exp, nn_src_t mod, nn_src_t r, nn_src_t r_square, word_t mpinv)
175 {
176 	return _nn_exp_monty_ladder_ltr(out, base, exp, mod, r, r_square, mpinv);
177 }
178 
179 /*
180  * NOT constant time with regard to the bitlength of exp.
181  *
182  * Reduces the base modulo mod if it is not already reduced,
183  * which is also a small divergence wrt constant time leaking
184  * the information that base <= mod or not: please use with care
185  * in the callers if this information is sensitive.
186  *
187  * Aliasing is supported. Expects caller to check parameters
188  * have been initialized. This is an internal helper.
189  *
190  * Compute (base ** exp) mod (mod) using a Montgomery Ladder algorithm.
191  * This function works for all values of "mod", but is slower that the one
192  * using Montgomery multiplication (which only works for odd "mod"). Hence,
193  * it is only used on even "mod" by upper layers.
194  *
195  * Returns 0 on success, -1 on error.
196  */
_nn_mod_pow(nn_t out,nn_src_t base,nn_src_t exp,nn_src_t mod)197 ATTRIBUTE_WARN_UNUSED_RET static int _nn_mod_pow(nn_t out, nn_src_t base, nn_src_t exp, nn_src_t mod)
198 {
199 	int ret;
200 
201 	if ((out == base) || (out == exp) || (out == mod)) {
202 		nn _out;
203 		_out.magic = WORD(0);
204 
205 		ret = nn_init(&_out, 0); EG(ret, err);
206 		ret = _nn_exp_monty_ladder_ltr(&_out, base, exp, mod, NULL, NULL, WORD(0)); EG(ret, err);
207 		ret = nn_copy(out, &_out);
208 	}
209 	else{
210 		ret = _nn_exp_monty_ladder_ltr(out, base, exp, mod, NULL, NULL, WORD(0));
211 	}
212 
213 err:
214 	return ret;
215 }
216 
217 /*
218  * Same purpose as above but handles aliasing.
219  * Expects caller to check parameters
220  * have been initialized. This is an internal helper.
221  */
_nn_mod_pow_redc_aliased(nn_t out,nn_src_t base,nn_src_t exp,nn_src_t mod,nn_src_t r,nn_src_t r_square,word_t mpinv)222 ATTRIBUTE_WARN_UNUSED_RET static int _nn_mod_pow_redc_aliased(nn_t out, nn_src_t base, nn_src_t exp, nn_src_t mod, nn_src_t r, nn_src_t r_square, word_t mpinv)
223 {
224 	nn _out;
225 	int ret;
226 	_out.magic = WORD(0);
227 
228 	ret = nn_init(&_out, 0); EG(ret, err);
229 	ret = _nn_mod_pow_redc(&_out, base, exp, mod, r, r_square, mpinv); EG(ret, err);
230 	ret = nn_copy(out, &_out);
231 
232 err:
233 	nn_uninit(&_out);
234 
235 	return ret;
236 }
237 
238 /* Aliased version of previous one.
239  * NOTE: our nn_mod_pow_redc primitives suppose that the modulo is odd for Montgomery multiplication
240  * primitives to provide consistent results.
241  *
242  * Aliasing is supported.
243  */
nn_mod_pow_redc(nn_t out,nn_src_t base,nn_src_t exp,nn_src_t mod,nn_src_t r,nn_src_t r_square,word_t mpinv)244 int nn_mod_pow_redc(nn_t out, nn_src_t base, nn_src_t exp, nn_src_t mod, nn_src_t r, nn_src_t r_square, word_t mpinv)
245 {
246 	int ret, isodd;
247 
248 	ret = nn_check_initialized(base); EG(ret, err);
249 	ret = nn_check_initialized(exp); EG(ret, err);
250 	ret = nn_check_initialized(mod); EG(ret, err);
251 	ret = nn_check_initialized(r); EG(ret, err);
252 	ret = nn_check_initialized(r_square); EG(ret, err);
253 
254 	/* Check that we have an odd number for our modulo */
255 	ret = nn_isodd(mod, &isodd); EG(ret, err);
256 	MUST_HAVE(isodd, ret, err);
257 
258 	/* Handle the case where our prime is less than two words size.
259 	 * We need it to be >= 2 words size for the Montgomery multiplication to be
260 	 * usable.
261 	 */
262 	if(mod->wlen < 2){
263 		/* Local copy our modulo */
264 		nn _mod;
265 		_mod.magic = WORD(0);
266 
267 		/* And set its length accordingly */
268 		ret = nn_copy(&_mod, mod); EG(ret, err1);
269 		ret = nn_set_wlen(&_mod, 2); EG(ret, err1);
270 		/* Handle output aliasing */
271 		if ((out == base) || (out == exp) || (out == mod) || (out == r) || (out == r_square)) {
272 			ret = _nn_mod_pow_redc_aliased(out, base, exp, &_mod, r, r_square, mpinv); EG(ret, err1);
273 		} else {
274 			ret = _nn_mod_pow_redc(out, base, exp, &_mod, r, r_square, mpinv); EG(ret, err1);
275 		}
276 err1:
277 		nn_uninit(&_mod);
278 		EG(ret, err);
279 	}
280 	else{
281 		/* Handle output aliasing */
282 		if ((out == base) || (out == exp) || (out == mod) || (out == r) || (out == r_square)) {
283 			ret = _nn_mod_pow_redc_aliased(out, base, exp, mod, r, r_square, mpinv);
284 		} else {
285 			ret = _nn_mod_pow_redc(out, base, exp, mod, r, r_square, mpinv);
286 		}
287 	}
288 
289 err:
290 	return ret;
291 }
292 
293 
294 /*
295  * NOT constant time with regard to the bitlength of exp.
296  * Aliasing is supported.
297  *
298  * Compute (base ** exp) mod (mod) using a Montgomery Ladder algorithm.
299  * Internally, this computes Montgomery coefficients and uses the redc
300  * function.
301  *
302  * Returns 0 on success, -1 on error.
303  */
nn_mod_pow(nn_t out,nn_src_t base,nn_src_t exp,nn_src_t mod)304 int nn_mod_pow(nn_t out, nn_src_t base, nn_src_t exp, nn_src_t mod)
305 {
306 	nn r, r_square;
307 	word_t mpinv;
308 	int ret, isodd;
309 	r.magic = r_square.magic = WORD(0);
310 
311 	/* Handle the case where our modulo is even: in this case, we cannot
312 	 * use our Montgomery multiplication primitives as they are only suitable
313 	 * for odd numbers. We fallback to less efficient regular modular exponentiation.
314 	 */
315 	ret = nn_isodd(mod, &isodd); EG(ret, err);
316 	if(!isodd){
317 		/* mod is even: use the regular unoptimized modular exponentiation */
318 		ret = _nn_mod_pow(out, base, exp, mod);
319 	}
320 	else{
321 		/* mod is odd */
322 		/* Compute the Montgomery coefficients */
323 		ret = nn_compute_redc1_coefs(&r, &r_square, mod, &mpinv); EG(ret, err);
324 
325 		/* Now use the redc version */
326 		ret = nn_mod_pow_redc(out, base, exp, mod, &r, &r_square, mpinv);
327 	}
328 
329 err:
330 	nn_uninit(&r);
331 	nn_uninit(&r_square);
332 
333 	return ret;
334 }
335