crypto/bn/bn_gcd.c

   1 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
   2  * All rights reserved.
   3  *
   4  * This package is an SSL implementation written
   5  * by Eric Young (eay@cryptsoft.com).
   6  * The implementation was written so as to conform with Netscapes SSL.
   7  *
   8  * This library is free for commercial and non-commercial use as long as
   9  * the following conditions are aheared to.  The following conditions
  10  * apply to all code found in this distribution, be it the RC4, RSA,
  11  * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
  12  * included with this distribution is covered by the same copyright terms
  13  * except that the holder is Tim Hudson (tjh@cryptsoft.com).
  14  *
  15  * Copyright remains Eric Young's, and as such any Copyright notices in
  16  * the code are not to be removed.
  17  * If this package is used in a product, Eric Young should be given attribution
  18  * as the author of the parts of the library used.
  19  * This can be in the form of a textual message at program startup or
  20  * in documentation (online or textual) provided with the package.
  21  *
  22  * Redistribution and use in source and binary forms, with or without
  23  * modification, are permitted provided that the following conditions
  24  * are met:
  25  * 1. Redistributions of source code must retain the copyright
  26  *    notice, this list of conditions and the following disclaimer.
  27  * 2. Redistributions in binary form must reproduce the above copyright
  28  *    notice, this list of conditions and the following disclaimer in the
  29  *    documentation and/or other materials provided with the distribution.
  30  * 3. All advertising materials mentioning features or use of this software
  31  *    must display the following acknowledgement:
  32  *    "This product includes cryptographic software written by
  33  *     Eric Young (eay@cryptsoft.com)"
  34  *    The word 'cryptographic' can be left out if the rouines from the library
  35  *    being used are not cryptographic related :-).
  36  * 4. If you include any Windows specific code (or a derivative thereof) from
  37  *    the apps directory (application code) you must include an acknowledgement:
  38  *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
  39  *
  40  * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND
  41  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  42  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  43  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  44  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  45  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  46  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  47  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  48  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  49  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  50  * SUCH DAMAGE.
  51  *
  52  * The licence and distribution terms for any publically available version or
  53  * derivative of this code cannot be changed.  i.e. this code cannot simply be
  54  * copied and put under another distribution licence
  55  * [including the GNU Public Licence.]
  56  */
  57 /* ====================================================================
  58  * Copyright (c) 1998-2001 The OpenSSL Project.  All rights reserved.
  59  *
  60  * Redistribution and use in source and binary forms, with or without
  61  * modification, are permitted provided that the following conditions
  62  * are met:
  63  *
  64  * 1. Redistributions of source code must retain the above copyright
  65  *    notice, this list of conditions and the following disclaimer.
  66  *
  67  * 2. Redistributions in binary form must reproduce the above copyright
  68  *    notice, this list of conditions and the following disclaimer in
  69  *    the documentation and/or other materials provided with the
  70  *    distribution.
  71  *
  72  * 3. All advertising materials mentioning features or use of this
  73  *    software must display the following acknowledgment:
  74  *    "This product includes software developed by the OpenSSL Project
  75  *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
  76  *
  77  * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
  78  *    endorse or promote products derived from this software without
  79  *    prior written permission. For written permission, please contact
  80  *    openssl-core@openssl.org.
  81  *
  82  * 5. Products derived from this software may not be called "OpenSSL"
  83  *    nor may "OpenSSL" appear in their names without prior written
  84  *    permission of the OpenSSL Project.
  85  *
  86  * 6. Redistributions of any form whatsoever must retain the following
  87  *    acknowledgment:
  88  *    "This product includes software developed by the OpenSSL Project
  89  *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
  90  *
  91  * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
  92  * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  93  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  94  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
  95  * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  96  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  97  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  98  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  99  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 100  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 101  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
 102  * OF THE POSSIBILITY OF SUCH DAMAGE.
 103  * ====================================================================
 104  *
 105  * This product includes cryptographic software written by Eric Young
 106  * (eay@cryptsoft.com).  This product includes software written by Tim
 107  * Hudson (tjh@cryptsoft.com).
 108  *
 109  */
 110
 111 #include "internal/cryptlib.h"
 112 #include "bn_lcl.h"
 113
 114 static BIGNUM *euclid(BIGNUM *a, BIGNUM *b);
 115
 116 int BN_gcd(BIGNUM *r, const BIGNUM *in_a, const BIGNUM *in_b, BN_CTX *ctx)
 117 {
 118     BIGNUM *a, *b, *t;
 119     int ret = 0;
 120
 121     bn_check_top(in_a);
 122     bn_check_top(in_b);
 123
 124     BN_CTX_start(ctx);
 125     a = BN_CTX_get(ctx);
 126     b = BN_CTX_get(ctx);
 127     if (a == NULL || b == NULL)
 128         goto err;
 129
 130     if (BN_copy(a, in_a) == NULL)
 131         goto err;
 132     if (BN_copy(b, in_b) == NULL)
 133         goto err;
 134     a->neg = 0;
 135     b->neg = 0;
 136
 137     if (BN_cmp(a, b) < 0) {
 138         t = a;
 139         a = b;
 140         b = t;
 141     }
 142     t = euclid(a, b);
 143     if (t == NULL)
 144         goto err;
 145
 146     if (BN_copy(r, t) == NULL)
 147         goto err;
 148     ret = 1;
 149  err:
 150     BN_CTX_end(ctx);
 151     bn_check_top(r);
 152     return (ret);
 153 }
 154
 155 static BIGNUM *euclid(BIGNUM *a, BIGNUM *b)
 156 {
 157     BIGNUM *t;
 158     int shifts = 0;
 159
 160     bn_check_top(a);
 161     bn_check_top(b);
 162
 163     /* 0 <= b <= a */
 164     while (!BN_is_zero(b)) {
 165         /* 0 < b <= a */
 166
 167         if (BN_is_odd(a)) {
 168             if (BN_is_odd(b)) {
 169                 if (!BN_sub(a, a, b))
 170                     goto err;
 171                 if (!BN_rshift1(a, a))
 172                     goto err;
 173                 if (BN_cmp(a, b) < 0) {
 174                     t = a;
 175                     a = b;
 176                     b = t;
 177                 }
 178             } else {            /* a odd - b even */
 179
 180                 if (!BN_rshift1(b, b))
 181                     goto err;
 182                 if (BN_cmp(a, b) < 0) {
 183                     t = a;
 184                     a = b;
 185                     b = t;
 186                 }
 187             }
 188         } else {                /* a is even */
 189
 190             if (BN_is_odd(b)) {
 191                 if (!BN_rshift1(a, a))
 192                     goto err;
 193                 if (BN_cmp(a, b) < 0) {
 194                     t = a;
 195                     a = b;
 196                     b = t;
 197                 }
 198             } else {            /* a even - b even */
 199
 200                 if (!BN_rshift1(a, a))
 201                     goto err;
 202                 if (!BN_rshift1(b, b))
 203                     goto err;
 204                 shifts++;
 205             }
 206         }
 207         /* 0 <= b <= a */
 208     }
 209
 210     if (shifts) {
 211         if (!BN_lshift(a, a, shifts))
 212             goto err;
 213     }
 214     bn_check_top(a);
 215     return (a);
 216  err:
 217     return (NULL);
 218 }
 219
 220 /* solves ax == 1 (mod n) */
 221 static BIGNUM *BN_mod_inverse_no_branch(BIGNUM *in,
 222                                         const BIGNUM *a, const BIGNUM *n,
 223                                         BN_CTX *ctx);
 224
 225 BIGNUM *BN_mod_inverse(BIGNUM *in,
 226                        const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
 227 {
 228     BIGNUM *rv;
 229     int noinv;
 230     rv = int_bn_mod_inverse(in, a, n, ctx, &noinv);
 231     if (noinv)
 232         BNerr(BN_F_BN_MOD_INVERSE, BN_R_NO_INVERSE);
 233     return rv;
 234 }
 235
 236 BIGNUM *int_bn_mod_inverse(BIGNUM *in,
 237                            const BIGNUM *a, const BIGNUM *n, BN_CTX *ctx,
 238                            int *pnoinv)
 239 {
 240     BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
 241     BIGNUM *ret = NULL;
 242     int sign;
 243
 244     if (pnoinv)
 245         *pnoinv = 0;
 246
 247     if ((BN_get_flags(a, BN_FLG_CONSTTIME) != 0)
 248         || (BN_get_flags(n, BN_FLG_CONSTTIME) != 0)) {
 249         return BN_mod_inverse_no_branch(in, a, n, ctx);
 250     }
 251
 252     bn_check_top(a);
 253     bn_check_top(n);
 254
 255     BN_CTX_start(ctx);
 256     A = BN_CTX_get(ctx);
 257     B = BN_CTX_get(ctx);
 258     X = BN_CTX_get(ctx);
 259     D = BN_CTX_get(ctx);
 260     M = BN_CTX_get(ctx);
 261     Y = BN_CTX_get(ctx);
 262     T = BN_CTX_get(ctx);
 263     if (T == NULL)
 264         goto err;
 265
 266     if (in == NULL)
 267         R = BN_new();
 268     else
 269         R = in;
 270     if (R == NULL)
 271         goto err;
 272
 273     BN_one(X);
 274     BN_zero(Y);
 275     if (BN_copy(B, a) == NULL)
 276         goto err;
 277     if (BN_copy(A, n) == NULL)
 278         goto err;
 279     A->neg = 0;
 280     if (B->neg || (BN_ucmp(B, A) >= 0)) {
 281         if (!BN_nnmod(B, B, A, ctx))
 282             goto err;
 283     }
 284     sign = -1;
 285     /*-
 286      * From  B = a mod |n|,  A = |n|  it follows that
 287      *
 288      *      0 <= B < A,
 289      *     -sign*X*a  ==  B   (mod |n|),
 290      *      sign*Y*a  ==  A   (mod |n|).
 291      */
 292
 293     if (BN_is_odd(n) && (BN_num_bits(n) <= 2048)) {
 294         /*
 295          * Binary inversion algorithm; requires odd modulus. This is faster
 296          * than the general algorithm if the modulus is sufficiently small
 297          * (about 400 .. 500 bits on 32-bit systems, but much more on 64-bit
 298          * systems)
 299          */
 300         int shift;
 301
 302         while (!BN_is_zero(B)) {
 303             /*-
 304              *      0 < B < |n|,
 305              *      0 < A <= |n|,
 306              * (1) -sign*X*a  ==  B   (mod |n|),
 307              * (2)  sign*Y*a  ==  A   (mod |n|)
 308              */
 309
 310             /*
 311              * Now divide B by the maximum possible power of two in the
 312              * integers, and divide X by the same value mod |n|. When we're
 313              * done, (1) still holds.
 314              */
 315             shift = 0;
 316             while (!BN_is_bit_set(B, shift)) { /* note that 0 < B */
 317                 shift++;
 318
 319                 if (BN_is_odd(X)) {
 320                     if (!BN_uadd(X, X, n))
 321                         goto err;
 322                 }
 323                 /*
 324                  * now X is even, so we can easily divide it by two
 325                  */
 326                 if (!BN_rshift1(X, X))
 327                     goto err;
 328             }
 329             if (shift > 0) {
 330                 if (!BN_rshift(B, B, shift))
 331                     goto err;
 332             }
 333
 334             /*
 335              * Same for A and Y.  Afterwards, (2) still holds.
 336              */
 337             shift = 0;
 338             while (!BN_is_bit_set(A, shift)) { /* note that 0 < A */
 339                 shift++;
 340
 341                 if (BN_is_odd(Y)) {
 342                     if (!BN_uadd(Y, Y, n))
 343                         goto err;
 344                 }
 345                 /* now Y is even */
 346                 if (!BN_rshift1(Y, Y))
 347                     goto err;
 348             }
 349             if (shift > 0) {
 350                 if (!BN_rshift(A, A, shift))
 351                     goto err;
 352             }
 353
 354             /*-
 355              * We still have (1) and (2).
 356              * Both  A  and  B  are odd.
 357              * The following computations ensure that
 358              *
 359              *     0 <= B < |n|,
 360              *      0 < A < |n|,
 361              * (1) -sign*X*a  ==  B   (mod |n|),
 362              * (2)  sign*Y*a  ==  A   (mod |n|),
 363              *
 364              * and that either  A  or  B  is even in the next iteration.
 365              */
 366             if (BN_ucmp(B, A) >= 0) {
 367                 /* -sign*(X + Y)*a == B - A  (mod |n|) */
 368                 if (!BN_uadd(X, X, Y))
 369                     goto err;
 370                 /*
 371                  * NB: we could use BN_mod_add_quick(X, X, Y, n), but that
 372                  * actually makes the algorithm slower
 373                  */
 374                 if (!BN_usub(B, B, A))
 375                     goto err;
 376             } else {
 377                 /*  sign*(X + Y)*a == A - B  (mod |n|) */
 378                 if (!BN_uadd(Y, Y, X))
 379                     goto err;
 380                 /*
 381                  * as above, BN_mod_add_quick(Y, Y, X, n) would slow things
 382                  * down
 383                  */
 384                 if (!BN_usub(A, A, B))
 385                     goto err;
 386             }
 387         }
 388     } else {
 389         /* general inversion algorithm */
 390
 391         while (!BN_is_zero(B)) {
 392             BIGNUM *tmp;
 393
 394             /*-
 395              *      0 < B < A,
 396              * (*) -sign*X*a  ==  B   (mod |n|),
 397              *      sign*Y*a  ==  A   (mod |n|)
 398              */
 399
 400             /* (D, M) := (A/B, A%B) ... */
 401             if (BN_num_bits(A) == BN_num_bits(B)) {
 402                 if (!BN_one(D))
 403                     goto err;
 404                 if (!BN_sub(M, A, B))
 405                     goto err;
 406             } else if (BN_num_bits(A) == BN_num_bits(B) + 1) {
 407                 /* A/B is 1, 2, or 3 */
 408                 if (!BN_lshift1(T, B))
 409                     goto err;
 410                 if (BN_ucmp(A, T) < 0) {
 411                     /* A < 2*B, so D=1 */
 412                     if (!BN_one(D))
 413                         goto err;
 414                     if (!BN_sub(M, A, B))
 415                         goto err;
 416                 } else {
 417                     /* A >= 2*B, so D=2 or D=3 */
 418                     if (!BN_sub(M, A, T))
 419                         goto err;
 420                     if (!BN_add(D, T, B))
 421                         goto err; /* use D (:= 3*B) as temp */
 422                     if (BN_ucmp(A, D) < 0) {
 423                         /* A < 3*B, so D=2 */
 424                         if (!BN_set_word(D, 2))
 425                             goto err;
 426                         /*
 427                          * M (= A - 2*B) already has the correct value
 428                          */
 429                     } else {
 430                         /* only D=3 remains */
 431                         if (!BN_set_word(D, 3))
 432                             goto err;
 433                         /*
 434                          * currently M = A - 2*B, but we need M = A - 3*B
 435                          */
 436                         if (!BN_sub(M, M, B))
 437                             goto err;
 438                     }
 439                 }
 440             } else {
 441                 if (!BN_div(D, M, A, B, ctx))
 442                     goto err;
 443             }
 444
 445             /*-
 446              * Now
 447              *      A = D*B + M;
 448              * thus we have
 449              * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
 450              */
 451
 452             tmp = A;            /* keep the BIGNUM object, the value does not
 453                                  * matter */
 454
 455             /* (A, B) := (B, A mod B) ... */
 456             A = B;
 457             B = M;
 458             /* ... so we have  0 <= B < A  again */
 459
 460             /*-
 461              * Since the former  M  is now  B  and the former  B  is now  A,
 462              * (**) translates into
 463              *       sign*Y*a  ==  D*A + B    (mod |n|),
 464              * i.e.
 465              *       sign*Y*a - D*A  ==  B    (mod |n|).
 466              * Similarly, (*) translates into
 467              *      -sign*X*a  ==  A          (mod |n|).
 468              *
 469              * Thus,
 470              *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
 471              * i.e.
 472              *        sign*(Y + D*X)*a  ==  B  (mod |n|).
 473              *
 474              * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
 475              *      -sign*X*a  ==  B   (mod |n|),
 476              *       sign*Y*a  ==  A   (mod |n|).
 477              * Note that  X  and  Y  stay non-negative all the time.
 478              */
 479
 480             /*
 481              * most of the time D is very small, so we can optimize tmp :=
 482              * D*X+Y
 483              */
 484             if (BN_is_one(D)) {
 485                 if (!BN_add(tmp, X, Y))
 486                     goto err;
 487             } else {
 488                 if (BN_is_word(D, 2)) {
 489                     if (!BN_lshift1(tmp, X))
 490                         goto err;
 491                 } else if (BN_is_word(D, 4)) {
 492                     if (!BN_lshift(tmp, X, 2))
 493                         goto err;
 494                 } else if (D->top == 1) {
 495                     if (!BN_copy(tmp, X))
 496                         goto err;
 497                     if (!BN_mul_word(tmp, D->d[0]))
 498                         goto err;
 499                 } else {
 500                     if (!BN_mul(tmp, D, X, ctx))
 501                         goto err;
 502                 }
 503                 if (!BN_add(tmp, tmp, Y))
 504                     goto err;
 505             }
 506
 507             M = Y;              /* keep the BIGNUM object, the value does not
 508                                  * matter */
 509             Y = X;
 510             X = tmp;
 511             sign = -sign;
 512         }
 513     }
 514
 515     /*-
 516      * The while loop (Euclid's algorithm) ends when
 517      *      A == gcd(a,n);
 518      * we have
 519      *       sign*Y*a  ==  A  (mod |n|),
 520      * where  Y  is non-negative.
 521      */
 522
 523     if (sign < 0) {
 524         if (!BN_sub(Y, n, Y))
 525             goto err;
 526     }
 527     /* Now  Y*a  ==  A  (mod |n|).  */
 528
 529     if (BN_is_one(A)) {
 530         /* Y*a == 1  (mod |n|) */
 531         if (!Y->neg && BN_ucmp(Y, n) < 0) {
 532             if (!BN_copy(R, Y))
 533                 goto err;
 534         } else {
 535             if (!BN_nnmod(R, Y, n, ctx))
 536                 goto err;
 537         }
 538     } else {
 539         if (pnoinv)
 540             *pnoinv = 1;
 541         goto err;
 542     }
 543     ret = R;
 544  err:
 545     if ((ret == NULL) && (in == NULL))
 546         BN_free(R);
 547     BN_CTX_end(ctx);
 548     bn_check_top(ret);
 549     return (ret);
 550 }
 551
 552 /*
 553  * BN_mod_inverse_no_branch is a special version of BN_mod_inverse. It does
 554  * not contain branches that may leak sensitive information.
 555  */
 556 static BIGNUM *BN_mod_inverse_no_branch(BIGNUM *in,
 557                                         const BIGNUM *a, const BIGNUM *n,
 558                                         BN_CTX *ctx)
 559 {
 560     BIGNUM *A, *B, *X, *Y, *M, *D, *T, *R = NULL;
 561     BIGNUM *ret = NULL;
 562     int sign;
 563
 564     bn_check_top(a);
 565     bn_check_top(n);
 566
 567     BN_CTX_start(ctx);
 568     A = BN_CTX_get(ctx);
 569     B = BN_CTX_get(ctx);
 570     X = BN_CTX_get(ctx);
 571     D = BN_CTX_get(ctx);
 572     M = BN_CTX_get(ctx);
 573     Y = BN_CTX_get(ctx);
 574     T = BN_CTX_get(ctx);
 575     if (T == NULL)
 576         goto err;
 577
 578     if (in == NULL)
 579         R = BN_new();
 580     else
 581         R = in;
 582     if (R == NULL)
 583         goto err;
 584
 585     BN_one(X);
 586     BN_zero(Y);
 587     if (BN_copy(B, a) == NULL)
 588         goto err;
 589     if (BN_copy(A, n) == NULL)
 590         goto err;
 591     A->neg = 0;
 592
 593     if (B->neg || (BN_ucmp(B, A) >= 0)) {
 594         /*
 595          * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
 596          * BN_div_no_branch will be called eventually.
 597          */
 598          {
 599             BIGNUM local_B;
 600             bn_init(&local_B);
 601             BN_with_flags(&local_B, B, BN_FLG_CONSTTIME);
 602             if (!BN_nnmod(B, &local_B, A, ctx))
 603                 goto err;
 604             /* Ensure local_B goes out of scope before any further use of B */
 605         }
 606     }
 607     sign = -1;
 608     /*-
 609      * From  B = a mod |n|,  A = |n|  it follows that
 610      *
 611      *      0 <= B < A,
 612      *     -sign*X*a  ==  B   (mod |n|),
 613      *      sign*Y*a  ==  A   (mod |n|).
 614      */
 615
 616     while (!BN_is_zero(B)) {
 617         BIGNUM *tmp;
 618
 619         /*-
 620          *      0 < B < A,
 621          * (*) -sign*X*a  ==  B   (mod |n|),
 622          *      sign*Y*a  ==  A   (mod |n|)
 623          */
 624
 625         /*
 626          * Turn BN_FLG_CONSTTIME flag on, so that when BN_div is invoked,
 627          * BN_div_no_branch will be called eventually.
 628          */
 629         {
 630             BIGNUM local_A;
 631             bn_init(&local_A);
 632             BN_with_flags(&local_A, A, BN_FLG_CONSTTIME);
 633
 634             /* (D, M) := (A/B, A%B) ... */
 635             if (!BN_div(D, M, &local_A, B, ctx))
 636                 goto err;
 637             /* Ensure local_A goes out of scope before any further use of A */
 638         }
 639
 640         /*-
 641          * Now
 642          *      A = D*B + M;
 643          * thus we have
 644          * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
 645          */
 646
 647         tmp = A;                /* keep the BIGNUM object, the value does not
 648                                  * matter */
 649
 650         /* (A, B) := (B, A mod B) ... */
 651         A = B;
 652         B = M;
 653         /* ... so we have  0 <= B < A  again */
 654
 655         /*-
 656          * Since the former  M  is now  B  and the former  B  is now  A,
 657          * (**) translates into
 658          *       sign*Y*a  ==  D*A + B    (mod |n|),
 659          * i.e.
 660          *       sign*Y*a - D*A  ==  B    (mod |n|).
 661          * Similarly, (*) translates into
 662          *      -sign*X*a  ==  A          (mod |n|).
 663          *
 664          * Thus,
 665          *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
 666          * i.e.
 667          *        sign*(Y + D*X)*a  ==  B  (mod |n|).
 668          *
 669          * So if we set  (X, Y, sign) := (Y + D*X, X, -sign),  we arrive back at
 670          *      -sign*X*a  ==  B   (mod |n|),
 671          *       sign*Y*a  ==  A   (mod |n|).
 672          * Note that  X  and  Y  stay non-negative all the time.
 673          */
 674
 675         if (!BN_mul(tmp, D, X, ctx))
 676             goto err;
 677         if (!BN_add(tmp, tmp, Y))
 678             goto err;
 679
 680         M = Y;                  /* keep the BIGNUM object, the value does not
 681                                  * matter */
 682         Y = X;
 683         X = tmp;
 684         sign = -sign;
 685     }
 686
 687     /*-
 688      * The while loop (Euclid's algorithm) ends when
 689      *      A == gcd(a,n);
 690      * we have
 691      *       sign*Y*a  ==  A  (mod |n|),
 692      * where  Y  is non-negative.
 693      */
 694
 695     if (sign < 0) {
 696         if (!BN_sub(Y, n, Y))
 697             goto err;
 698     }
 699     /* Now  Y*a  ==  A  (mod |n|).  */
 700
 701     if (BN_is_one(A)) {
 702         /* Y*a == 1  (mod |n|) */
 703         if (!Y->neg && BN_ucmp(Y, n) < 0) {
 704             if (!BN_copy(R, Y))
 705                 goto err;
 706         } else {
 707             if (!BN_nnmod(R, Y, n, ctx))
 708                 goto err;
 709         }
 710     } else {
 711         BNerr(BN_F_BN_MOD_INVERSE_NO_BRANCH, BN_R_NO_INVERSE);
 712         goto err;
 713     }
 714     ret = R;
 715  err:
 716     if ((ret == NULL) && (in == NULL))
 717         BN_free(R);
 718     BN_CTX_end(ctx);
 719     bn_check_top(ret);
 720     return (ret);
 721 }