crypto/bn/bn_mul.c

   1 /* crypto/bn/bn_mul.c */
   2 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
   3  * All rights reserved.
   4  *
   5  * This package is an SSL implementation written
   6  * by Eric Young (eay@cryptsoft.com).
   7  * The implementation was written so as to conform with Netscapes SSL.
   8  *
   9  * This library is free for commercial and non-commercial use as long as
  10  * the following conditions are aheared to.  The following conditions
  11  * apply to all code found in this distribution, be it the RC4, RSA,
  12  * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
  13  * included with this distribution is covered by the same copyright terms
  14  * except that the holder is Tim Hudson (tjh@cryptsoft.com).
  15  *
  16  * Copyright remains Eric Young's, and as such any Copyright notices in
  17  * the code are not to be removed.
  18  * If this package is used in a product, Eric Young should be given attribution
  19  * as the author of the parts of the library used.
  20  * This can be in the form of a textual message at program startup or
  21  * in documentation (online or textual) provided with the package.
  22  *
  23  * Redistribution and use in source and binary forms, with or without
  24  * modification, are permitted provided that the following conditions
  25  * are met:
  26  * 1. Redistributions of source code must retain the copyright
  27  *    notice, this list of conditions and the following disclaimer.
  28  * 2. Redistributions in binary form must reproduce the above copyright
  29  *    notice, this list of conditions and the following disclaimer in the
  30  *    documentation and/or other materials provided with the distribution.
  31  * 3. All advertising materials mentioning features or use of this software
  32  *    must display the following acknowledgement:
  33  *    "This product includes cryptographic software written by
  34  *     Eric Young (eay@cryptsoft.com)"
  35  *    The word 'cryptographic' can be left out if the rouines from the library
  36  *    being used are not cryptographic related :-).
  37  * 4. If you include any Windows specific code (or a derivative thereof) from
  38  *    the apps directory (application code) you must include an acknowledgement:
  39  *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
  40  *
  41  * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND
  42  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  43  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  44  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  45  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  46  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  47  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  48  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  49  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  50  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  51  * SUCH DAMAGE.
  52  *
  53  * The licence and distribution terms for any publically available version or
  54  * derivative of this code cannot be changed.  i.e. this code cannot simply be
  55  * copied and put under another distribution licence
  56  * [including the GNU Public Licence.]
  57  */
  58
  59 #ifndef BN_DEBUG
  60 # undef NDEBUG                  /* avoid conflicting definitions */
  61 # define NDEBUG
  62 #endif
  63
  64 #include <assert.h>
  65 #include "internal/cryptlib.h"
  66 #include "bn_lcl.h"
  67
  68 #if defined(OPENSSL_NO_ASM) || !defined(OPENSSL_BN_ASM_PART_WORDS)
  69 /*
  70  * Here follows specialised variants of bn_add_words() and bn_sub_words().
  71  * They have the property performing operations on arrays of different sizes.
  72  * The sizes of those arrays is expressed through cl, which is the common
  73  * length ( basicall, min(len(a),len(b)) ), and dl, which is the delta
  74  * between the two lengths, calculated as len(a)-len(b). All lengths are the
  75  * number of BN_ULONGs...  For the operations that require a result array as
  76  * parameter, it must have the length cl+abs(dl). These functions should
  77  * probably end up in bn_asm.c as soon as there are assembler counterparts
  78  * for the systems that use assembler files.
  79  */
  80
  81 BN_ULONG bn_sub_part_words(BN_ULONG *r,
  82                            const BN_ULONG *a, const BN_ULONG *b,
  83                            int cl, int dl)
  84 {
  85     BN_ULONG c, t;
  86
  87     assert(cl >= 0);
  88     c = bn_sub_words(r, a, b, cl);
  89
  90     if (dl == 0)
  91         return c;
  92
  93     r += cl;
  94     a += cl;
  95     b += cl;
  96
  97     if (dl < 0) {
  98         for (;;) {
  99             t = b[0];
 100             r[0] = (0 - t - c) & BN_MASK2;
 101             if (t != 0)
 102                 c = 1;
 103             if (++dl >= 0)
 104                 break;
 105
 106             t = b[1];
 107             r[1] = (0 - t - c) & BN_MASK2;
 108             if (t != 0)
 109                 c = 1;
 110             if (++dl >= 0)
 111                 break;
 112
 113             t = b[2];
 114             r[2] = (0 - t - c) & BN_MASK2;
 115             if (t != 0)
 116                 c = 1;
 117             if (++dl >= 0)
 118                 break;
 119
 120             t = b[3];
 121             r[3] = (0 - t - c) & BN_MASK2;
 122             if (t != 0)
 123                 c = 1;
 124             if (++dl >= 0)
 125                 break;
 126
 127             b += 4;
 128             r += 4;
 129         }
 130     } else {
 131         int save_dl = dl;
 132         while (c) {
 133             t = a[0];
 134             r[0] = (t - c) & BN_MASK2;
 135             if (t != 0)
 136                 c = 0;
 137             if (--dl <= 0)
 138                 break;
 139
 140             t = a[1];
 141             r[1] = (t - c) & BN_MASK2;
 142             if (t != 0)
 143                 c = 0;
 144             if (--dl <= 0)
 145                 break;
 146
 147             t = a[2];
 148             r[2] = (t - c) & BN_MASK2;
 149             if (t != 0)
 150                 c = 0;
 151             if (--dl <= 0)
 152                 break;
 153
 154             t = a[3];
 155             r[3] = (t - c) & BN_MASK2;
 156             if (t != 0)
 157                 c = 0;
 158             if (--dl <= 0)
 159                 break;
 160
 161             save_dl = dl;
 162             a += 4;
 163             r += 4;
 164         }
 165         if (dl > 0) {
 166             if (save_dl > dl) {
 167                 switch (save_dl - dl) {
 168                 case 1:
 169                     r[1] = a[1];
 170                     if (--dl <= 0)
 171                         break;
 172                 case 2:
 173                     r[2] = a[2];
 174                     if (--dl <= 0)
 175                         break;
 176                 case 3:
 177                     r[3] = a[3];
 178                     if (--dl <= 0)
 179                         break;
 180                 }
 181                 a += 4;
 182                 r += 4;
 183             }
 184         }
 185         if (dl > 0) {
 186             for (;;) {
 187                 r[0] = a[0];
 188                 if (--dl <= 0)
 189                     break;
 190                 r[1] = a[1];
 191                 if (--dl <= 0)
 192                     break;
 193                 r[2] = a[2];
 194                 if (--dl <= 0)
 195                     break;
 196                 r[3] = a[3];
 197                 if (--dl <= 0)
 198                     break;
 199
 200                 a += 4;
 201                 r += 4;
 202             }
 203         }
 204     }
 205     return c;
 206 }
 207 #endif
 208
 209 BN_ULONG bn_add_part_words(BN_ULONG *r,
 210                            const BN_ULONG *a, const BN_ULONG *b,
 211                            int cl, int dl)
 212 {
 213     BN_ULONG c, l, t;
 214
 215     assert(cl >= 0);
 216     c = bn_add_words(r, a, b, cl);
 217
 218     if (dl == 0)
 219         return c;
 220
 221     r += cl;
 222     a += cl;
 223     b += cl;
 224
 225     if (dl < 0) {
 226         int save_dl = dl;
 227         while (c) {
 228             l = (c + b[0]) & BN_MASK2;
 229             c = (l < c);
 230             r[0] = l;
 231             if (++dl >= 0)
 232                 break;
 233
 234             l = (c + b[1]) & BN_MASK2;
 235             c = (l < c);
 236             r[1] = l;
 237             if (++dl >= 0)
 238                 break;
 239
 240             l = (c + b[2]) & BN_MASK2;
 241             c = (l < c);
 242             r[2] = l;
 243             if (++dl >= 0)
 244                 break;
 245
 246             l = (c + b[3]) & BN_MASK2;
 247             c = (l < c);
 248             r[3] = l;
 249             if (++dl >= 0)
 250                 break;
 251
 252             save_dl = dl;
 253             b += 4;
 254             r += 4;
 255         }
 256         if (dl < 0) {
 257             if (save_dl < dl) {
 258                 switch (dl - save_dl) {
 259                 case 1:
 260                     r[1] = b[1];
 261                     if (++dl >= 0)
 262                         break;
 263                 case 2:
 264                     r[2] = b[2];
 265                     if (++dl >= 0)
 266                         break;
 267                 case 3:
 268                     r[3] = b[3];
 269                     if (++dl >= 0)
 270                         break;
 271                 }
 272                 b += 4;
 273                 r += 4;
 274             }
 275         }
 276         if (dl < 0) {
 277             for (;;) {
 278                 r[0] = b[0];
 279                 if (++dl >= 0)
 280                     break;
 281                 r[1] = b[1];
 282                 if (++dl >= 0)
 283                     break;
 284                 r[2] = b[2];
 285                 if (++dl >= 0)
 286                     break;
 287                 r[3] = b[3];
 288                 if (++dl >= 0)
 289                     break;
 290
 291                 b += 4;
 292                 r += 4;
 293             }
 294         }
 295     } else {
 296         int save_dl = dl;
 297         while (c) {
 298             t = (a[0] + c) & BN_MASK2;
 299             c = (t < c);
 300             r[0] = t;
 301             if (--dl <= 0)
 302                 break;
 303
 304             t = (a[1] + c) & BN_MASK2;
 305             c = (t < c);
 306             r[1] = t;
 307             if (--dl <= 0)
 308                 break;
 309
 310             t = (a[2] + c) & BN_MASK2;
 311             c = (t < c);
 312             r[2] = t;
 313             if (--dl <= 0)
 314                 break;
 315
 316             t = (a[3] + c) & BN_MASK2;
 317             c = (t < c);
 318             r[3] = t;
 319             if (--dl <= 0)
 320                 break;
 321
 322             save_dl = dl;
 323             a += 4;
 324             r += 4;
 325         }
 326         if (dl > 0) {
 327             if (save_dl > dl) {
 328                 switch (save_dl - dl) {
 329                 case 1:
 330                     r[1] = a[1];
 331                     if (--dl <= 0)
 332                         break;
 333                 case 2:
 334                     r[2] = a[2];
 335                     if (--dl <= 0)
 336                         break;
 337                 case 3:
 338                     r[3] = a[3];
 339                     if (--dl <= 0)
 340                         break;
 341                 }
 342                 a += 4;
 343                 r += 4;
 344             }
 345         }
 346         if (dl > 0) {
 347             for (;;) {
 348                 r[0] = a[0];
 349                 if (--dl <= 0)
 350                     break;
 351                 r[1] = a[1];
 352                 if (--dl <= 0)
 353                     break;
 354                 r[2] = a[2];
 355                 if (--dl <= 0)
 356                     break;
 357                 r[3] = a[3];
 358                 if (--dl <= 0)
 359                     break;
 360
 361                 a += 4;
 362                 r += 4;
 363             }
 364         }
 365     }
 366     return c;
 367 }
 368
 369 #ifdef BN_RECURSION
 370 /*
 371  * Karatsuba recursive multiplication algorithm (cf. Knuth, The Art of
 372  * Computer Programming, Vol. 2)
 373  */
 374
 375 /*-
 376  * r is 2*n2 words in size,
 377  * a and b are both n2 words in size.
 378  * n2 must be a power of 2.
 379  * We multiply and return the result.
 380  * t must be 2*n2 words in size
 381  * We calculate
 382  * a[0]*b[0]
 383  * a[0]*b[0]+a[1]*b[1]+(a[0]-a[1])*(b[1]-b[0])
 384  * a[1]*b[1]
 385  */
 386 /* dnX may not be positive, but n2/2+dnX has to be */
 387 void bn_mul_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
 388                       int dna, int dnb, BN_ULONG *t)
 389 {
 390     int n = n2 / 2, c1, c2;
 391     int tna = n + dna, tnb = n + dnb;
 392     unsigned int neg, zero;
 393     BN_ULONG ln, lo, *p;
 394
 395 # ifdef BN_MUL_COMBA
 396 #  if 0
 397     if (n2 == 4) {
 398         bn_mul_comba4(r, a, b);
 399         return;
 400     }
 401 #  endif
 402     /*
 403      * Only call bn_mul_comba 8 if n2 == 8 and the two arrays are complete
 404      * [steve]
 405      */
 406     if (n2 == 8 && dna == 0 && dnb == 0) {
 407         bn_mul_comba8(r, a, b);
 408         return;
 409     }
 410 # endif                         /* BN_MUL_COMBA */
 411     /* Else do normal multiply */
 412     if (n2 < BN_MUL_RECURSIVE_SIZE_NORMAL) {
 413         bn_mul_normal(r, a, n2 + dna, b, n2 + dnb);
 414         if ((dna + dnb) < 0)
 415             memset(&r[2 * n2 + dna + dnb], 0,
 416                    sizeof(BN_ULONG) * -(dna + dnb));
 417         return;
 418     }
 419     /* r=(a[0]-a[1])*(b[1]-b[0]) */
 420     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
 421     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
 422     zero = neg = 0;
 423     switch (c1 * 3 + c2) {
 424     case -4:
 425         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
 426         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
 427         break;
 428     case -3:
 429         zero = 1;
 430         break;
 431     case -2:
 432         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
 433         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
 434         neg = 1;
 435         break;
 436     case -1:
 437     case 0:
 438     case 1:
 439         zero = 1;
 440         break;
 441     case 2:
 442         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
 443         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
 444         neg = 1;
 445         break;
 446     case 3:
 447         zero = 1;
 448         break;
 449     case 4:
 450         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
 451         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
 452         break;
 453     }
 454
 455 # ifdef BN_MUL_COMBA
 456     if (n == 4 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba4 could take
 457                                            * extra args to do this well */
 458         if (!zero)
 459             bn_mul_comba4(&(t[n2]), t, &(t[n]));
 460         else
 461             memset(&t[n2], 0, sizeof(*t) * 8);
 462
 463         bn_mul_comba4(r, a, b);
 464         bn_mul_comba4(&(r[n2]), &(a[n]), &(b[n]));
 465     } else if (n == 8 && dna == 0 && dnb == 0) { /* XXX: bn_mul_comba8 could
 466                                                   * take extra args to do
 467                                                   * this well */
 468         if (!zero)
 469             bn_mul_comba8(&(t[n2]), t, &(t[n]));
 470         else
 471             memset(&t[n2], 0, sizeof(*t) * 16);
 472
 473         bn_mul_comba8(r, a, b);
 474         bn_mul_comba8(&(r[n2]), &(a[n]), &(b[n]));
 475     } else
 476 # endif                         /* BN_MUL_COMBA */
 477     {
 478         p = &(t[n2 * 2]);
 479         if (!zero)
 480             bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
 481         else
 482             memset(&t[n2], 0, sizeof(*t) * n2);
 483         bn_mul_recursive(r, a, b, n, 0, 0, p);
 484         bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]), n, dna, dnb, p);
 485     }
 486
 487     /*-
 488      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
 489      * r[10] holds (a[0]*b[0])
 490      * r[32] holds (b[1]*b[1])
 491      */
 492
 493     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
 494
 495     if (neg) {                  /* if t[32] is negative */
 496         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
 497     } else {
 498         /* Might have a carry */
 499         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
 500     }
 501
 502     /*-
 503      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
 504      * r[10] holds (a[0]*b[0])
 505      * r[32] holds (b[1]*b[1])
 506      * c1 holds the carry bits
 507      */
 508     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
 509     if (c1) {
 510         p = &(r[n + n2]);
 511         lo = *p;
 512         ln = (lo + c1) & BN_MASK2;
 513         *p = ln;
 514
 515         /*
 516          * The overflow will stop before we over write words we should not
 517          * overwrite
 518          */
 519         if (ln < (BN_ULONG)c1) {
 520             do {
 521                 p++;
 522                 lo = *p;
 523                 ln = (lo + 1) & BN_MASK2;
 524                 *p = ln;
 525             } while (ln == 0);
 526         }
 527     }
 528 }
 529
 530 /*
 531  * n+tn is the word length t needs to be n*4 is size, as does r
 532  */
 533 /* tnX may not be negative but less than n */
 534 void bn_mul_part_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n,
 535                            int tna, int tnb, BN_ULONG *t)
 536 {
 537     int i, j, n2 = n * 2;
 538     int c1, c2, neg;
 539     BN_ULONG ln, lo, *p;
 540
 541     if (n < 8) {
 542         bn_mul_normal(r, a, n + tna, b, n + tnb);
 543         return;
 544     }
 545
 546     /* r=(a[0]-a[1])*(b[1]-b[0]) */
 547     c1 = bn_cmp_part_words(a, &(a[n]), tna, n - tna);
 548     c2 = bn_cmp_part_words(&(b[n]), b, tnb, tnb - n);
 549     neg = 0;
 550     switch (c1 * 3 + c2) {
 551     case -4:
 552         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
 553         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
 554         break;
 555     case -3:
 556         /* break; */
 557     case -2:
 558         bn_sub_part_words(t, &(a[n]), a, tna, tna - n); /* - */
 559         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n); /* + */
 560         neg = 1;
 561         break;
 562     case -1:
 563     case 0:
 564     case 1:
 565         /* break; */
 566     case 2:
 567         bn_sub_part_words(t, a, &(a[n]), tna, n - tna); /* + */
 568         bn_sub_part_words(&(t[n]), b, &(b[n]), tnb, n - tnb); /* - */
 569         neg = 1;
 570         break;
 571     case 3:
 572         /* break; */
 573     case 4:
 574         bn_sub_part_words(t, a, &(a[n]), tna, n - tna);
 575         bn_sub_part_words(&(t[n]), &(b[n]), b, tnb, tnb - n);
 576         break;
 577     }
 578     /*
 579      * The zero case isn't yet implemented here. The speedup would probably
 580      * be negligible.
 581      */
 582 # if 0
 583     if (n == 4) {
 584         bn_mul_comba4(&(t[n2]), t, &(t[n]));
 585         bn_mul_comba4(r, a, b);
 586         bn_mul_normal(&(r[n2]), &(a[n]), tn, &(b[n]), tn);
 587         memset(&r[n2 + tn * 2], 0, sizeof(*r) * (n2 - tn * 2));
 588     } else
 589 # endif
 590     if (n == 8) {
 591         bn_mul_comba8(&(t[n2]), t, &(t[n]));
 592         bn_mul_comba8(r, a, b);
 593         bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
 594         memset(&r[n2 + tna + tnb], 0, sizeof(*r) * (n2 - tna - tnb));
 595     } else {
 596         p = &(t[n2 * 2]);
 597         bn_mul_recursive(&(t[n2]), t, &(t[n]), n, 0, 0, p);
 598         bn_mul_recursive(r, a, b, n, 0, 0, p);
 599         i = n / 2;
 600         /*
 601          * If there is only a bottom half to the number, just do it
 602          */
 603         if (tna > tnb)
 604             j = tna - i;
 605         else
 606             j = tnb - i;
 607         if (j == 0) {
 608             bn_mul_recursive(&(r[n2]), &(a[n]), &(b[n]),
 609                              i, tna - i, tnb - i, p);
 610             memset(&r[n2 + i * 2], 0, sizeof(*r) * (n2 - i * 2));
 611         } else if (j > 0) {     /* eg, n == 16, i == 8 and tn == 11 */
 612             bn_mul_part_recursive(&(r[n2]), &(a[n]), &(b[n]),
 613                                   i, tna - i, tnb - i, p);
 614             memset(&(r[n2 + tna + tnb]), 0,
 615                    sizeof(BN_ULONG) * (n2 - tna - tnb));
 616         } else {                /* (j < 0) eg, n == 16, i == 8 and tn == 5 */
 617
 618             memset(&r[n2], 0, sizeof(*r) * n2);
 619             if (tna < BN_MUL_RECURSIVE_SIZE_NORMAL
 620                 && tnb < BN_MUL_RECURSIVE_SIZE_NORMAL) {
 621                 bn_mul_normal(&(r[n2]), &(a[n]), tna, &(b[n]), tnb);
 622             } else {
 623                 for (;;) {
 624                     i /= 2;
 625                     /*
 626                      * these simplified conditions work exclusively because
 627                      * difference between tna and tnb is 1 or 0
 628                      */
 629                     if (i < tna || i < tnb) {
 630                         bn_mul_part_recursive(&(r[n2]),
 631                                               &(a[n]), &(b[n]),
 632                                               i, tna - i, tnb - i, p);
 633                         break;
 634                     } else if (i == tna || i == tnb) {
 635                         bn_mul_recursive(&(r[n2]),
 636                                          &(a[n]), &(b[n]),
 637                                          i, tna - i, tnb - i, p);
 638                         break;
 639                     }
 640                 }
 641             }
 642         }
 643     }
 644
 645     /*-
 646      * t[32] holds (a[0]-a[1])*(b[1]-b[0]), c1 is the sign
 647      * r[10] holds (a[0]*b[0])
 648      * r[32] holds (b[1]*b[1])
 649      */
 650
 651     c1 = (int)(bn_add_words(t, r, &(r[n2]), n2));
 652
 653     if (neg) {                  /* if t[32] is negative */
 654         c1 -= (int)(bn_sub_words(&(t[n2]), t, &(t[n2]), n2));
 655     } else {
 656         /* Might have a carry */
 657         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), t, n2));
 658     }
 659
 660     /*-
 661      * t[32] holds (a[0]-a[1])*(b[1]-b[0])+(a[0]*b[0])+(a[1]*b[1])
 662      * r[10] holds (a[0]*b[0])
 663      * r[32] holds (b[1]*b[1])
 664      * c1 holds the carry bits
 665      */
 666     c1 += (int)(bn_add_words(&(r[n]), &(r[n]), &(t[n2]), n2));
 667     if (c1) {
 668         p = &(r[n + n2]);
 669         lo = *p;
 670         ln = (lo + c1) & BN_MASK2;
 671         *p = ln;
 672
 673         /*
 674          * The overflow will stop before we over write words we should not
 675          * overwrite
 676          */
 677         if (ln < (BN_ULONG)c1) {
 678             do {
 679                 p++;
 680                 lo = *p;
 681                 ln = (lo + 1) & BN_MASK2;
 682                 *p = ln;
 683             } while (ln == 0);
 684         }
 685     }
 686 }
 687
 688 /*-
 689  * a and b must be the same size, which is n2.
 690  * r needs to be n2 words and t needs to be n2*2
 691  */
 692 void bn_mul_low_recursive(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n2,
 693                           BN_ULONG *t)
 694 {
 695     int n = n2 / 2;
 696
 697     bn_mul_recursive(r, a, b, n, 0, 0, &(t[0]));
 698     if (n >= BN_MUL_LOW_RECURSIVE_SIZE_NORMAL) {
 699         bn_mul_low_recursive(&(t[0]), &(a[0]), &(b[n]), n, &(t[n2]));
 700         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
 701         bn_mul_low_recursive(&(t[0]), &(a[n]), &(b[0]), n, &(t[n2]));
 702         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
 703     } else {
 704         bn_mul_low_normal(&(t[0]), &(a[0]), &(b[n]), n);
 705         bn_mul_low_normal(&(t[n]), &(a[n]), &(b[0]), n);
 706         bn_add_words(&(r[n]), &(r[n]), &(t[0]), n);
 707         bn_add_words(&(r[n]), &(r[n]), &(t[n]), n);
 708     }
 709 }
 710
 711 /*-
 712  * a and b must be the same size, which is n2.
 713  * r needs to be n2 words and t needs to be n2*2
 714  * l is the low words of the output.
 715  * t needs to be n2*3
 716  */
 717 void bn_mul_high(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, BN_ULONG *l, int n2,
 718                  BN_ULONG *t)
 719 {
 720     int i, n;
 721     int c1, c2;
 722     int neg, oneg, zero;
 723     BN_ULONG ll, lc, *lp, *mp;
 724
 725     n = n2 / 2;
 726
 727     /* Calculate (al-ah)*(bh-bl) */
 728     neg = zero = 0;
 729     c1 = bn_cmp_words(&(a[0]), &(a[n]), n);
 730     c2 = bn_cmp_words(&(b[n]), &(b[0]), n);
 731     switch (c1 * 3 + c2) {
 732     case -4:
 733         bn_sub_words(&(r[0]), &(a[n]), &(a[0]), n);
 734         bn_sub_words(&(r[n]), &(b[0]), &(b[n]), n);
 735         break;
 736     case -3:
 737         zero = 1;
 738         break;
 739     case -2:
 740         bn_sub_words(&(r[0]), &(a[n]), &(a[0]), n);
 741         bn_sub_words(&(r[n]), &(b[n]), &(b[0]), n);
 742         neg = 1;
 743         break;
 744     case -1:
 745     case 0:
 746     case 1:
 747         zero = 1;
 748         break;
 749     case 2:
 750         bn_sub_words(&(r[0]), &(a[0]), &(a[n]), n);
 751         bn_sub_words(&(r[n]), &(b[0]), &(b[n]), n);
 752         neg = 1;
 753         break;
 754     case 3:
 755         zero = 1;
 756         break;
 757     case 4:
 758         bn_sub_words(&(r[0]), &(a[0]), &(a[n]), n);
 759         bn_sub_words(&(r[n]), &(b[n]), &(b[0]), n);
 760         break;
 761     }
 762
 763     oneg = neg;
 764     /* t[10] = (a[0]-a[1])*(b[1]-b[0]) */
 765     /* r[10] = (a[1]*b[1]) */
 766 # ifdef BN_MUL_COMBA
 767     if (n == 8) {
 768         bn_mul_comba8(&(t[0]), &(r[0]), &(r[n]));
 769         bn_mul_comba8(r, &(a[n]), &(b[n]));
 770     } else
 771 # endif
 772     {
 773         bn_mul_recursive(&(t[0]), &(r[0]), &(r[n]), n, 0, 0, &(t[n2]));
 774         bn_mul_recursive(r, &(a[n]), &(b[n]), n, 0, 0, &(t[n2]));
 775     }
 776
 777     /*-
 778      * s0 == low(al*bl)
 779      * s1 == low(ah*bh)+low((al-ah)*(bh-bl))+low(al*bl)+high(al*bl)
 780      * We know s0 and s1 so the only unknown is high(al*bl)
 781      * high(al*bl) == s1 - low(ah*bh+s0+(al-ah)*(bh-bl))
 782      * high(al*bl) == s1 - (r[0]+l[0]+t[0])
 783      */
 784     if (l != NULL) {
 785         lp = &(t[n2 + n]);
 786         c1 = (int)(bn_add_words(lp, &(r[0]), &(l[0]), n));
 787     } else {
 788         c1 = 0;
 789         lp = &(r[0]);
 790     }
 791
 792     if (neg)
 793         neg = (int)(bn_sub_words(&(t[n2]), lp, &(t[0]), n));
 794     else {
 795         bn_add_words(&(t[n2]), lp, &(t[0]), n);
 796         neg = 0;
 797     }
 798
 799     if (l != NULL) {
 800         bn_sub_words(&(t[n2 + n]), &(l[n]), &(t[n2]), n);
 801     } else {
 802         lp = &(t[n2 + n]);
 803         mp = &(t[n2]);
 804         for (i = 0; i < n; i++)
 805             lp[i] = ((~mp[i]) + 1) & BN_MASK2;
 806     }
 807
 808     /*-
 809      * s[0] = low(al*bl)
 810      * t[3] = high(al*bl)
 811      * t[10] = (a[0]-a[1])*(b[1]-b[0]) neg is the sign
 812      * r[10] = (a[1]*b[1])
 813      */
 814     /*-
 815      * R[10] = al*bl
 816      * R[21] = al*bl + ah*bh + (a[0]-a[1])*(b[1]-b[0])
 817      * R[32] = ah*bh
 818      */
 819     /*-
 820      * R[1]=t[3]+l[0]+r[0](+-)t[0] (have carry/borrow)
 821      * R[2]=r[0]+t[3]+r[1](+-)t[1] (have carry/borrow)
 822      * R[3]=r[1]+(carry/borrow)
 823      */
 824     if (l != NULL) {
 825         lp = &(t[n2]);
 826         c1 = (int)(bn_add_words(lp, &(t[n2 + n]), &(l[0]), n));
 827     } else {
 828         lp = &(t[n2 + n]);
 829         c1 = 0;
 830     }
 831     c1 += (int)(bn_add_words(&(t[n2]), lp, &(r[0]), n));
 832     if (oneg)
 833         c1 -= (int)(bn_sub_words(&(t[n2]), &(t[n2]), &(t[0]), n));
 834     else
 835         c1 += (int)(bn_add_words(&(t[n2]), &(t[n2]), &(t[0]), n));
 836
 837     c2 = (int)(bn_add_words(&(r[0]), &(r[0]), &(t[n2 + n]), n));
 838     c2 += (int)(bn_add_words(&(r[0]), &(r[0]), &(r[n]), n));
 839     if (oneg)
 840         c2 -= (int)(bn_sub_words(&(r[0]), &(r[0]), &(t[n]), n));
 841     else
 842         c2 += (int)(bn_add_words(&(r[0]), &(r[0]), &(t[n]), n));
 843
 844     if (c1 != 0) {              /* Add starting at r[0], could be +ve or -ve */
 845         i = 0;
 846         if (c1 > 0) {
 847             lc = c1;
 848             do {
 849                 ll = (r[i] + lc) & BN_MASK2;
 850                 r[i++] = ll;
 851                 lc = (lc > ll);
 852             } while (lc);
 853         } else {
 854             lc = -c1;
 855             do {
 856                 ll = r[i];
 857                 r[i++] = (ll - lc) & BN_MASK2;
 858                 lc = (lc > ll);
 859             } while (lc);
 860         }
 861     }
 862     if (c2 != 0) {              /* Add starting at r[1] */
 863         i = n;
 864         if (c2 > 0) {
 865             lc = c2;
 866             do {
 867                 ll = (r[i] + lc) & BN_MASK2;
 868                 r[i++] = ll;
 869                 lc = (lc > ll);
 870             } while (lc);
 871         } else {
 872             lc = -c2;
 873             do {
 874                 ll = r[i];
 875                 r[i++] = (ll - lc) & BN_MASK2;
 876                 lc = (lc > ll);
 877             } while (lc);
 878         }
 879     }
 880 }
 881 #endif                          /* BN_RECURSION */
 882
 883 int BN_mul(BIGNUM *r, const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
 884 {
 885     int ret = 0;
 886     int top, al, bl;
 887     BIGNUM *rr;
 888 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
 889     int i;
 890 #endif
 891 #ifdef BN_RECURSION
 892     BIGNUM *t = NULL;
 893     int j = 0, k;
 894 #endif
 895
 896     bn_check_top(a);
 897     bn_check_top(b);
 898     bn_check_top(r);
 899
 900     al = a->top;
 901     bl = b->top;
 902
 903     if ((al == 0) || (bl == 0)) {
 904         BN_zero(r);
 905         return (1);
 906     }
 907     top = al + bl;
 908
 909     BN_CTX_start(ctx);
 910     if ((r == a) || (r == b)) {
 911         if ((rr = BN_CTX_get(ctx)) == NULL)
 912             goto err;
 913     } else
 914         rr = r;
 915     rr->neg = a->neg ^ b->neg;
 916
 917 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
 918     i = al - bl;
 919 #endif
 920 #ifdef BN_MUL_COMBA
 921     if (i == 0) {
 922 # if 0
 923         if (al == 4) {
 924             if (bn_wexpand(rr, 8) == NULL)
 925                 goto err;
 926             rr->top = 8;
 927             bn_mul_comba4(rr->d, a->d, b->d);
 928             goto end;
 929         }
 930 # endif
 931         if (al == 8) {
 932             if (bn_wexpand(rr, 16) == NULL)
 933                 goto err;
 934             rr->top = 16;
 935             bn_mul_comba8(rr->d, a->d, b->d);
 936             goto end;
 937         }
 938     }
 939 #endif                          /* BN_MUL_COMBA */
 940 #ifdef BN_RECURSION
 941     if ((al >= BN_MULL_SIZE_NORMAL) && (bl >= BN_MULL_SIZE_NORMAL)) {
 942         if (i >= -1 && i <= 1) {
 943             /*
 944              * Find out the power of two lower or equal to the longest of the
 945              * two numbers
 946              */
 947             if (i >= 0) {
 948                 j = BN_num_bits_word((BN_ULONG)al);
 949             }
 950             if (i == -1) {
 951                 j = BN_num_bits_word((BN_ULONG)bl);
 952             }
 953             j = 1 << (j - 1);
 954             assert(j <= al || j <= bl);
 955             k = j + j;
 956             t = BN_CTX_get(ctx);
 957             if (t == NULL)
 958                 goto err;
 959             if (al > j || bl > j) {
 960                 if (bn_wexpand(t, k * 4) == NULL)
 961                     goto err;
 962                 if (bn_wexpand(rr, k * 4) == NULL)
 963                     goto err;
 964                 bn_mul_part_recursive(rr->d, a->d, b->d,
 965                                       j, al - j, bl - j, t->d);
 966             } else {            /* al <= j || bl <= j */
 967
 968                 if (bn_wexpand(t, k * 2) == NULL)
 969                     goto err;
 970                 if (bn_wexpand(rr, k * 2) == NULL)
 971                     goto err;
 972                 bn_mul_recursive(rr->d, a->d, b->d, j, al - j, bl - j, t->d);
 973             }
 974             rr->top = top;
 975             goto end;
 976         }
 977 # if 0
 978         if (i == 1 && !BN_get_flags(b, BN_FLG_STATIC_DATA)) {
 979             BIGNUM *tmp_bn = (BIGNUM *)b;
 980             if (bn_wexpand(tmp_bn, al) == NULL)
 981                 goto err;
 982             tmp_bn->d[bl] = 0;
 983             bl++;
 984             i--;
 985         } else if (i == -1 && !BN_get_flags(a, BN_FLG_STATIC_DATA)) {
 986             BIGNUM *tmp_bn = (BIGNUM *)a;
 987             if (bn_wexpand(tmp_bn, bl) == NULL)
 988                 goto err;
 989             tmp_bn->d[al] = 0;
 990             al++;
 991             i++;
 992         }
 993         if (i == 0) {
 994             /* symmetric and > 4 */
 995             /* 16 or larger */
 996             j = BN_num_bits_word((BN_ULONG)al);
 997             j = 1 << (j - 1);
 998             k = j + j;
 999             t = BN_CTX_get(ctx);
1000             if (al == j) {      /* exact multiple */
1001                 if (bn_wexpand(t, k * 2) == NULL)
1002                     goto err;
1003                 if (bn_wexpand(rr, k * 2) == NULL)
1004                     goto err;
1005                 bn_mul_recursive(rr->d, a->d, b->d, al, t->d);
1006             } else {
1007                 if (bn_wexpand(t, k * 4) == NULL)
1008                     goto err;
1009                 if (bn_wexpand(rr, k * 4) == NULL)
1010                     goto err;
1011                 bn_mul_part_recursive(rr->d, a->d, b->d, al - j, j, t->d);
1012             }
1013             rr->top = top;
1014             goto end;
1015         }
1016 # endif
1017     }
1018 #endif                          /* BN_RECURSION */
1019     if (bn_wexpand(rr, top) == NULL)
1020         goto err;
1021     rr->top = top;
1022     bn_mul_normal(rr->d, a->d, al, b->d, bl);
1023
1024 #if defined(BN_MUL_COMBA) || defined(BN_RECURSION)
1025  end:
1026 #endif
1027     bn_correct_top(rr);
1028     if (r != rr)
1029         BN_copy(r, rr);
1030     ret = 1;
1031  err:
1032     bn_check_top(r);
1033     BN_CTX_end(ctx);
1034     return (ret);
1035 }
1036
1037 void bn_mul_normal(BN_ULONG *r, BN_ULONG *a, int na, BN_ULONG *b, int nb)
1038 {
1039     BN_ULONG *rr;
1040
1041     if (na < nb) {
1042         int itmp;
1043         BN_ULONG *ltmp;
1044
1045         itmp = na;
1046         na = nb;
1047         nb = itmp;
1048         ltmp = a;
1049         a = b;
1050         b = ltmp;
1051
1052     }
1053     rr = &(r[na]);
1054     if (nb <= 0) {
1055         (void)bn_mul_words(r, a, na, 0);
1056         return;
1057     } else
1058         rr[0] = bn_mul_words(r, a, na, b[0]);
1059
1060     for (;;) {
1061         if (--nb <= 0)
1062             return;
1063         rr[1] = bn_mul_add_words(&(r[1]), a, na, b[1]);
1064         if (--nb <= 0)
1065             return;
1066         rr[2] = bn_mul_add_words(&(r[2]), a, na, b[2]);
1067         if (--nb <= 0)
1068             return;
1069         rr[3] = bn_mul_add_words(&(r[3]), a, na, b[3]);
1070         if (--nb <= 0)
1071             return;
1072         rr[4] = bn_mul_add_words(&(r[4]), a, na, b[4]);
1073         rr += 4;
1074         r += 4;
1075         b += 4;
1076     }
1077 }
1078
1079 void bn_mul_low_normal(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n)
1080 {
1081     bn_mul_words(r, a, n, b[0]);
1082
1083     for (;;) {
1084         if (--n <= 0)
1085             return;
1086         bn_mul_add_words(&(r[1]), a, n, b[1]);
1087         if (--n <= 0)
1088             return;
1089         bn_mul_add_words(&(r[2]), a, n, b[2]);
1090         if (--n <= 0)
1091             return;
1092         bn_mul_add_words(&(r[3]), a, n, b[3]);
1093         if (--n <= 0)
1094             return;
1095         bn_mul_add_words(&(r[4]), a, n, b[4]);
1096         r += 4;
1097         b += 4;
1098     }
1099 }