1 /*
2 * Copyright (C) 2017 - 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_modinv.h>
17 #include <libecc/nn/nn_div_public.h>
18 #include <libecc/nn/nn_logical.h>
19 #include <libecc/nn/nn_add.h>
20 #include <libecc/nn/nn_mod_pow.h>
21 #include <libecc/nn/nn.h>
22 /* Include the "internal" header as we use non public API here */
23 #include "../nn/nn_mul.h"
24
25 /*
26 * Compute out = x^-1 mod m, i.e. out such that (out * x) = 1 mod m
27 * out is initialized by the function, i.e. caller needs not initialize
28 * it; only provide the associated storage space. Done in *constant
29 * time* if underlying routines are.
30 *
31 * Asserts that m is odd and that x is smaller than m.
32 * This second condition is not strictly necessary,
33 * but it allows to perform all operations on nn's of the same length,
34 * namely the length of m.
35 *
36 * Uses a binary xgcd algorithm,
37 * only keeps track of coefficient for inverting x,
38 * and performs reduction modulo m at each step.
39 *
40 * This does not normalize out on return.
41 *
42 * 0 is returned on success (everything went ok and x has reciprocal), -1
43 * on error or if x has no reciprocal. On error, out is not meaningful.
44 *
45 * The function is an internal helper: caller MUST check params have been
46 * initialized, i.e. this is not done by the function.
47 *
48 */
_nn_modinv_odd(nn_t out,nn_src_t x,nn_src_t m)49 ATTRIBUTE_WARN_UNUSED_RET static int _nn_modinv_odd(nn_t out, nn_src_t x, nn_src_t m)
50 {
51 int isodd, swap, smaller, ret, cmp, iszero, tmp_isodd;
52 nn a, b, u, tmp, mp1d2;
53 nn_t uu = out;
54 bitcnt_t cnt;
55 a.magic = b.magic = u.magic = tmp.magic = mp1d2.magic = WORD(0);
56
57 ret = nn_init(out, 0); EG(ret, err);
58 ret = nn_init(&a, (u16)(m->wlen * WORD_BYTES)); EG(ret, err);
59 ret = nn_init(&b, (u16)(m->wlen * WORD_BYTES)); EG(ret, err);
60 ret = nn_init(&u, (u16)(m->wlen * WORD_BYTES)); EG(ret, err);
61 ret = nn_init(&mp1d2, (u16)(m->wlen * WORD_BYTES)); EG(ret, err);
62 /*
63 * Temporary space needed to only deal with positive stuff.
64 */
65 ret = nn_init(&tmp, (u16)(m->wlen * WORD_BYTES)); EG(ret, err);
66
67 MUST_HAVE((!nn_isodd(m, &isodd)) && isodd, ret, err);
68 MUST_HAVE((!nn_cmp(x, m, &cmp)) && (cmp < 0), ret, err);
69 MUST_HAVE((!nn_iszero(x, &iszero)) && (!iszero), ret, err);
70
71 /*
72 * Maintain:
73 *
74 * a = u * x (mod m)
75 * b = uu * x (mod m)
76 *
77 * and b odd at all times. Initially,
78 *
79 * a = x, u = 1
80 * b = m, uu = 0
81 */
82 ret = nn_copy(&a, x); EG(ret, err);
83 ret = nn_set_wlen(&a, m->wlen); EG(ret, err);
84 ret = nn_copy(&b, m); EG(ret, err);
85 ret = nn_one(&u); EG(ret, err);
86 ret = nn_zero(uu); EG(ret, err);
87
88 /*
89 * The lengths of u and uu should not affect constant timeness but it
90 * does not hurt to set them already.
91 * They will always be strictly smaller than m.
92 */
93 ret = nn_set_wlen(&u, m->wlen); EG(ret, err);
94 ret = nn_set_wlen(uu, m->wlen); EG(ret, err);
95
96 /*
97 * Precompute inverse of 2 mod m:
98 * 2^-1 = (m+1)/2
99 * computed as (m >> 1) + 1.
100 */
101 ret = nn_rshift_fixedlen(&mp1d2, m, 1); EG(ret, err);
102
103 ret = nn_inc(&mp1d2, &mp1d2); EG(ret, err); /* no carry can occur here
104 because of prev. shift */
105
106 cnt = (bitcnt_t)((a.wlen + b.wlen) * WORD_BITS);
107 while (cnt > 0) {
108 cnt = (bitcnt_t)(cnt - 1);
109 /*
110 * Always maintain b odd. The logic of the iteration is as
111 * follows.
112 */
113
114 /*
115 * For a, b:
116 *
117 * odd = a & 1
118 * swap = odd & (a < b)
119 * if (swap)
120 * swap(a, b)
121 * if (odd)
122 * a -= b
123 * a /= 2
124 */
125
126 MUST_HAVE((!nn_isodd(&b, &tmp_isodd)) && tmp_isodd, ret, err);
127
128 ret = nn_isodd(&a, &isodd); EG(ret, err);
129 ret = nn_cmp(&a, &b, &cmp); EG(ret, err);
130 swap = isodd & (cmp == -1);
131
132 ret = nn_cnd_swap(swap, &a, &b); EG(ret, err);
133 ret = nn_cnd_sub(isodd, &a, &a, &b); EG(ret, err);
134
135 MUST_HAVE((!nn_isodd(&a, &tmp_isodd)) && (!tmp_isodd), ret, err); /* a is now even */
136
137 ret = nn_rshift_fixedlen(&a, &a, 1); EG(ret, err);/* division by 2 */
138
139 /*
140 * For u, uu:
141 *
142 * if (swap)
143 * swap u, uu
144 * smaller = (u < uu)
145 * if (odd)
146 * if (smaller)
147 * u += m - uu
148 * else
149 * u -= uu
150 * u >>= 1
151 * if (u was odd)
152 * u += (m+1)/2
153 */
154 ret = nn_cnd_swap(swap, &u, uu); EG(ret, err);
155
156 /* This parameter is used to avoid handling negative numbers. */
157 ret = nn_cmp(&u, uu, &cmp); EG(ret, err);
158 smaller = (cmp == -1);
159
160 /* Computation of 'm - uu' can always be performed. */
161 ret = nn_sub(&tmp, m, uu); EG(ret, err);
162
163 /* Selection btw 'm-uu' and '-uu' is made by the following function calls. */
164 ret = nn_cnd_add(isodd & smaller, &u, &u, &tmp); EG(ret, err); /* no carry can occur as 'u+(m-uu) = m-(uu-u) < m' */
165 ret = nn_cnd_sub(isodd & (!smaller), &u, &u, uu); EG(ret, err);
166
167 /* Divide u by 2 */
168 ret = nn_isodd(&u, &isodd); EG(ret, err);
169 ret = nn_rshift_fixedlen(&u, &u, 1); EG(ret, err);
170 ret = nn_cnd_add(isodd, &u, &u, &mp1d2); EG(ret, err); /* no carry can occur as u=1+u' with u'<m-1 and u' even so u'/2+(m+1)/2<(m-1)/2+(m+1)/2=m */
171
172 MUST_HAVE((!nn_cmp(&u, m, &cmp)) && (cmp < 0), ret, err);
173 MUST_HAVE((!nn_cmp(uu, m, &cmp)) && (cmp < 0), ret, err);
174
175 /*
176 * As long as a > 0, the quantity
177 * (bitsize of a) + (bitsize of b)
178 * is reduced by at least one bit per iteration,
179 * hence after (bitsize of x) + (bitsize of m) - 1
180 * iterations we surely have a = 0. Then b = gcd(x, m)
181 * and if b = 1 then also uu = x^{-1} (mod m).
182 */
183 }
184
185 MUST_HAVE((!nn_iszero(&a, &iszero)) && iszero, ret, err);
186
187 /* Check that gcd is one. */
188 ret = nn_cmp_word(&b, WORD(1), &cmp); EG(ret, err);
189
190 /* If not, set "inverse" to zero. */
191 ret = nn_cnd_sub(cmp != 0, uu, uu, uu); EG(ret, err);
192
193 ret = cmp ? -1 : 0;
194
195 err:
196 nn_uninit(&a);
197 nn_uninit(&b);
198 nn_uninit(&u);
199 nn_uninit(&mp1d2);
200 nn_uninit(&tmp);
201
202 PTR_NULLIFY(uu);
203
204 return ret;
205 }
206
207 /*
208 * Same as above without restriction on m.
209 * No attempt to make it constant time.
210 * Uses the above constant-time binary xgcd when m is odd
211 * and a not constant time plain Euclidean xgcd when m is even.
212 *
213 * _out parameter need not be initialized; this will be done by the function.
214 * x and m must be initialized nn.
215 *
216 * Return -1 on error or if if x has no reciprocal modulo m. out is zeroed.
217 * Return 0 if x has reciprocal modulo m.
218 *
219 * The function supports aliasing.
220 */
nn_modinv(nn_t _out,nn_src_t x,nn_src_t m)221 int nn_modinv(nn_t _out, nn_src_t x, nn_src_t m)
222 {
223 int sign, ret, cmp, isodd, isone;
224 nn_t x_mod_m;
225 nn u, v, out; /* Out to support aliasing */
226 out.magic = u.magic = v.magic = WORD(0);
227
228 ret = nn_check_initialized(x); EG(ret, err);
229 ret = nn_check_initialized(m); EG(ret, err);
230
231 /* Initialize out */
232 ret = nn_init(&out, 0); EG(ret, err);
233 ret = nn_isodd(m, &isodd); EG(ret, err);
234 if (isodd) {
235 ret = nn_cmp(x, m, &cmp); EG(ret, err);
236 if (cmp >= 0) {
237 /*
238 * If x >= m, (x^-1) mod m = ((x mod m)^-1) mod m
239 * Hence, compute x mod m. In order to avoid
240 * additional stack usage, we use 'u' (not
241 * already useful at that point in the function).
242 */
243 x_mod_m = &u;
244 ret = nn_mod(x_mod_m, x, m); EG(ret, err);
245 ret = _nn_modinv_odd(&out, x_mod_m, m); EG(ret, err);
246 } else {
247 ret = _nn_modinv_odd(&out, x, m); EG(ret, err);
248 }
249 ret = nn_copy(_out, &out);
250 goto err;
251 }
252
253 /* Now m is even */
254 ret = nn_isodd(x, &isodd); EG(ret, err);
255 MUST_HAVE(isodd, ret, err);
256
257 ret = nn_init(&u, 0); EG(ret, err);
258 ret = nn_init(&v, 0); EG(ret, err);
259 ret = nn_xgcd(&out, &u, &v, x, m, &sign); EG(ret, err);
260 ret = nn_isone(&out, &isone); EG(ret, err);
261 MUST_HAVE(isone, ret, err);
262
263 ret = nn_mod(&out, &u, m); EG(ret, err);
264 if (sign == -1) {
265 ret = nn_sub(&out, m, &out); EG(ret, err);
266 }
267 ret = nn_copy(_out, &out);
268
269 err:
270 nn_uninit(&out);
271 nn_uninit(&u);
272 nn_uninit(&v);
273
274 PTR_NULLIFY(x_mod_m);
275
276 return ret;
277 }
278
279 /*
280 * Compute (A - B) % 2^(storagebitsizeof(B) + 1). A and B must be initialized nn.
281 * the function is an internal helper and does not verify params have been
282 * initialized; this must be done by the caller. No assumption on A and B values
283 * such as A >= B. Done in *constant time. Returns 0 on success, -1 on error.
284 */
_nn_sub_mod_2exp(nn_t A,nn_src_t B)285 ATTRIBUTE_WARN_UNUSED_RET static inline int _nn_sub_mod_2exp(nn_t A, nn_src_t B)
286 {
287 u8 Awlen = A->wlen;
288 int ret;
289
290 ret = nn_set_wlen(A, (u8)(Awlen + 1)); EG(ret, err);
291
292 /* Make sure A > B */
293 /* NOTE: A->wlen - 1 is not an issue here thant to the nn_set_wlen above */
294 A->val[A->wlen - 1] = WORD(1);
295 ret = nn_sub(A, A, B); EG(ret, err);
296
297 /* The artificial word will be cleared in the following function call */
298 ret = nn_set_wlen(A, Awlen);
299
300 err:
301 return ret;
302 }
303
304 /*
305 * Invert x modulo 2^exp using Hensel lifting. Returns 0 on success, -1 on
306 * error. On success, x_isodd is 1 if x is odd, 0 if it is even.
307 * Please note that the result is correct (inverse of x) only when x is prime
308 * to 2^exp, i.e. x is odd (x_odd is 1).
309 *
310 * Operations are done in *constant time*.
311 *
312 * Aliasing is supported.
313 */
nn_modinv_2exp(nn_t _out,nn_src_t x,bitcnt_t exp,int * x_isodd)314 int nn_modinv_2exp(nn_t _out, nn_src_t x, bitcnt_t exp, int *x_isodd)
315 {
316 bitcnt_t cnt;
317 u8 exp_wlen = (u8)BIT_LEN_WORDS(exp);
318 bitcnt_t exp_cnt = exp % WORD_BITS;
319 word_t mask = (word_t)((exp_cnt == 0) ? WORD_MASK : (word_t)((WORD(1) << exp_cnt) - WORD(1)));
320 nn tmp_sqr, tmp_mul;
321 /* for aliasing */
322 int isodd, ret;
323 nn out;
324 out.magic = tmp_sqr.magic = tmp_mul.magic = WORD(0);
325
326 MUST_HAVE((x_isodd != NULL), ret, err);
327 ret = nn_check_initialized(x); EG(ret, err);
328 ret = nn_check_initialized(_out); EG(ret, err);
329
330 ret = nn_init(&out, 0); EG(ret, err);
331 ret = nn_init(&tmp_sqr, 0); EG(ret, err);
332 ret = nn_init(&tmp_mul, 0); EG(ret, err);
333 ret = nn_isodd(x, &isodd); EG(ret, err);
334 if (exp == (bitcnt_t)0){
335 /* Specific case of zero exponent, output 0 */
336 (*x_isodd) = isodd;
337 goto err;
338 }
339 if (!isodd) {
340 ret = nn_zero(_out); EG(ret, err);
341 (*x_isodd) = 0;
342 goto err;
343 }
344
345 /*
346 * Inverse modulo 2.
347 */
348 cnt = 1;
349 ret = nn_one(&out); EG(ret, err);
350
351 /*
352 * Inverse modulo 2^(2^i) <= 2^WORD_BITS.
353 * Assumes WORD_BITS is a power of two.
354 */
355 for (; cnt < WORD_MIN(WORD_BITS, exp); cnt = (bitcnt_t)(cnt << 1)) {
356 ret = nn_sqr_low(&tmp_sqr, &out, out.wlen); EG(ret, err);
357 ret = nn_mul_low(&tmp_mul, &tmp_sqr, x, out.wlen); EG(ret, err);
358 ret = nn_lshift_fixedlen(&out, &out, 1); EG(ret, err);
359
360 /*
361 * Allowing "negative" results for a subtraction modulo
362 * a power of two would allow to use directly:
363 * nn_sub(out, out, tmp_mul)
364 * which is always negative in ZZ except when x is one.
365 *
366 * Another solution is to add the opposite of tmp_mul.
367 * nn_modopp_2exp(tmp_mul, tmp_mul);
368 * nn_add(out, out, tmp_mul);
369 *
370 * The current solution is to add a sufficiently large power
371 * of two to out unconditionally to absorb the potential
372 * borrow. The result modulo 2^(2^i) is correct whether the
373 * borrow occurs or not.
374 */
375 ret = _nn_sub_mod_2exp(&out, &tmp_mul); EG(ret, err);
376 }
377
378 /*
379 * Inverse modulo 2^WORD_BITS < 2^(2^i) < 2^exp.
380 */
381 for (; cnt < ((exp + 1) >> 1); cnt = (bitcnt_t)(cnt << 1)) {
382 ret = nn_set_wlen(&out, (u8)(2 * out.wlen)); EG(ret, err);
383 ret = nn_sqr_low(&tmp_sqr, &out, out.wlen); EG(ret, err);
384 ret = nn_mul_low(&tmp_mul, &tmp_sqr, x, out.wlen); EG(ret, err);
385 ret = nn_lshift_fixedlen(&out, &out, 1); EG(ret, err);
386 ret = _nn_sub_mod_2exp(&out, &tmp_mul); EG(ret, err);
387 }
388
389 /*
390 * Inverse modulo 2^(2^i + j) >= 2^exp.
391 */
392 if (exp > WORD_BITS) {
393 ret = nn_set_wlen(&out, exp_wlen); EG(ret, err);
394 ret = nn_sqr_low(&tmp_sqr, &out, out.wlen); EG(ret, err);
395 ret = nn_mul_low(&tmp_mul, &tmp_sqr, x, out.wlen); EG(ret, err);
396 ret = nn_lshift_fixedlen(&out, &out, 1); EG(ret, err);
397 ret = _nn_sub_mod_2exp(&out, &tmp_mul); EG(ret, err);
398 }
399
400 /*
401 * Inverse modulo 2^exp.
402 */
403 out.val[exp_wlen - 1] &= mask;
404
405 ret = nn_copy(_out, &out); EG(ret, err);
406
407 (*x_isodd) = 1;
408
409 err:
410 nn_uninit(&out);
411 nn_uninit(&tmp_sqr);
412 nn_uninit(&tmp_mul);
413
414 return ret;
415 }
416
417 /*
418 * Invert word w modulo m.
419 *
420 * The function supports aliasing.
421 */
nn_modinv_word(nn_t out,word_t w,nn_src_t m)422 int nn_modinv_word(nn_t out, word_t w, nn_src_t m)
423 {
424 nn nn_tmp;
425 int ret;
426 nn_tmp.magic = WORD(0);
427
428 ret = nn_init(&nn_tmp, 0); EG(ret, err);
429 ret = nn_set_word_value(&nn_tmp, w); EG(ret, err);
430 ret = nn_modinv(out, &nn_tmp, m);
431
432 err:
433 nn_uninit(&nn_tmp);
434
435 return ret;
436 }
437
438
439 /*
440 * Internal function for nn_modinv_fermat and nn_modinv_fermat_redc used
441 * hereafter.
442 */
_nn_modinv_fermat_common(nn_t out,nn_src_t x,nn_src_t p,nn_t p_minus_two,int * lesstwo)443 ATTRIBUTE_WARN_UNUSED_RET static int _nn_modinv_fermat_common(nn_t out, nn_src_t x, nn_src_t p, nn_t p_minus_two, int *lesstwo)
444 {
445 int ret, cmp, isodd;
446 nn two;
447 two.magic = WORD(0);
448
449 /* Sanity checks on inputs */
450 ret = nn_check_initialized(x); EG(ret, err);
451 ret = nn_check_initialized(p); EG(ret, err);
452 /* NOTE: since this is an internal function, we are ensured that p_minus_two,
453 * two and regular are OK.
454 */
455
456 /* 0 is not invertible in any case */
457 ret = nn_iszero(x, &cmp); EG(ret, err);
458 if(cmp){
459 /* Zero the output and return an error */
460 ret = nn_init(out, 0); EG(ret, err);
461 ret = nn_zero(out); EG(ret, err);
462 ret = -1;
463 goto err;
464 }
465
466 /* For p <= 2, p being prime either p = 1 or p = 2.
467 * When p = 2, only 1 has an inverse, if p = 1 no one has an inverse.
468 */
469 (*lesstwo) = 0;
470 ret = nn_cmp_word(p, WORD(2), &cmp); EG(ret, err);
471 if(cmp == 0){
472 /* This is the p = 2 case, parity of x provides the result */
473 ret = nn_isodd(x, &isodd); EG(ret, err);
474 if(isodd){
475 /* x is odd, 1 is its inverse */
476 ret = nn_init(out, 0); EG(ret, err);
477 ret = nn_one(out); EG(ret, err);
478 ret = 0;
479 }
480 else{
481 /* x is even, no inverse. Zero the output */
482 ret = nn_init(out, 0); EG(ret, err);
483 ret = nn_zero(out); EG(ret, err);
484 ret = -1;
485 }
486 (*lesstwo) = 1;
487 goto err;
488 } else if (cmp < 0){
489 /* This is the p = 1 case, no inverse here: hence return an error */
490 /* Zero the output */
491 ret = nn_init(out, 0); EG(ret, err);
492 ret = nn_zero(out); EG(ret, err);
493 ret = -1;
494 (*lesstwo) = 1;
495 goto err;
496 }
497
498 /* Else we compute (p-2) for the upper layer */
499 if(p != p_minus_two){
500 /* Handle aliasing of p and p_minus_two */
501 ret = nn_init(p_minus_two, 0); EG(ret, err);
502 }
503
504 ret = nn_init(&two, 0); EG(ret, err);
505 ret = nn_set_word_value(&two, WORD(2)); EG(ret, err);
506 ret = nn_sub(p_minus_two, p, &two);
507
508 err:
509 nn_uninit(&two);
510
511 return ret;
512 }
513
514 /*
515 * Invert NN x modulo p using Fermat's little theorem for our inversion:
516 *
517 * p prime means that:
518 * x^(p-1) = 1 mod (p)
519 * which means that x^(p-2) mod(p) is the modular inverse of x mod (p)
520 * for x != 0
521 *
522 * NOTE: the input hypothesis is that p is prime.
523 * XXX WARNING: using this function with p not prime will produce wrong
524 * results without triggering an error!
525 *
526 * The function returns 0 on success, -1 on error
527 * (e.g. if x has no inverse modulo p, i.e. x = 0).
528 *
529 * Aliasing is supported.
530 */
nn_modinv_fermat(nn_t out,nn_src_t x,nn_src_t p)531 int nn_modinv_fermat(nn_t out, nn_src_t x, nn_src_t p)
532 {
533 int ret, lesstwo;
534 nn p_minus_two;
535 p_minus_two.magic = WORD(0);
536
537 /* Call our helper.
538 * NOTE: "marginal" cases where x = 0 and p <= 2 should be caught in this helper.
539 */
540 ret = _nn_modinv_fermat_common(out, x, p, &p_minus_two, &lesstwo); EG(ret, err);
541
542 if(!lesstwo){
543 /* Compute x^(p-2) mod (p) */
544 ret = nn_mod_pow(out, x, &p_minus_two, p);
545 }
546
547 err:
548 nn_uninit(&p_minus_two);
549
550 return ret;
551 }
552
553 /*
554 * Invert NN x modulo m using Fermat's little theorem for our inversion.
555 *
556 * This is a version with already (pre)computed Montgomery coefficients.
557 *
558 * NOTE: the input hypothesis is that p is prime.
559 * XXX WARNING: using this function with p not prime will produce wrong
560 * results without triggering an error!
561 *
562 * The function returns 0 on success, -1 on error
563 * (e.g. if x has no inverse modulo p, i.e. x = 0).
564 *
565 * Aliasing is supported.
566 */
nn_modinv_fermat_redc(nn_t out,nn_src_t x,nn_src_t p,nn_src_t r,nn_src_t r_square,word_t mpinv)567 int nn_modinv_fermat_redc(nn_t out, nn_src_t x, nn_src_t p, nn_src_t r, nn_src_t r_square, word_t mpinv)
568 {
569 int ret, lesstwo;
570 nn p_minus_two;
571 p_minus_two.magic = WORD(0);
572
573 /* Call our helper.
574 * NOTE: "marginal" cases where x = 0 and p <= 2 should be caught in this helper.
575 */
576 ret = _nn_modinv_fermat_common(out, x, p, &p_minus_two, &lesstwo); EG(ret, err);
577
578 if(!lesstwo){
579 /* Compute x^(p-2) mod (p) using precomputed Montgomery coefficients as input */
580 ret = nn_mod_pow_redc(out, x, &p_minus_two, p, r, r_square, mpinv);
581 }
582
583 err:
584 nn_uninit(&p_minus_two);
585
586 return ret;
587 }
588