crypto/bn/bn_mul.c

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