crypto/bn/bn_gcd.c

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