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