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