crypto/bn/bn_gcd.c

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