cscrypt/bn_exp.c

   1 #include "bn.h"
   2
   3 #ifndef WITH_LIBCRYPTO
   4 //FIXME Not checked on threadsafety yet; after checking please remove this line
   5 /* crypto/bn/bn_exp.c */
   6 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
   7  * All rights reserved.
   8  *
   9  * This package is an SSL implementation written
  10  * by Eric Young (eay@cryptsoft.com).
  11  * The implementation was written so as to conform with Netscapes SSL.
  12  *
  13  * This library is free for commercial and non-commercial use as long as
  14  * the following conditions are aheared to.  The following conditions
  15  * apply to all code found in this distribution, be it the RC4, RSA,
  16  * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
  17  * included with this distribution is covered by the same copyright terms
  18  * except that the holder is Tim Hudson (tjh@cryptsoft.com).
  19  *
  20  * Copyright remains Eric Young's, and as such any Copyright notices in
  21  * the code are not to be removed.
  22  * If this package is used in a product, Eric Young should be given attribution
  23  * as the author of the parts of the library used.
  24  * This can be in the form of a textual message at program startup or
  25  * in documentation (online or textual) provided with the package.
  26  *
  27  * Redistribution and use in source and binary forms, with or without
  28  * modification, are permitted provided that the following conditions
  29  * are met:
  30  * 1. Redistributions of source code must retain the copyright
  31  *    notice, this list of conditions and the following disclaimer.
  32  * 2. Redistributions in binary form must reproduce the above copyright
  33  *    notice, this list of conditions and the following disclaimer in the
  34  *    documentation and/or other materials provided with the distribution.
  35  * 3. All advertising materials mentioning features or use of this software
  36  *    must display the following acknowledgement:
  37  *    "This product includes cryptographic software written by
  38  *     Eric Young (eay@cryptsoft.com)"
  39  *    The word 'cryptographic' can be left out if the rouines from the library
  40  *    being used are not cryptographic related :-).
  41  * 4. If you include any Windows specific code (or a derivative thereof) from
  42  *    the apps directory (application code) you must include an acknowledgement:
  43  *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
  44  *
  45  * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND
  46  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  47  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  48  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  49  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  50  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  51  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  52  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  53  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  54  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  55  * SUCH DAMAGE.
  56  *
  57  * The license and distribution terms for any publically available version or
  58  * derivative of this code cannot be changed.  i.e. this code cannot simply be
  59  * copied and put under another distribution license
  60  * [including the GNU Public License.]
  61  */
  62 /* ====================================================================
  63  * Copyright (c) 1998-2000 The OpenSSL Project.  All rights reserved.
  64  *
  65  * Redistribution and use in source and binary forms, with or without
  66  * modification, are permitted provided that the following conditions
  67  * are met:
  68  *
  69  * 1. Redistributions of source code must retain the above copyright
  70  *    notice, this list of conditions and the following disclaimer.
  71  *
  72  * 2. Redistributions in binary form must reproduce the above copyright
  73  *    notice, this list of conditions and the following disclaimer in
  74  *    the documentation and/or other materials provided with the
  75  *    distribution.
  76  *
  77  * 3. All advertising materials mentioning features or use of this
  78  *    software must display the following acknowledgment:
  79  *    "This product includes software developed by the OpenSSL Project
  80  *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
  81  *
  82  * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
  83  *    endorse or promote products derived from this software without
  84  *    prior written permission. For written permission, please contact
  85  *    openssl-core@openssl.org.
  86  *
  87  * 5. Products derived from this software may not be called "OpenSSL"
  88  *    nor may "OpenSSL" appear in their names without prior written
  89  *    permission of the OpenSSL Project.
  90  *
  91  * 6. Redistributions of any form whatsoever must retain the following
  92  *    acknowledgment:
  93  *    "This product includes software developed by the OpenSSL Project
  94  *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
  95  *
  96  * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
  97  * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  98  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  99  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
 100  * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 101  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 102  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 103  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 104  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 105  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 106  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
 107  * OF THE POSSIBILITY OF SUCH DAMAGE.
 108  * ====================================================================
 109  *
 110  * This product includes cryptographic software written by Eric Young
 111  * (eay@cryptsoft.com).  This product includes software written by Tim
 112  * Hudson (tjh@cryptsoft.com).
 113  *
 114  */
 115
 116
 117 #include <stdio.h>
 118 #include "bn_lcl.h"
 119
 120 #define TABLE_SIZE  32
 121
 122 /* slow but works */
 123 int BN_mod_mul(BIGNUM *ret, BIGNUM *a, BIGNUM *b, const BIGNUM *m, BN_CTX *ctx)
 124 {
 125         BIGNUM *t;
 126         int r = 0;
 127
 128         bn_check_top(a);
 129         bn_check_top(b);
 130         bn_check_top(m);
 131
 132         BN_CTX_start(ctx);
 133         if((t = BN_CTX_get(ctx)) == NULL) { goto err; }
 134         if(a == b)
 135         {
 136                 if(!BN_sqr(t, a, ctx)) { goto err; }
 137         }
 138         else
 139         {
 140                 if(!BN_mul(t, a, b, ctx)) { goto err; }
 141         }
 142         if(!BN_mod(ret, t, m, ctx)) { goto err; }
 143         r = 1;
 144 err:
 145         BN_CTX_end(ctx);
 146         return (r);
 147 }
 148
 149
 150 int BN_mod_exp(BIGNUM *r, BIGNUM *a, const BIGNUM *p, const BIGNUM *m,
 151                            BN_CTX *ctx)
 152 {
 153         int ret;
 154
 155         bn_check_top(a);
 156         bn_check_top(p);
 157         bn_check_top(m);
 158
 159         ret = BN_mod_exp_simple(r, a, p, m, ctx);
 160
 161         return (ret);
 162 }
 163
 164
 165
 166
 167
 168 /* The old fallback, simple version :-) */
 169 int BN_mod_exp_simple(BIGNUM *r, BIGNUM *a, const BIGNUM *p, const BIGNUM *m,
 170                                           BN_CTX *ctx)
 171 {
 172         int i, j = 0, bits, ret = 0, wstart = 0, wend = 0, window, wvalue = 0, ts = 0;
 173         int start = 1;
 174         BIGNUM *d;
 175         BIGNUM val[TABLE_SIZE];
 176
 177         bits = BN_num_bits(p);
 178
 179         if(bits == 0)
 180         {
 181                 BN_one(r);
 182                 return (1);
 183         }
 184
 185         BN_CTX_start(ctx);
 186         if((d = BN_CTX_get(ctx)) == NULL) { goto err; }
 187
 188         BN_init(&(val[0]));
 189         ts = 1;
 190         if(!BN_mod(&(val[0]), a, m, ctx)) { goto err; }     /* 1 */
 191
 192         window = BN_window_bits_for_exponent_size(bits);
 193         if(window > 1)
 194         {
 195                 if(!BN_mod_mul(d, &(val[0]), &(val[0]), m, ctx))
 196                         { goto err; }               /* 2 */
 197                 j = 1 << (window - 1);
 198                 for(i = 1; i < j; i++)
 199                 {
 200                         BN_init(&(val[i]));
 201                         if(!BN_mod_mul(&(val[i]), &(val[i - 1]), d, m, ctx))
 202                                 { goto err; }
 203                 }
 204                 ts = i;
 205         }
 206
 207         start = 1;    /* This is used to avoid multiplication etc
 208              * when there is only the value '1' in the
 209              * buffer. */
 210         wstart = bits - 1; /* The top bit of the window */
 211
 212         if(!BN_one(r)) { goto err; }
 213
 214         for(;;)
 215         {
 216                 if(BN_is_bit_set(p, wstart) == 0)
 217                 {
 218                         if(!start)
 219                                 if(!BN_mod_mul(r, r, r, m, ctx))
 220                                         { goto err; }
 221                         if(wstart == 0) { break; }
 222                         wstart--;
 223                         continue;
 224                 }
 225                 /* We now have wstart on a 'set' bit, we now need to work out
 226                  * how bit a window to do.  To do this we need to scan
 227                  * forward until the last set bit before the end of the
 228                  * window */
 229                 j = wstart;
 230                 wvalue = 1;
 231                 wend = 0;
 232                 for(i = 1; i < window; i++)
 233                 {
 234                         if(wstart - i < 0) { break; }
 235                         if(BN_is_bit_set(p, wstart - i))
 236                         {
 237                                 wvalue <<= (i - wend);
 238                                 wvalue |= 1;
 239                                 wend = i;
 240                         }
 241                 }
 242
 243                 /* wend is the size of the current window */
 244                 j = wend + 1;
 245                 /* add the 'bytes above' */
 246                 if(!start)
 247                         for(i = 0; i < j; i++)
 248                         {
 249                                 if(!BN_mod_mul(r, r, r, m, ctx))
 250                                         { goto err; }
 251                         }
 252
 253                 /* wvalue will be an odd number < 2^window */
 254                 if(!BN_mod_mul(r, r, &(val[wvalue >> 1]), m, ctx))
 255                         { goto err; }
 256
 257                 /* move the 'window' down further */
 258                 wstart -= wend + 1;
 259                 wvalue = 0;
 260                 start = 0;
 261                 if(wstart < 0) { break; }
 262         }
 263         ret = 1;
 264 err:
 265         BN_CTX_end(ctx);
 266         for(i = 0; i < ts; i++)
 267                 { BN_clear_free(&(val[i])); }
 268         return (ret);
 269 }
 270
 271 int
 272 BN_nnmod(BIGNUM *r, const BIGNUM *m, const BIGNUM *d, BN_CTX *ctx)
 273 {
 274         /* like BN_mod, but returns non-negative remainder
 275          * (i.e., 0 <= r < |d|  always holds)
 276          */
 277
 278         if (!(BN_mod(r, m, d, ctx)))
 279                 return 0;
 280         if (!r->neg)
 281                 return 1;
 282         /* now   -|d| < r < 0, so we have to set  r : = r + |d| */
 283         return (d->neg ? BN_sub : BN_add)(r, r, d);
 284 }
 285
 286 /* solves ax == 1 (mod n) */
 287 BIGNUM *
 288 BN_mod_inverse(BIGNUM *ret, BIGNUM *a, const BIGNUM *n, BN_CTX *ctx)
 289 {
 290         BIGNUM *A, *B, *X, *Y, *M, *D, *T = NULL;
 291         int sign;
 292
 293         bn_check_top(a);
 294         bn_check_top(n);
 295
 296         BN_CTX_start(ctx);
 297         A = BN_CTX_get(ctx);
 298         B = BN_CTX_get(ctx);
 299         X = BN_CTX_get(ctx);
 300         D = BN_CTX_get(ctx);
 301         M = BN_CTX_get(ctx);
 302         Y = BN_CTX_get(ctx);
 303         T = BN_CTX_get(ctx);
 304         if (T == NULL) goto err;
 305
 306         if (ret == NULL) goto err;
 307
 308         BN_one(X);
 309         BN_zero(Y);
 310         if (BN_copy(B, a) == NULL) goto err;
 311         if (BN_copy(A, n) == NULL) goto err;
 312         A->neg = 0;
 313         if (B->neg || (BN_ucmp(B, A) >= 0)) {
 314                 if (!BN_nnmod(B, B, A, ctx)) goto err;
 315         }
 316         sign = -1;
 317         /* From  B = a mod |n|, A = |n|  it follows that
 318          *
 319          *      0 <= B < A,
 320          *     -sign*X*a  ==  B   (mod |n|),
 321          *      sign*Y*a  ==  A   (mod |n|).
 322          */
 323
 324         if (BN_is_odd(n) && (BN_num_bits(n) <= (BN_BITS <= 32 ? 450 : 2048))) {
 325                 /* Binary inversion algorithm; requires odd modulus.
 326                  * This is faster than the general algorithm if the modulus
 327                  * is sufficiently small (about 400 .. 500 bits on 32-bit
 328                  * sytems, but much more on 64-bit systems)
 329                  */
 330                 int shift;
 331
 332                 while (!BN_is_zero(B)) {
 333                         /*
 334                          *      0 < B < |n|,
 335                          *      0 < A <= |n|,
 336                          * (1) -sign*X*a  ==  B   (mod |n|),
 337                          * (2)  sign*Y*a  ==  A   (mod |n|)
 338                          */
 339
 340                         /* Now divide  B  by the maximum possible power of two in the integers,
 341                          * and divide  X  by the same value mod |n|.
 342                          * When we're done, (1) still holds.
 343                          */
 344                         shift = 0;
 345                         while (!BN_is_bit_set(B, shift)) /* note that 0 < B */ {
 346                                 shift++;
 347
 348                                 if (BN_is_odd(X)) {
 349                                         if (!BN_uadd(X, X, n)) goto err;
 350                                 }
 351                                 /* now X is even, so we can easily divide it by two */
 352                                 if (!BN_rshift1(X, X)) goto err;
 353                         }
 354                         if (shift > 0) {
 355                                 if (!BN_rshift(B, B, shift)) goto err;
 356                         }
 357
 358
 359                         /* Same for  A  and  Y.  Afterwards, (2) still holds. */
 360                         shift = 0;
 361                         while (!BN_is_bit_set(A, shift)) /* note that 0 < A */ {
 362                                 shift++;
 363
 364                                 if (BN_is_odd(Y)) {
 365                                         if (!BN_uadd(Y, Y, n)) goto err;
 366                                         }
 367                                 /* now Y is even */
 368                                 if (!BN_rshift1(Y, Y)) goto err;
 369                         }
 370                         if (shift > 0) {
 371                                 if (!BN_rshift(A, A, shift)) goto err;
 372                         }
 373
 374                         /* We still have (1) and (2).
 375                          * Both  A  and  B  are odd.
 376                          * The following computations ensure that
 377                          *
 378                          *     0 <= B < |n|,
 379                          *      0 < A < |n|,
 380                          * (1) -sign*X*a  ==  B   (mod |n|),
 381                          * (2)  sign*Y*a  ==  A   (mod |n|),
 382                          *
 383                          * and that either  A  or  B  is even in the next iteration.
 384                          */
 385                         if (BN_ucmp(B, A) >= 0) {
 386                                 /* -sign*(X + Y)*a == B - A  (mod |n|) */
 387                                 if (!BN_uadd(X, X, Y)) goto err;
 388                                 /* NB: we could use BN_mod_add_quick(X, X, Y, n), but that
 389                                  * actually makes the algorithm slower
 390                                  */
 391                                 if (!BN_usub(B, B, A)) goto err;
 392                         } else {
 393                                 /*  sign*(X + Y)*a == A - B  (mod |n|) */
 394                                 if (!BN_uadd(Y, Y, X)) goto err;
 395                                 /* as above, BN_mod_add_quick(Y, Y, X, n) would slow things down */
 396                                 if (!BN_usub(A, A, B)) goto err;
 397                         }
 398                 }
 399         } else {
 400                 /* general inversion algorithm */
 401
 402                 while (!BN_is_zero(B)) {
 403                         BIGNUM *tmp;
 404
 405                         /*
 406                          *      0 < B < A,
 407                          * (*) -sign*X*a  ==  B   (mod |n|),
 408                          *      sign*Y*a  ==  A   (mod |n|)
 409                          */
 410
 411                         /* (D, M) : = (A/B, A%B) ... */
 412                         if (BN_num_bits(A) == BN_num_bits(B)) {
 413                                 if (!BN_one(D)) goto err;
 414                                 if (!BN_sub(M, A, B)) goto err;
 415                         } else if (BN_num_bits(A) == BN_num_bits(B) + 1) {
 416                                 /* A/B is 1, 2, or 3 */
 417                                 if (!BN_lshift1(T, B)) goto err;
 418                                 if (BN_ucmp(A, T) < 0) {
 419                                         /* A < 2*B, so D = 1 */
 420                                         if (!BN_one(D)) goto err;
 421                                         if (!BN_sub(M, A, B)) goto err;
 422                                 } else {
 423                                         /* A >= 2*B, so D = 2 or D = 3 */
 424                                         if (!BN_sub(M, A, T)) goto err;
 425                                         if (!BN_add(D, T, B)) goto err;
 426                                         /* use D ( := 3 * B) as temp */
 427                                         if (BN_ucmp(A, D) < 0) {
 428                                                 /* A < 3*B, so D = 2 */
 429                                                 if (!BN_set_word(D, 2)) goto err;
 430                                                 /* M ( = A - 2*B) already has the correct value */
 431                                         } else {
 432                                                 /* only D = 3 remains */
 433                                                 if (!BN_set_word(D, 3)) goto err;
 434                                                 /* currently  M = A - 2 * B,
 435                                                  * but we need  M = A - 3 * B
 436                                                  */
 437                                                 if (!BN_sub(M, M, B)) goto err;
 438                                         }
 439                                 }
 440                         } else {
 441                                 if (!BN_div(D, M, A, B, ctx)) goto err;
 442                         }
 443
 444                         /* Now
 445                          *      A = D*B + M;
 446                          * thus we have
 447                          * (**)  sign*Y*a  ==  D*B + M   (mod |n|).
 448                          */
 449
 450                         tmp = A; /* keep the BIGNUM object, the value does not matter */
 451
 452                         /* (A, B) : = (B, A mod B) ... */
 453                         A = B;
 454                         B = M;
 455                         /* ... so we have  0 <= B < A  again */
 456
 457                         /* Since the former  M  is now  B  and the former  B  is now  A,
 458                          * (**) translates into
 459                          *       sign*Y*a  ==  D*A + B    (mod |n|),
 460                          * i.e.
 461                          *       sign*Y*a - D*A  ==  B    (mod |n|).
 462                          * Similarly, (*) translates into
 463                          *      -sign*X*a  ==  A          (mod |n|).
 464                          *
 465                          * Thus,
 466                          *   sign*Y*a + D*sign*X*a  ==  B  (mod |n|),
 467                          * i.e.
 468                          *        sign*(Y + D*X)*a  ==  B  (mod |n|).
 469                          *
 470                          * So if we set  (X, Y, sign) : = (Y + D*X, X, -sign), we arrive back at
 471                          *      -sign*X*a  ==  B   (mod |n|),
 472                          *       sign*Y*a  ==  A   (mod |n|).
 473                          * Note that  X  and  Y  stay non-negative all the time.
 474                          */
 475
 476                         /* most of the time D is very small, so we can optimize tmp : = D*X+Y */
 477                         if (BN_is_one(D)) {
 478                                 if (!BN_add(tmp, X, Y)) goto err;
 479                         } else {
 480                                 if (BN_is_word(D, 2)) {
 481                                         if (!BN_lshift1(tmp, X)) goto err;
 482                                 } else if (BN_is_word(D, 4)) {
 483                                         if (!BN_lshift(tmp, X, 2)) goto err;
 484                                 } else if (D->top == 1) {
 485                                         if (!BN_copy(tmp, X)) goto err;
 486                                         if (!BN_mul_word(tmp, D->d[0])) goto err;
 487                                 } else {
 488                                         if (!BN_mul(tmp, D, X, ctx)) goto err;
 489                                 }
 490                                 if (!BN_add(tmp, tmp, Y)) goto err;
 491                         }
 492
 493                         M = Y; /* keep the BIGNUM object, the value does not matter */
 494                         Y = X;
 495                         X = tmp;
 496                         sign = -sign;
 497                 }
 498         }
 499
 500         /*
 501          * The while loop (Euclid's algorithm) ends when
 502          *      A == gcd(a, n);
 503          * we have
 504          *       sign*Y*a  ==  A  (mod |n|),
 505          * where  Y  is non-negative.
 506          */
 507
 508         if (sign < 0) {
 509                 if (!BN_sub(Y, n, Y)) goto err;
 510         }
 511         /* Now  Y*a  ==  A  (mod |n|).  */
 512
 513
 514         if (BN_is_one(A)) {
 515                 /* Y*a == 1  (mod |n|) */
 516                 if (!Y->neg && BN_ucmp(Y, n) < 0) {
 517                         if (!BN_copy(ret, Y)) goto err;
 518                 } else {
 519                         if (!BN_nnmod(ret, Y, n, ctx)) goto err;
 520                 }
 521         } else {
 522                 goto err;
 523         }
 524 err:
 525         BN_CTX_end(ctx);
 526         return (ret);
 527 }
 528
 529 #endif