nn_modinv.c
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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 }