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