crypto/bn/bn_gcd.c

   1 /*
   2  * Copyright 1995-2020 The OpenSSL Project Authors. All Rights Reserved.
   3  *
   4  * Licensed under the Apache License 2.0 (the "License").  You may not use
   5  * this file except in compliance with the License.  You can obtain a copy
   6  * in the file LICENSE in the source distribution or at
   7  * https://www.openssl.org/source/license.html
   8  */
   9
  10 #include "internal/cryptlib.h"
  11 #include "bn_local.h"
  12
  13 /*
  14  * bn_mod_inverse_no_branch is a special version of BN_mod_inverse. It does
  15  * not contain branches that may leak sensitive information.
  16  *
  17  * This is a static function, we ensure all callers in this file pass valid
  18  * arguments: all passed pointers here are non-NULL.
  19  */
  20 static ossl_inline
  21 BIGNUM *bn_mod_inverse_no_branch(BIGNUM *in,
  22                                  const BIGNUM *a, const BIGNUM *n,
  23                                  BN_CTX *ctx, int *pnoinv)
  24 {
  25     BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
  26     BIGNUM *ret = NULL;
  27     int sign;
  28
  29     bn_check_top(a);
  30     bn_check_top(n);
  31
  32     BN_CTX_start(ctx);
  33     A = BN_CTX_get(ctx);
  34     B = BN_CTX_get(ctx);
  35     X = BN_CTX_get(ctx);
  36     D = BN_CTX_get(ctx);
  37     M = BN_CTX_get(ctx);
  38     Y = BN_CTX_get(ctx);
  39     T = BN_CTX_get(ctx);
  40     if (T == NULL)
  41         goto err;
  42
  43     if (in == NULL)
  44         R = BN_new();
  45     else
  46         R = in;
  47     if (R == NULL)
  48         goto err;
  49
  50     BN_one(X);
  51     BN_zero(Y);
  52     if (BN_copy(B, a) == NULL)
  53         goto err;
  54     if (BN_copy(A, n) == NULL)
  55         goto err;
  56     A->neg = 0;
  57
  58     if (B->neg || (BN_ucmp(B, A) >= 0)) {
  59         /*
  60          * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
  61          * BN_div_no_branch will be called eventually.
  62          */
  63          {
  64             BIGNUM local_B;
  65             bn_init(&local_B);
  66             BN_with_flags(&local_B, B, BN_FLG_CONSTTIME);
  67             if (!BN_nnmod(B, &local_B, A, ctx))
  68                 goto err;
  69             /* Ensure local_B goes out of scope before any further use of B */
  70         }
  71     }
  72     sign = -1;
  73     /*-
  74      * From  B = a mod |n|,  A = |n|  it follows that
  75      *
  76      *      0 <= B < A,
  77      *     -sign*X*a  ==  B   (mod |n|),
  78      *      sign*Y*a  ==  A   (mod |n|).
  79      */
  80
  81     while (!BN_is_zero(B)) {
  82         BIGNUM *tmp;
  83
  84         /*-
  85          *      0 < B < A,
  86          * (*) -sign*X*a  ==  B   (mod |n|),
  87          *      sign*Y*a  ==  A   (mod |n|)
  88          */
  89
  90         /*
  91          * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
  92          * BN_div_no_branch will be called eventually.
  93          */
  94         {
  95             BIGNUM local_A;
  96             bn_init(&local_A);
  97             BN_with_flags(&local_A, A, BN_FLG_CONSTTIME);
  98
  99             /* (D, M) := (A/B, A%B) ... */
 100             if (!BN_div(D, M, &local_A, B, ctx))
 101                 goto err;
 102             /* Ensure local_A goes out of scope before any further use of A */
 103         }
 104
 105         /*-
 106          * Now
 107          *      A = D*B + M;
 108          * thus we have
 109          * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
 110          */
 111
 112         tmp = A;                /* keep the BIGNUM object, the value does not
 113                                  * matter */
 114
 115         /* (A, B) := (B, A mod B) ... */
 116         A = B;
 117         B = M;
 118         /* ... so we have  0 <= B < A  again */
 119
 120         /*-
 121          * Since the former  M  is now  B  and the former  B  is now  A,
 122          * (**) translates into
 123          *       sign*Y*a  ==  D*A + B    (mod |n|),
 124          * i.e.
 125          *       sign*Y*a - D*A  ==  B    (mod |n|).
 126          * Similarly, (*) translates into
 127          *      -sign*X*a  ==  A          (mod |n|).
 128          *
 129          * Thus,
 130          *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
 131          * i.e.
 132          *        sign*(Y + D*X)*a  ==  B  (mod |n|).
 133          *
 134          * So if we set  (X, Y, sign) := (Y + D*X, X, -sign), we arrive back at
 135          *      -sign*X*a  ==  B   (mod |n|),
 136          *       sign*Y*a  ==  A   (mod |n|).
 137          * Note that  X  and  Y  stay non-negative all the time.
 138          */
 139
 140         if (!BN_mul(tmp, D, X, ctx))
 141             goto err;
 142         if (!BN_add(tmp, tmp, Y))
 143             goto err;
 144
 145         M = Y;                  /* keep the BIGNUM object, the value does not
 146                                  * matter */
 147         Y = X;
 148         X = tmp;
 149         sign = -sign;
 150     }
 151
 152     /*-
 153      * The while loop (Euclid's algorithm) ends when
 154      *      A == gcd(a,n);
 155      * we have
 156      *       sign*Y*a  ==  A  (mod |n|),
 157      * where  Y  is non-negative.
 158      */
 159
 160     if (sign < 0) {
 161         if (!BN_sub(Y, n, Y))
 162             goto err;
 163     }
 164     /* Now  Y*a  ==  A  (mod |n|).  */
 165
 166     if (BN_is_one(A)) {
 167         /* Y*a == 1  (mod |n|) */
 168         if (!Y->neg && BN_ucmp(Y, n) < 0) {
 169             if (!BN_copy(R, Y))
 170                 goto err;
 171         } else {
 172             if (!BN_nnmod(R, Y, n, ctx))
 173                 goto err;
 174         }
 175     } else {
 176         *pnoinv = 1;
 177         /* caller sets the BN_R_NO_INVERSE error */
 178         goto err;
 179     }
 180
 181     ret = R;
 182     *pnoinv = 0;
 183
 184  err:
 185     if ((ret == NULL) && (in == NULL))
 186         BN_free(R);
 187     BN_CTX_end(ctx);
 188     bn_check_top(ret);
 189     return ret;
 190 }
 191
 192 /*
 193  * This is an internal function, we assume all callers pass valid arguments:
 194  * all pointers passed here are assumed non-NULL.
 195  */
 196 BIGNUM *int_bn_mod_inverse(BIGNUM *in,
 197                            const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx,
 198                            int *pnoinv)
 199 {
 200     BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
 201     BIGNUM *ret = NULL;
 202     int sign;
 203
 204     /* This is invalid input so we don't worry about constant time here */
 205     if (BN_abs_is_word(n, 1) || BN_is_zero(n)) {
 206         *pnoinv = 1;
 207         return NULL;
 208     }
 209
 210     *pnoinv = 0;
 211
 212     if ((BN_get_flags(a, BN_FLG_CONSTTIME) != 0)
 213         || (BN_get_flags(n, BN_FLG_CONSTTIME) != 0)) {
 214         return bn_mod_inverse_no_branch(in, a, n, ctx, pnoinv);
 215     }
 216
 217     bn_check_top(a);
 218     bn_check_top(n);
 219
 220     BN_CTX_start(ctx);
 221     A = BN_CTX_get(ctx);
 222     B = BN_CTX_get(ctx);
 223     X = BN_CTX_get(ctx);
 224     D = BN_CTX_get(ctx);
 225     M = BN_CTX_get(ctx);
 226     Y = BN_CTX_get(ctx);
 227     T = BN_CTX_get(ctx);
 228     if (T == NULL)
 229         goto err;
 230
 231     if (in == NULL)
 232         R = BN_new();
 233     else
 234         R = in;
 235     if (R == NULL)
 236         goto err;
 237
 238     BN_one(X);
 239     BN_zero(Y);
 240     if (BN_copy(B, a) == NULL)
 241         goto err;
 242     if (BN_copy(A, n) == NULL)
 243         goto err;
 244     A->neg = 0;
 245     if (B->neg || (BN_ucmp(B, A) >= 0)) {
 246         if (!BN_nnmod(B, B, A, ctx))
 247             goto err;
 248     }
 249     sign = -1;
 250     /*-
 251      * From  B = a mod |n|,  A = |n|  it follows that
 252      *
 253      *      0 <= B < A,
 254      *     -sign*X*a  ==  B   (mod |n|),
 255      *      sign*Y*a  ==  A   (mod |n|).
 256      */
 257
 258     if (BN_is_odd(n) && (BN_num_bits(n) <= 2048)) {
 259         /*
 260          * Binary inversion algorithm; requires odd modulus. This is faster
 261          * than the general algorithm if the modulus is sufficiently small
 262          * (about 400 .. 500 bits on 32-bit systems, but much more on 64-bit
 263          * systems)
 264          */
 265         int shift;
 266
 267         while (!BN_is_zero(B)) {
 268             /*-
 269              *      0 < B < |n|,
 270              *      0 < A <= |n|,
 271              * (1) -sign*X*a  ==  B   (mod |n|),
 272              * (2)  sign*Y*a  ==  A   (mod |n|)
 273              */
 274
 275             /*
 276              * Now divide B by the maximum possible power of two in the
 277              * integers, and divide X by the same value mod |n|. When we're
 278              * done, (1) still holds.
 279              */
 280             shift = 0;
 281             while (!BN_is_bit_set(B, shift)) { /* note that 0 < B */
 282                 shift++;
 283
 284                 if (BN_is_odd(X)) {
 285                     if (!BN_uadd(X, X, n))
 286                         goto err;
 287                 }
 288                 /*
 289                  * now X is even, so we can easily divide it by two
 290                  */
 291                 if (!BN_rshift1(X, X))
 292                     goto err;
 293             }
 294             if (shift > 0) {
 295                 if (!BN_rshift(B, B, shift))
 296                     goto err;
 297             }
 298
 299             /*
 300              * Same for A and Y.  Afterwards, (2) still holds.
 301              */
 302             shift = 0;
 303             while (!BN_is_bit_set(A, shift)) { /* note that 0 < A */
 304                 shift++;
 305
 306                 if (BN_is_odd(Y)) {
 307                     if (!BN_uadd(Y, Y, n))
 308                         goto err;
 309                 }
 310                 /* now Y is even */
 311                 if (!BN_rshift1(Y, Y))
 312                     goto err;
 313             }
 314             if (shift > 0) {
 315                 if (!BN_rshift(A, A, shift))
 316                     goto err;
 317             }
 318
 319             /*-
 320              * We still have (1) and (2).
 321              * Both  A  and  B  are odd.
 322              * The following computations ensure that
 323              *
 324              *     0 <= B < |n|,
 325              *      0 < A < |n|,
 326              * (1) -sign*X*a  ==  B   (mod |n|),
 327              * (2)  sign*Y*a  ==  A   (mod |n|),
 328              *
 329              * and that either  A  or  B  is even in the next iteration.
 330              */
 331             if (BN_ucmp(B, A) >= 0) {
 332                 /* -sign*(X + Y)*a == B - A  (mod |n|) */
 333                 if (!BN_uadd(X, X, Y))
 334                     goto err;
 335                 /*
 336                  * NB: we could use BN_mod_add_quick(X, X, Y, n), but that
 337                  * actually makes the algorithm slower
 338                  */
 339                 if (!BN_usub(B, B, A))
 340                     goto err;
 341             } else {
 342                 /*  sign*(X + Y)*a == A - B  (mod |n|) */
 343                 if (!BN_uadd(Y, Y, X))
 344                     goto err;
 345                 /*
 346                  * as above, BN_mod_add_quick(Y, Y, X, n) would slow things down
 347                  */
 348                 if (!BN_usub(A, A, B))
 349                     goto err;
 350             }
 351         }
 352     } else {
 353         /* general inversion algorithm */
 354
 355         while (!BN_is_zero(B)) {
 356             BIGNUM *tmp;
 357
 358             /*-
 359              *      0 < B < A,
 360              * (*) -sign*X*a  ==  B   (mod |n|),
 361              *      sign*Y*a  ==  A   (mod |n|)
 362              */
 363
 364             /* (D, M) := (A/B, A%B) ... */
 365             if (BN_num_bits(A) == BN_num_bits(B)) {
 366                 if (!BN_one(D))
 367                     goto err;
 368                 if (!BN_sub(M, A, B))
 369                     goto err;
 370             } else if (BN_num_bits(A) == BN_num_bits(B) + 1) {
 371                 /* A/B is 1, 2, or 3 */
 372                 if (!BN_lshift1(T, B))
 373                     goto err;
 374                 if (BN_ucmp(A, T) < 0) {
 375                     /* A < 2*B, so D=1 */
 376                     if (!BN_one(D))
 377                         goto err;
 378                     if (!BN_sub(M, A, B))
 379                         goto err;
 380                 } else {
 381                     /* A >= 2*B, so D=2 or D=3 */
 382                     if (!BN_sub(M, A, T))
 383                         goto err;
 384                     if (!BN_add(D, T, B))
 385                         goto err; /* use D (:= 3*B) as temp */
 386                     if (BN_ucmp(A, D) < 0) {
 387                         /* A < 3*B, so D=2 */
 388                         if (!BN_set_word(D, 2))
 389                             goto err;
 390                         /*
 391                          * M (= A - 2*B) already has the correct value
 392                          */
 393                     } else {
 394                         /* only D=3 remains */
 395                         if (!BN_set_word(D, 3))
 396                             goto err;
 397                         /*
 398                          * currently M = A - 2*B, but we need M = A - 3*B
 399                          */
 400                         if (!BN_sub(M, M, B))
 401                             goto err;
 402                     }
 403                 }
 404             } else {
 405                 if (!BN_div(D, M, A, B, ctx))
 406                     goto err;
 407             }
 408
 409             /*-
 410              * Now
 411              *      A = D*B + M;
 412              * thus we have
 413              * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
 414              */
 415
 416             tmp = A;    /* keep the BIGNUM object, the value does not matter */
 417
 418             /* (A, B) := (B, A mod B) ... */
 419             A = B;
 420             B = M;
 421             /* ... so we have  0 <= B < A  again */
 422
 423             /*-
 424              * Since the former  M  is now  B  and the former  B  is now  A,
 425              * (**) translates into
 426              *       sign*Y*a  ==  D*A + B    (mod |n|),
 427              * i.e.
 428              *       sign*Y*a - D*A  ==  B    (mod |n|).
 429              * Similarly, (*) translates into
 430              *      -sign*X*a  ==  A          (mod |n|).
 431              *
 432              * Thus,
 433              *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
 434              * i.e.
 435              *        sign*(Y + D*X)*a  ==  B  (mod |n|).
 436              *
 437              * So if we set  (X, Y, sign) := (Y + D*X, X, -sign), we arrive back at
 438              *      -sign*X*a  ==  B   (mod |n|),
 439              *       sign*Y*a  ==  A   (mod |n|).
 440              * Note that  X  and  Y  stay non-negative all the time.
 441              */
 442
 443             /*
 444              * most of the time D is very small, so we can optimize tmp := D*X+Y
 445              */
 446             if (BN_is_one(D)) {
 447                 if (!BN_add(tmp, X, Y))
 448                     goto err;
 449             } else {
 450                 if (BN_is_word(D, 2)) {
 451                     if (!BN_lshift1(tmp, X))
 452                         goto err;
 453                 } else if (BN_is_word(D, 4)) {
 454                     if (!BN_lshift(tmp, X, 2))
 455                         goto err;
 456                 } else if (D->top == 1) {
 457                     if (!BN_copy(tmp, X))
 458                         goto err;
 459                     if (!BN_mul_word(tmp, D->d[0]))
 460                         goto err;
 461                 } else {
 462                     if (!BN_mul(tmp, D, X, ctx))
 463                         goto err;
 464                 }
 465                 if (!BN_add(tmp, tmp, Y))
 466                     goto err;
 467             }
 468
 469             M = Y;      /* keep the BIGNUM object, the value does not matter */
 470             Y = X;
 471             X = tmp;
 472             sign = -sign;
 473         }
 474     }
 475
 476     /*-
 477      * The while loop (Euclid's algorithm) ends when
 478      *      A == gcd(a,n);
 479      * we have
 480      *       sign*Y*a  ==  A  (mod |n|),
 481      * where  Y  is non-negative.
 482      */
 483
 484     if (sign < 0) {
 485         if (!BN_sub(Y, n, Y))
 486             goto err;
 487     }
 488     /* Now  Y*a  ==  A  (mod |n|).  */
 489
 490     if (BN_is_one(A)) {
 491         /* Y*a == 1  (mod |n|) */
 492         if (!Y->neg && BN_ucmp(Y, n) < 0) {
 493             if (!BN_copy(R, Y))
 494                 goto err;
 495         } else {
 496             if (!BN_nnmod(R, Y, n, ctx))
 497                 goto err;
 498         }
 499     } else {
 500         *pnoinv = 1;
 501         goto err;
 502     }
 503     ret = R;
 504  err:
 505     if ((ret == NULL) && (in == NULL))
 506         BN_free(R);
 507     BN_CTX_end(ctx);
 508     bn_check_top(ret);
 509     return ret;
 510 }
 511
 512 /* solves ax == 1 (mod n) */
 513 BIGNUM *BN_mod_inverse(BIGNUM *in,
 514                        const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
 515 {
 516     BN_CTX *new_ctx = NULL;
 517     BIGNUM *rv;
 518     int noinv = 0;
 519
 520     if (ctx == NULL) {
 521         ctx = new_ctx = BN_CTX_new_ex(NULL);
 522         if (ctx == NULL) {
 523             ERR_raise(ERR_LIB_BN, ERR_R_MALLOC_FAILURE);
 524             return NULL;
 525         }
 526     }
 527
 528     rv = int_bn_mod_inverse(in, a, n, ctx, &noinv);
 529     if (noinv)
 530         ERR_raise(ERR_LIB_BN, BN_R_NO_INVERSE);
 531     BN_CTX_free(new_ctx);
 532     return rv;
 533 }
 534
 535 /*-
 536  * This function is based on the constant-time GCD work by Bernstein and Yang:
 537  * https://eprint.iacr.org/2019/266
 538  * Generalized fast GCD function to allow even inputs.
 539  * The algorithm first finds the shared powers of 2 between
 540  * the inputs, and removes them, reducing at least one of the
 541  * inputs to an odd value. Then it proceeds to calculate the GCD.
 542  * Before returning the resulting GCD, we take care of adding
 543  * back the powers of two removed at the beginning.
 544  * Note 1: we assume the bit length of both inputs is public information,
 545  * since access to top potentially leaks this information.
 546  */
 547 int BN_gcd(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
 548 {
 549     BIGNUM *g, *temp = NULL;
 550     BN_ULONG mask = 0;
 551     int i, j, top, rlen, glen, m, bit = 1, delta = 1, cond = 0, shifts = 0, ret = 0;
 552
 553     /* Note 2: zero input corner cases are not constant-time since they are
 554      * handled immediately. An attacker can run an attack under this
 555      * assumption without the need of side-channel information. */
 556     if (BN_is_zero(in_b)) {
 557         ret = BN_copy(r, in_a) != NULL;
 558         r->neg = 0;
 559         return ret;
 560     }
 561     if (BN_is_zero(in_a)) {
 562         ret = BN_copy(r, in_b) != NULL;
 563         r->neg = 0;
 564         return ret;
 565     }
 566
 567     bn_check_top(in_a);
 568     bn_check_top(in_b);
 569
 570     BN_CTX_start(ctx);
 571     temp = BN_CTX_get(ctx);
 572     g = BN_CTX_get(ctx);
 573
 574     /* make r != 0, g != 0 even, so BN_rshift is not a potential nop */
 575     if (g == NULL
 576         || !BN_lshift1(g, in_b)
 577         || !BN_lshift1(r, in_a))
 578         goto err;
 579
 580     /* find shared powers of two, i.e. "shifts" >= 1 */
 581     for (i = 0; i < r->dmax && i < g->dmax; i++) {
 582         mask = ~(r->d[i] | g->d[i]);
 583         for (j = 0; j < BN_BITS2; j++) {
 584             bit &= mask;
 585             shifts += bit;
 586             mask >>= 1;
 587         }
 588     }
 589
 590     /* subtract shared powers of two; shifts >= 1 */
 591     if (!BN_rshift(r, r, shifts)
 592         || !BN_rshift(g, g, shifts))
 593         goto err;
 594
 595     /* expand to biggest nword, with room for a possible extra word */
 596     top = 1 + ((r->top >= g->top) ? r->top : g->top);
 597     if (bn_wexpand(r, top) == NULL
 598         || bn_wexpand(g, top) == NULL
 599         || bn_wexpand(temp, top) == NULL)
 600         goto err;
 601
 602     /* re arrange inputs s.t. r is odd */
 603     BN_consttime_swap((~r->d[0]) & 1, r, g, top);
 604
 605     /* compute the number of iterations */
 606     rlen = BN_num_bits(r);
 607     glen = BN_num_bits(g);
 608     m = 4 + 3 * ((rlen >= glen) ? rlen : glen);
 609
 610     for (i = 0; i < m; i++) {
 611         /* conditionally flip signs if delta is positive and g is odd */
 612         cond = (-delta >> (8 * sizeof(delta) - 1)) & g->d[0] & 1
 613             /* make sure g->top > 0 (i.e. if top == 0 then g == 0 always) */
 614             & (~((g->top - 1) >> (sizeof(g->top) * 8 - 1)));
 615         delta = (-cond & -delta) | ((cond - 1) & delta);
 616         r->neg ^= cond;
 617         /* swap */
 618         BN_consttime_swap(cond, r, g, top);
 619
 620         /* elimination step */
 621         delta++;
 622         if (!BN_add(temp, g, r))
 623             goto err;
 624         BN_consttime_swap(g->d[0] & 1 /* g is odd */
 625                 /* make sure g->top > 0 (i.e. if top == 0 then g == 0 always) */
 626                 & (~((g->top - 1) >> (sizeof(g->top) * 8 - 1))),
 627                 g, temp, top);
 628         if (!BN_rshift1(g, g))
 629             goto err;
 630     }
 631
 632     /* remove possible negative sign */
 633     r->neg = 0;
 634     /* add powers of 2 removed, then correct the artificial shift */
 635     if (!BN_lshift(r, r, shifts)
 636         || !BN_rshift1(r, r))
 637         goto err;
 638
 639     ret = 1;
 640
 641  err:
 642     BN_CTX_end(ctx);
 643     bn_check_top(r);
 644     return ret;
 645 }