1 /*
   2  * Copyright (c) 1999, 2013, Oracle and/or its affiliates. All rights reserved.
   3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   4  *
   5  * This code is free software; you can redistribute it and/or modify it
   6  * under the terms of the GNU General Public License version 2 only, as
   7  * published by the Free Software Foundation.  Oracle designates this
   8  * particular file as subject to the "Classpath" exception as provided
   9  * by Oracle in the LICENSE file that accompanied this code.
  10  *
  11  * This code is distributed in the hope that it will be useful, but WITHOUT
  12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14  * version 2 for more details (a copy is included in the LICENSE file that
  15  * accompanied this code).
  16  *
  17  * You should have received a copy of the GNU General Public License version
  18  * 2 along with this work; if not, write to the Free Software Foundation,
  19  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  20  *
  21  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
  22  * or visit www.oracle.com if you need additional information or have any
  23  * questions.
  24  */
  25 
  26 package java.math;
  27 
  28 /**
  29  * A class used to represent multiprecision integers that makes efficient
  30  * use of allocated space by allowing a number to occupy only part of
  31  * an array so that the arrays do not have to be reallocated as often.
  32  * When performing an operation with many iterations the array used to
  33  * hold a number is only reallocated when necessary and does not have to
  34  * be the same size as the number it represents. A mutable number allows
  35  * calculations to occur on the same number without having to create
  36  * a new number for every step of the calculation as occurs with
  37  * BigIntegers.
  38  *
  39  * @see     BigInteger
  40  * @author  Michael McCloskey
  41  * @author  Timothy Buktu
  42  * @since   1.3
  43  */
  44 
  45 import static java.math.BigDecimal.INFLATED;
  46 import static java.math.BigInteger.LONG_MASK;
  47 import java.util.Arrays;
  48 
  49 class MutableBigInteger {
  50     /**
  51      * Holds the magnitude of this MutableBigInteger in big endian order.
  52      * The magnitude may start at an offset into the value array, and it may
  53      * end before the length of the value array.
  54      */
  55     int[] value;
  56 
  57     /**
  58      * The number of ints of the value array that are currently used
  59      * to hold the magnitude of this MutableBigInteger. The magnitude starts
  60      * at an offset and offset + intLen may be less than value.length.
  61      */
  62     int intLen;
  63 
  64     /**
  65      * The offset into the value array where the magnitude of this
  66      * MutableBigInteger begins.
  67      */
  68     int offset = 0;
  69 
  70     // Constants
  71     /**
  72      * MutableBigInteger with one element value array with the value 1. Used by
  73      * BigDecimal divideAndRound to increment the quotient. Use this constant
  74      * only when the method is not going to modify this object.
  75      */
  76     static final MutableBigInteger ONE = new MutableBigInteger(1);
  77 
  78     /**
  79      * The minimum {@code intLen} for cancelling powers of two before
  80      * dividing.
  81      * If the number of ints is less than this threshold,
  82      * {@code divideKnuth} does not eliminate common powers of two from
  83      * the dividend and divisor.
  84      */
  85     static final int KNUTH_POW2_THRESH_LEN = 6;
  86 
  87     /**
  88      * The minimum number of trailing zero ints for cancelling powers of two
  89      * before dividing.
  90      * If the dividend and divisor don't share at least this many zero ints
  91      * at the end, {@code divideKnuth} does not eliminate common powers
  92      * of two from the dividend and divisor.
  93      */
  94     static final int KNUTH_POW2_THRESH_ZEROS = 3;
  95 
  96     // Constructors
  97 
  98     /**
  99      * The default constructor. An empty MutableBigInteger is created with
 100      * a one word capacity.
 101      */
 102     MutableBigInteger() {
 103         value = new int[1];
 104         intLen = 0;
 105     }
 106 
 107     /**
 108      * Construct a new MutableBigInteger with a magnitude specified by
 109      * the int val.
 110      */
 111     MutableBigInteger(int val) {
 112         value = new int[1];
 113         intLen = 1;
 114         value[0] = val;
 115     }
 116 
 117     /**
 118      * Construct a new MutableBigInteger with the specified value array
 119      * up to the length of the array supplied.
 120      */
 121     MutableBigInteger(int[] val) {
 122         value = val;
 123         intLen = val.length;
 124     }
 125 
 126     /**
 127      * Construct a new MutableBigInteger with a magnitude equal to the
 128      * specified BigInteger.
 129      */
 130     MutableBigInteger(BigInteger b) {
 131         intLen = b.mag.length;
 132         value = Arrays.copyOf(b.mag, intLen);
 133     }
 134 
 135     /**
 136      * Construct a new MutableBigInteger with a magnitude equal to the
 137      * specified MutableBigInteger.
 138      */
 139     MutableBigInteger(MutableBigInteger val) {
 140         intLen = val.intLen;
 141         value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
 142     }
 143 
 144     /**
 145      * Makes this number an {@code n}-int number all of whose bits are ones.
 146      * Used by Burnikel-Ziegler division.
 147      * @param n number of ints in the {@code value} array
 148      * @return a number equal to {@code ((1<<(32*n)))-1}
 149      */
 150     private void ones(int n) {
 151         if (n > value.length)
 152             value = new int[n];
 153         Arrays.fill(value, -1);
 154         offset = 0;
 155         intLen = n;
 156     }
 157 
 158     /**
 159      * Internal helper method to return the magnitude array. The caller is not
 160      * supposed to modify the returned array.
 161      */
 162     private int[] getMagnitudeArray() {
 163         if (offset > 0 || value.length != intLen)
 164             return Arrays.copyOfRange(value, offset, offset + intLen);
 165         return value;
 166     }
 167 
 168     /**
 169      * Convert this MutableBigInteger to a long value. The caller has to make
 170      * sure this MutableBigInteger can be fit into long.
 171      */
 172     private long toLong() {
 173         assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
 174         if (intLen == 0)
 175             return 0;
 176         long d = value[offset] & LONG_MASK;
 177         return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
 178     }
 179 
 180     /**
 181      * Convert this MutableBigInteger to a BigInteger object.
 182      */
 183     BigInteger toBigInteger(int sign) {
 184         if (intLen == 0 || sign == 0)
 185             return BigInteger.ZERO;
 186         return new BigInteger(getMagnitudeArray(), sign);
 187     }
 188 
 189     /**
 190      * Converts this number to a nonnegative {@code BigInteger}.
 191      */
 192     BigInteger toBigInteger() {
 193         normalize();
 194         return toBigInteger(isZero() ? 0 : 1);
 195     }
 196 
 197     /**
 198      * Convert this MutableBigInteger to BigDecimal object with the specified sign
 199      * and scale.
 200      */
 201     BigDecimal toBigDecimal(int sign, int scale) {
 202         if (intLen == 0 || sign == 0)
 203             return BigDecimal.zeroValueOf(scale);
 204         int[] mag = getMagnitudeArray();
 205         int len = mag.length;
 206         int d = mag[0];
 207         // If this MutableBigInteger can't be fit into long, we need to
 208         // make a BigInteger object for the resultant BigDecimal object.
 209         if (len > 2 || (d < 0 && len == 2))
 210             return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
 211         long v = (len == 2) ?
 212             ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
 213             d & LONG_MASK;
 214         return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
 215     }
 216 
 217     /**
 218      * This is for internal use in converting from a MutableBigInteger
 219      * object into a long value given a specified sign.
 220      * returns INFLATED if value is not fit into long
 221      */
 222     long toCompactValue(int sign) {
 223         if (intLen == 0 || sign == 0)
 224             return 0L;
 225         int[] mag = getMagnitudeArray();
 226         int len = mag.length;
 227         int d = mag[0];
 228         // If this MutableBigInteger can not be fitted into long, we need to
 229         // make a BigInteger object for the resultant BigDecimal object.
 230         if (len > 2 || (d < 0 && len == 2))
 231             return INFLATED;
 232         long v = (len == 2) ?
 233             ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
 234             d & LONG_MASK;
 235         return sign == -1 ? -v : v;
 236     }
 237 
 238     /**
 239      * Clear out a MutableBigInteger for reuse.
 240      */
 241     void clear() {
 242         offset = intLen = 0;
 243         for (int index=0, n=value.length; index < n; index++)
 244             value[index] = 0;
 245     }
 246 
 247     /**
 248      * Set a MutableBigInteger to zero, removing its offset.
 249      */
 250     void reset() {
 251         offset = intLen = 0;
 252     }
 253 
 254     /**
 255      * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
 256      * as this MutableBigInteger is numerically less than, equal to, or
 257      * greater than <tt>b</tt>.
 258      */
 259     final int compare(MutableBigInteger b) {
 260         int blen = b.intLen;
 261         if (intLen < blen)
 262             return -1;
 263         if (intLen > blen)
 264            return 1;
 265 
 266         // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
 267         // comparison.
 268         int[] bval = b.value;
 269         for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
 270             int b1 = value[i] + 0x80000000;
 271             int b2 = bval[j]  + 0x80000000;
 272             if (b1 < b2)
 273                 return -1;
 274             if (b1 > b2)
 275                 return 1;
 276         }
 277         return 0;
 278     }
 279 
 280     /**
 281      * Returns a value equal to what {@code b.leftShift(32*ints); return compare(b);}
 282      * would return, but doesn't change the value of {@code b}.
 283      */
 284     private int compareShifted(MutableBigInteger b, int ints) {
 285         int blen = b.intLen;
 286         int alen = intLen - ints;
 287         if (alen < blen)
 288             return -1;
 289         if (alen > blen)
 290            return 1;
 291 
 292         // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
 293         // comparison.
 294         int[] bval = b.value;
 295         for (int i = offset, j = b.offset; i < alen + offset; i++, j++) {
 296             int b1 = value[i] + 0x80000000;
 297             int b2 = bval[j]  + 0x80000000;
 298             if (b1 < b2)
 299                 return -1;
 300             if (b1 > b2)
 301                 return 1;
 302         }
 303         return 0;
 304     }
 305 
 306     /**
 307      * Compare this against half of a MutableBigInteger object (Needed for
 308      * remainder tests).
 309      * Assumes no leading unnecessary zeros, which holds for results
 310      * from divide().
 311      */
 312     final int compareHalf(MutableBigInteger b) {
 313         int blen = b.intLen;
 314         int len = intLen;
 315         if (len <= 0)
 316             return blen <= 0 ? 0 : -1;
 317         if (len > blen)
 318             return 1;
 319         if (len < blen - 1)
 320             return -1;
 321         int[] bval = b.value;
 322         int bstart = 0;
 323         int carry = 0;
 324         // Only 2 cases left:len == blen or len == blen - 1
 325         if (len != blen) { // len == blen - 1
 326             if (bval[bstart] == 1) {
 327                 ++bstart;
 328                 carry = 0x80000000;
 329             } else
 330                 return -1;
 331         }
 332         // compare values with right-shifted values of b,
 333         // carrying shifted-out bits across words
 334         int[] val = value;
 335         for (int i = offset, j = bstart; i < len + offset;) {
 336             int bv = bval[j++];
 337             long hb = ((bv >>> 1) + carry) & LONG_MASK;
 338             long v = val[i++] & LONG_MASK;
 339             if (v != hb)
 340                 return v < hb ? -1 : 1;
 341             carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
 342         }
 343         return carry == 0 ? 0 : -1;
 344     }
 345 
 346     /**
 347      * Return the index of the lowest set bit in this MutableBigInteger. If the
 348      * magnitude of this MutableBigInteger is zero, -1 is returned.
 349      */
 350     private final int getLowestSetBit() {
 351         if (intLen == 0)
 352             return -1;
 353         int j, b;
 354         for (j=intLen-1; (j > 0) && (value[j+offset] == 0); j--)
 355             ;
 356         b = value[j+offset];
 357         if (b == 0)
 358             return -1;
 359         return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
 360     }
 361 
 362     /**
 363      * Return the int in use in this MutableBigInteger at the specified
 364      * index. This method is not used because it is not inlined on all
 365      * platforms.
 366      */
 367     private final int getInt(int index) {
 368         return value[offset+index];
 369     }
 370 
 371     /**
 372      * Return a long which is equal to the unsigned value of the int in
 373      * use in this MutableBigInteger at the specified index. This method is
 374      * not used because it is not inlined on all platforms.
 375      */
 376     private final long getLong(int index) {
 377         return value[offset+index] & LONG_MASK;
 378     }
 379 
 380     /**
 381      * Ensure that the MutableBigInteger is in normal form, specifically
 382      * making sure that there are no leading zeros, and that if the
 383      * magnitude is zero, then intLen is zero.
 384      */
 385     final void normalize() {
 386         if (intLen == 0) {
 387             offset = 0;
 388             return;
 389         }
 390 
 391         int index = offset;
 392         if (value[index] != 0)
 393             return;
 394 
 395         int indexBound = index+intLen;
 396         do {
 397             index++;
 398         } while(index < indexBound && value[index] == 0);
 399 
 400         int numZeros = index - offset;
 401         intLen -= numZeros;
 402         offset = (intLen == 0 ?  0 : offset+numZeros);
 403     }
 404 
 405     /**
 406      * If this MutableBigInteger cannot hold len words, increase the size
 407      * of the value array to len words.
 408      */
 409     private final void ensureCapacity(int len) {
 410         if (value.length < len) {
 411             value = new int[len];
 412             offset = 0;
 413             intLen = len;
 414         }
 415     }
 416 
 417     /**
 418      * Convert this MutableBigInteger into an int array with no leading
 419      * zeros, of a length that is equal to this MutableBigInteger's intLen.
 420      */
 421     int[] toIntArray() {
 422         int[] result = new int[intLen];
 423         for(int i=0; i < intLen; i++)
 424             result[i] = value[offset+i];
 425         return result;
 426     }
 427 
 428     /**
 429      * Sets the int at index+offset in this MutableBigInteger to val.
 430      * This does not get inlined on all platforms so it is not used
 431      * as often as originally intended.
 432      */
 433     void setInt(int index, int val) {
 434         value[offset + index] = val;
 435     }
 436 
 437     /**
 438      * Sets this MutableBigInteger's value array to the specified array.
 439      * The intLen is set to the specified length.
 440      */
 441     void setValue(int[] val, int length) {
 442         value = val;
 443         intLen = length;
 444         offset = 0;
 445     }
 446 
 447     /**
 448      * Sets this MutableBigInteger's value array to a copy of the specified
 449      * array. The intLen is set to the length of the new array.
 450      */
 451     void copyValue(MutableBigInteger src) {
 452         int len = src.intLen;
 453         if (value.length < len)
 454             value = new int[len];
 455         System.arraycopy(src.value, src.offset, value, 0, len);
 456         intLen = len;
 457         offset = 0;
 458     }
 459 
 460     /**
 461      * Sets this MutableBigInteger's value array to a copy of the specified
 462      * array. The intLen is set to the length of the specified array.
 463      */
 464     void copyValue(int[] val) {
 465         int len = val.length;
 466         if (value.length < len)
 467             value = new int[len];
 468         System.arraycopy(val, 0, value, 0, len);
 469         intLen = len;
 470         offset = 0;
 471     }
 472 
 473     /**
 474      * Returns true iff this MutableBigInteger has a value of one.
 475      */
 476     boolean isOne() {
 477         return (intLen == 1) && (value[offset] == 1);
 478     }
 479 
 480     /**
 481      * Returns true iff this MutableBigInteger has a value of zero.
 482      */
 483     boolean isZero() {
 484         return (intLen == 0);
 485     }
 486 
 487     /**
 488      * Returns true iff this MutableBigInteger is even.
 489      */
 490     boolean isEven() {
 491         return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
 492     }
 493 
 494     /**
 495      * Returns true iff this MutableBigInteger is odd.
 496      */
 497     boolean isOdd() {
 498         return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
 499     }
 500 
 501     /**
 502      * Returns true iff this MutableBigInteger is in normal form. A
 503      * MutableBigInteger is in normal form if it has no leading zeros
 504      * after the offset, and intLen + offset <= value.length.
 505      */
 506     boolean isNormal() {
 507         if (intLen + offset > value.length)
 508             return false;
 509         if (intLen == 0)
 510             return true;
 511         return (value[offset] != 0);
 512     }
 513 
 514     /**
 515      * Returns a String representation of this MutableBigInteger in radix 10.
 516      */
 517     public String toString() {
 518         BigInteger b = toBigInteger(1);
 519         return b.toString();
 520     }
 521 
 522     /**
 523      * Like {@link #rightShift(int)} but {@code n} can be greater than the length of the number.
 524      */
 525     void safeRightShift(int n) {
 526         if (n/32 >= intLen) {
 527             reset();
 528         } else {
 529             rightShift(n);
 530         }
 531     }
 532 
 533     /**
 534      * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
 535      * in normal form.
 536      */
 537     void rightShift(int n) {
 538         if (intLen == 0)
 539             return;
 540         int nInts = n >>> 5;
 541         int nBits = n & 0x1F;
 542         this.intLen -= nInts;
 543         if (nBits == 0)
 544             return;
 545         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
 546         if (nBits >= bitsInHighWord) {
 547             this.primitiveLeftShift(32 - nBits);
 548             this.intLen--;
 549         } else {
 550             primitiveRightShift(nBits);
 551         }
 552     }
 553 
 554     /**
 555      * Like {@link #leftShift(int)} but {@code n} can be zero.
 556      */
 557     void safeLeftShift(int n) {
 558         if (n > 0) {
 559             leftShift(n);
 560         }
 561     }
 562 
 563     /**
 564      * Left shift this MutableBigInteger n bits.
 565      */
 566     void leftShift(int n) {
 567         /*
 568          * If there is enough storage space in this MutableBigInteger already
 569          * the available space will be used. Space to the right of the used
 570          * ints in the value array is faster to utilize, so the extra space
 571          * will be taken from the right if possible.
 572          */
 573         if (intLen == 0)
 574            return;
 575         int nInts = n >>> 5;
 576         int nBits = n&0x1F;
 577         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
 578 
 579         // If shift can be done without moving words, do so
 580         if (n <= (32-bitsInHighWord)) {
 581             primitiveLeftShift(nBits);
 582             return;
 583         }
 584 
 585         int newLen = intLen + nInts +1;
 586         if (nBits <= (32-bitsInHighWord))
 587             newLen--;
 588         if (value.length < newLen) {
 589             // The array must grow
 590             int[] result = new int[newLen];
 591             for (int i=0; i < intLen; i++)
 592                 result[i] = value[offset+i];
 593             setValue(result, newLen);
 594         } else if (value.length - offset >= newLen) {
 595             // Use space on right
 596             for(int i=0; i < newLen - intLen; i++)
 597                 value[offset+intLen+i] = 0;
 598         } else {
 599             // Must use space on left
 600             for (int i=0; i < intLen; i++)
 601                 value[i] = value[offset+i];
 602             for (int i=intLen; i < newLen; i++)
 603                 value[i] = 0;
 604             offset = 0;
 605         }
 606         intLen = newLen;
 607         if (nBits == 0)
 608             return;
 609         if (nBits <= (32-bitsInHighWord))
 610             primitiveLeftShift(nBits);
 611         else
 612             primitiveRightShift(32 -nBits);
 613     }
 614 
 615     /**
 616      * A primitive used for division. This method adds in one multiple of the
 617      * divisor a back to the dividend result at a specified offset. It is used
 618      * when qhat was estimated too large, and must be adjusted.
 619      */
 620     private int divadd(int[] a, int[] result, int offset) {
 621         long carry = 0;
 622 
 623         for (int j=a.length-1; j >= 0; j--) {
 624             long sum = (a[j] & LONG_MASK) +
 625                        (result[j+offset] & LONG_MASK) + carry;
 626             result[j+offset] = (int)sum;
 627             carry = sum >>> 32;
 628         }
 629         return (int)carry;
 630     }
 631 
 632     /**
 633      * This method is used for division. It multiplies an n word input a by one
 634      * word input x, and subtracts the n word product from q. This is needed
 635      * when subtracting qhat*divisor from dividend.
 636      */
 637     private int mulsub(int[] q, int[] a, int x, int len, int offset) {
 638         long xLong = x & LONG_MASK;
 639         long carry = 0;
 640         offset += len;
 641 
 642         for (int j=len-1; j >= 0; j--) {
 643             long product = (a[j] & LONG_MASK) * xLong + carry;
 644             long difference = q[offset] - product;
 645             q[offset--] = (int)difference;
 646             carry = (product >>> 32)
 647                      + (((difference & LONG_MASK) >
 648                          (((~(int)product) & LONG_MASK))) ? 1:0);
 649         }
 650         return (int)carry;
 651     }
 652 
 653     /**
 654      * The method is the same as mulsun, except the fact that q array is not
 655      * updated, the only result of the method is borrow flag.
 656      */
 657     private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
 658         long xLong = x & LONG_MASK;
 659         long carry = 0;
 660         offset += len;
 661         for (int j=len-1; j >= 0; j--) {
 662             long product = (a[j] & LONG_MASK) * xLong + carry;
 663             long difference = q[offset--] - product;
 664             carry = (product >>> 32)
 665                      + (((difference & LONG_MASK) >
 666                          (((~(int)product) & LONG_MASK))) ? 1:0);
 667         }
 668         return (int)carry;
 669     }
 670 
 671     /**
 672      * Right shift this MutableBigInteger n bits, where n is
 673      * less than 32.
 674      * Assumes that intLen > 0, n > 0 for speed
 675      */
 676     private final void primitiveRightShift(int n) {
 677         int[] val = value;
 678         int n2 = 32 - n;
 679         for (int i=offset+intLen-1, c=val[i]; i > offset; i--) {
 680             int b = c;
 681             c = val[i-1];
 682             val[i] = (c << n2) | (b >>> n);
 683         }
 684         val[offset] >>>= n;
 685     }
 686 
 687     /**
 688      * Left shift this MutableBigInteger n bits, where n is
 689      * less than 32.
 690      * Assumes that intLen > 0, n > 0 for speed
 691      */
 692     private final void primitiveLeftShift(int n) {
 693         int[] val = value;
 694         int n2 = 32 - n;
 695         for (int i=offset, c=val[i], m=i+intLen-1; i < m; i++) {
 696             int b = c;
 697             c = val[i+1];
 698             val[i] = (b << n) | (c >>> n2);
 699         }
 700         val[offset+intLen-1] <<= n;
 701     }
 702 
 703     /**
 704      * Returns a {@code BigInteger} equal to the {@code n}
 705      * low ints of this number.
 706      */
 707     private BigInteger getLower(int n) {
 708         if (isZero()) {
 709             return BigInteger.ZERO;
 710         } else if (intLen < n) {
 711             return toBigInteger(1);
 712         } else {
 713             // strip zeros
 714             int len = n;
 715             while (len > 0 && value[offset+intLen-len] == 0)
 716                 len--;
 717             int sign = len > 0 ? 1 : 0;
 718             return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
 719         }
 720     }
 721 
 722     /**
 723      * Discards all ints whose index is greater than {@code n}.
 724      */
 725     private void keepLower(int n) {
 726         if (intLen >= n) {
 727             offset += intLen - n;
 728             intLen = n;
 729         }
 730     }
 731 
 732     /**
 733      * Adds the contents of two MutableBigInteger objects.The result
 734      * is placed within this MutableBigInteger.
 735      * The contents of the addend are not changed.
 736      */
 737     void add(MutableBigInteger addend) {
 738         int x = intLen;
 739         int y = addend.intLen;
 740         int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
 741         int[] result = (value.length < resultLen ? new int[resultLen] : value);
 742 
 743         int rstart = result.length-1;
 744         long sum;
 745         long carry = 0;
 746 
 747         // Add common parts of both numbers
 748         while(x > 0 && y > 0) {
 749             x--; y--;
 750             sum = (value[x+offset] & LONG_MASK) +
 751                 (addend.value[y+addend.offset] & LONG_MASK) + carry;
 752             result[rstart--] = (int)sum;
 753             carry = sum >>> 32;
 754         }
 755 
 756         // Add remainder of the longer number
 757         while(x > 0) {
 758             x--;
 759             if (carry == 0 && result == value && rstart == (x + offset))
 760                 return;
 761             sum = (value[x+offset] & LONG_MASK) + carry;
 762             result[rstart--] = (int)sum;
 763             carry = sum >>> 32;
 764         }
 765         while(y > 0) {
 766             y--;
 767             sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
 768             result[rstart--] = (int)sum;
 769             carry = sum >>> 32;
 770         }
 771 
 772         if (carry > 0) { // Result must grow in length
 773             resultLen++;
 774             if (result.length < resultLen) {
 775                 int temp[] = new int[resultLen];
 776                 // Result one word longer from carry-out; copy low-order
 777                 // bits into new result.
 778                 System.arraycopy(result, 0, temp, 1, result.length);
 779                 temp[0] = 1;
 780                 result = temp;
 781             } else {
 782                 result[rstart--] = 1;
 783             }
 784         }
 785 
 786         value = result;
 787         intLen = resultLen;
 788         offset = result.length - resultLen;
 789     }
 790 
 791     /**
 792      * Adds the value of {@code addend} shifted {@code n} ints to the left.
 793      * Has the same effect as {@code addend.leftShift(32*ints); add(addend);}
 794      * but doesn't change the value of {@code addend}.
 795      */
 796     void addShifted(MutableBigInteger addend, int n) {
 797         if (addend.isZero()) {
 798             return;
 799         }
 800 
 801         int x = intLen;
 802         int y = addend.intLen + n;
 803         int resultLen = (intLen > y ? intLen : y);
 804         int[] result = (value.length < resultLen ? new int[resultLen] : value);
 805 
 806         int rstart = result.length-1;
 807         long sum;
 808         long carry = 0;
 809 
 810         // Add common parts of both numbers
 811         while (x > 0 && y > 0) {
 812             x--; y--;
 813             int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
 814             sum = (value[x+offset] & LONG_MASK) +
 815                 (bval & LONG_MASK) + carry;
 816             result[rstart--] = (int)sum;
 817             carry = sum >>> 32;
 818         }
 819 
 820         // Add remainder of the longer number
 821         while (x > 0) {
 822             x--;
 823             if (carry == 0 && result == value && rstart == (x + offset)) {
 824                 return;
 825             }
 826             sum = (value[x+offset] & LONG_MASK) + carry;
 827             result[rstart--] = (int)sum;
 828             carry = sum >>> 32;
 829         }
 830         while (y > 0) {
 831             y--;
 832             int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
 833             sum = (bval & LONG_MASK) + carry;
 834             result[rstart--] = (int)sum;
 835             carry = sum >>> 32;
 836         }
 837 
 838         if (carry > 0) { // Result must grow in length
 839             resultLen++;
 840             if (result.length < resultLen) {
 841                 int temp[] = new int[resultLen];
 842                 // Result one word longer from carry-out; copy low-order
 843                 // bits into new result.
 844                 System.arraycopy(result, 0, temp, 1, result.length);
 845                 temp[0] = 1;
 846                 result = temp;
 847             } else {
 848                 result[rstart--] = 1;
 849             }
 850         }
 851 
 852         value = result;
 853         intLen = resultLen;
 854         offset = result.length - resultLen;
 855     }
 856 
 857     /**
 858      * Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must
 859      * not be greater than {@code n}. In other words, concatenates {@code this}
 860      * and {@code addend}.
 861      */
 862     void addDisjoint(MutableBigInteger addend, int n) {
 863         if (addend.isZero())
 864             return;
 865 
 866         int x = intLen;
 867         int y = addend.intLen + n;
 868         int resultLen = (intLen > y ? intLen : y);
 869         int[] result;
 870         if (value.length < resultLen)
 871             result = new int[resultLen];
 872         else {
 873             result = value;
 874             Arrays.fill(value, offset+intLen, value.length, 0);
 875         }
 876 
 877         int rstart = result.length-1;
 878 
 879         // copy from this if needed
 880         System.arraycopy(value, offset, result, rstart+1-x, x);
 881         y -= x;
 882         rstart -= x;
 883 
 884         int len = Math.min(y, addend.value.length-addend.offset);
 885         System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
 886 
 887         // zero the gap
 888         for (int i=rstart+1-y+len; i < rstart+1; i++)
 889             result[i] = 0;
 890 
 891         value = result;
 892         intLen = resultLen;
 893         offset = result.length - resultLen;
 894     }
 895 
 896     /**
 897      * Adds the low {@code n} ints of {@code addend}.
 898      */
 899     void addLower(MutableBigInteger addend, int n) {
 900         MutableBigInteger a = new MutableBigInteger(addend);
 901         if (a.offset + a.intLen >= n) {
 902             a.offset = a.offset + a.intLen - n;
 903             a.intLen = n;
 904         }
 905         a.normalize();
 906         add(a);
 907     }
 908 
 909     /**
 910      * Subtracts the smaller of this and b from the larger and places the
 911      * result into this MutableBigInteger.
 912      */
 913     int subtract(MutableBigInteger b) {
 914         MutableBigInteger a = this;
 915 
 916         int[] result = value;
 917         int sign = a.compare(b);
 918 
 919         if (sign == 0) {
 920             reset();
 921             return 0;
 922         }
 923         if (sign < 0) {
 924             MutableBigInteger tmp = a;
 925             a = b;
 926             b = tmp;
 927         }
 928 
 929         int resultLen = a.intLen;
 930         if (result.length < resultLen)
 931             result = new int[resultLen];
 932 
 933         long diff = 0;
 934         int x = a.intLen;
 935         int y = b.intLen;
 936         int rstart = result.length - 1;
 937 
 938         // Subtract common parts of both numbers
 939         while (y > 0) {
 940             x--; y--;
 941 
 942             diff = (a.value[x+a.offset] & LONG_MASK) -
 943                    (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
 944             result[rstart--] = (int)diff;
 945         }
 946         // Subtract remainder of longer number
 947         while (x > 0) {
 948             x--;
 949             diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
 950             result[rstart--] = (int)diff;
 951         }
 952 
 953         value = result;
 954         intLen = resultLen;
 955         offset = value.length - resultLen;
 956         normalize();
 957         return sign;
 958     }
 959 
 960     /**
 961      * Subtracts the smaller of a and b from the larger and places the result
 962      * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
 963      * operation was performed.
 964      */
 965     private int difference(MutableBigInteger b) {
 966         MutableBigInteger a = this;
 967         int sign = a.compare(b);
 968         if (sign == 0)
 969             return 0;
 970         if (sign < 0) {
 971             MutableBigInteger tmp = a;
 972             a = b;
 973             b = tmp;
 974         }
 975 
 976         long diff = 0;
 977         int x = a.intLen;
 978         int y = b.intLen;
 979 
 980         // Subtract common parts of both numbers
 981         while (y > 0) {
 982             x--; y--;
 983             diff = (a.value[a.offset+ x] & LONG_MASK) -
 984                 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
 985             a.value[a.offset+x] = (int)diff;
 986         }
 987         // Subtract remainder of longer number
 988         while (x > 0) {
 989             x--;
 990             diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
 991             a.value[a.offset+x] = (int)diff;
 992         }
 993 
 994         a.normalize();
 995         return sign;
 996     }
 997 
 998     /**
 999      * Multiply the contents of two MutableBigInteger objects. The result is
1000      * placed into MutableBigInteger z. The contents of y are not changed.
1001      */
1002     void multiply(MutableBigInteger y, MutableBigInteger z) {
1003         int xLen = intLen;
1004         int yLen = y.intLen;
1005         int newLen = xLen + yLen;
1006 
1007         // Put z into an appropriate state to receive product
1008         if (z.value.length < newLen)
1009             z.value = new int[newLen];
1010         z.offset = 0;
1011         z.intLen = newLen;
1012 
1013         // The first iteration is hoisted out of the loop to avoid extra add
1014         long carry = 0;
1015         for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
1016                 long product = (y.value[j+y.offset] & LONG_MASK) *
1017                                (value[xLen-1+offset] & LONG_MASK) + carry;
1018                 z.value[k] = (int)product;
1019                 carry = product >>> 32;
1020         }
1021         z.value[xLen-1] = (int)carry;
1022 
1023         // Perform the multiplication word by word
1024         for (int i = xLen-2; i >= 0; i--) {
1025             carry = 0;
1026             for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
1027                 long product = (y.value[j+y.offset] & LONG_MASK) *
1028                                (value[i+offset] & LONG_MASK) +
1029                                (z.value[k] & LONG_MASK) + carry;
1030                 z.value[k] = (int)product;
1031                 carry = product >>> 32;
1032             }
1033             z.value[i] = (int)carry;
1034         }
1035 
1036         // Remove leading zeros from product
1037         z.normalize();
1038     }
1039 
1040     /**
1041      * Multiply the contents of this MutableBigInteger by the word y. The
1042      * result is placed into z.
1043      */
1044     void mul(int y, MutableBigInteger z) {
1045         if (y == 1) {
1046             z.copyValue(this);
1047             return;
1048         }
1049 
1050         if (y == 0) {
1051             z.clear();
1052             return;
1053         }
1054 
1055         // Perform the multiplication word by word
1056         long ylong = y & LONG_MASK;
1057         int[] zval = (z.value.length < intLen+1 ? new int[intLen + 1]
1058                                               : z.value);
1059         long carry = 0;
1060         for (int i = intLen-1; i >= 0; i--) {
1061             long product = ylong * (value[i+offset] & LONG_MASK) + carry;
1062             zval[i+1] = (int)product;
1063             carry = product >>> 32;
1064         }
1065 
1066         if (carry == 0) {
1067             z.offset = 1;
1068             z.intLen = intLen;
1069         } else {
1070             z.offset = 0;
1071             z.intLen = intLen + 1;
1072             zval[0] = (int)carry;
1073         }
1074         z.value = zval;
1075     }
1076 
1077      /**
1078      * This method is used for division of an n word dividend by a one word
1079      * divisor. The quotient is placed into quotient. The one word divisor is
1080      * specified by divisor.
1081      *
1082      * @return the remainder of the division is returned.
1083      *
1084      */
1085     int divideOneWord(int divisor, MutableBigInteger quotient) {
1086         long divisorLong = divisor & LONG_MASK;
1087 
1088         // Special case of one word dividend
1089         if (intLen == 1) {
1090             long dividendValue = value[offset] & LONG_MASK;
1091             int q = (int) (dividendValue / divisorLong);
1092             int r = (int) (dividendValue - q * divisorLong);
1093             quotient.value[0] = q;
1094             quotient.intLen = (q == 0) ? 0 : 1;
1095             quotient.offset = 0;
1096             return r;
1097         }
1098 
1099         if (quotient.value.length < intLen)
1100             quotient.value = new int[intLen];
1101         quotient.offset = 0;
1102         quotient.intLen = intLen;
1103 
1104         // Normalize the divisor
1105         int shift = Integer.numberOfLeadingZeros(divisor);
1106 
1107         int rem = value[offset];
1108         long remLong = rem & LONG_MASK;
1109         if (remLong < divisorLong) {
1110             quotient.value[0] = 0;
1111         } else {
1112             quotient.value[0] = (int)(remLong / divisorLong);
1113             rem = (int) (remLong - (quotient.value[0] * divisorLong));
1114             remLong = rem & LONG_MASK;
1115         }
1116         int xlen = intLen;
1117         while (--xlen > 0) {
1118             long dividendEstimate = (remLong << 32) |
1119                     (value[offset + intLen - xlen] & LONG_MASK);
1120             int q;
1121             if (dividendEstimate >= 0) {
1122                 q = (int) (dividendEstimate / divisorLong);
1123                 rem = (int) (dividendEstimate - q * divisorLong);
1124             } else {
1125                 long tmp = divWord(dividendEstimate, divisor);
1126                 q = (int) (tmp & LONG_MASK);
1127                 rem = (int) (tmp >>> 32);
1128             }
1129             quotient.value[intLen - xlen] = q;
1130             remLong = rem & LONG_MASK;
1131         }
1132 
1133         quotient.normalize();
1134         // Unnormalize
1135         if (shift > 0)
1136             return rem % divisor;
1137         else
1138             return rem;
1139     }
1140 
1141     /**
1142      * Calculates the quotient of this div b and places the quotient in the
1143      * provided MutableBigInteger objects and the remainder object is returned.
1144      *
1145      */
1146     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
1147         return divide(b,quotient,true);
1148     }
1149 
1150     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1151         if (b.intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD ||
1152                 intLen - b.intLen < BigInteger.BURNIKEL_ZIEGLER_OFFSET) {
1153             return divideKnuth(b, quotient, needRemainder);
1154         } else {
1155             return divideAndRemainderBurnikelZiegler(b, quotient);
1156         }
1157     }
1158 
1159     /**
1160      * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
1161      */
1162     MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
1163         return divideKnuth(b,quotient,true);
1164     }
1165 
1166     /**
1167      * Calculates the quotient of this div b and places the quotient in the
1168      * provided MutableBigInteger objects and the remainder object is returned.
1169      *
1170      * Uses Algorithm D in Knuth section 4.3.1.
1171      * Many optimizations to that algorithm have been adapted from the Colin
1172      * Plumb C library.
1173      * It special cases one word divisors for speed. The content of b is not
1174      * changed.
1175      *
1176      */
1177     MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1178         if (b.intLen == 0)
1179             throw new ArithmeticException("BigInteger divide by zero");
1180 
1181         // Dividend is zero
1182         if (intLen == 0) {
1183             quotient.intLen = quotient.offset = 0;
1184             return needRemainder ? new MutableBigInteger() : null;
1185         }
1186 
1187         int cmp = compare(b);
1188         // Dividend less than divisor
1189         if (cmp < 0) {
1190             quotient.intLen = quotient.offset = 0;
1191             return needRemainder ? new MutableBigInteger(this) : null;
1192         }
1193         // Dividend equal to divisor
1194         if (cmp == 0) {
1195             quotient.value[0] = quotient.intLen = 1;
1196             quotient.offset = 0;
1197             return needRemainder ? new MutableBigInteger() : null;
1198         }
1199 
1200         quotient.clear();
1201         // Special case one word divisor
1202         if (b.intLen == 1) {
1203             int r = divideOneWord(b.value[b.offset], quotient);
1204             if(needRemainder) {
1205                 if (r == 0)
1206                     return new MutableBigInteger();
1207                 return new MutableBigInteger(r);
1208             } else {
1209                 return null;
1210             }
1211         }
1212 
1213         // Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
1214         if (intLen >= KNUTH_POW2_THRESH_LEN) {
1215             int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit());
1216             if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS*32) {
1217                 MutableBigInteger a = new MutableBigInteger(this);
1218                 b = new MutableBigInteger(b);
1219                 a.rightShift(trailingZeroBits);
1220                 b.rightShift(trailingZeroBits);
1221                 MutableBigInteger r = a.divideKnuth(b, quotient);
1222                 r.leftShift(trailingZeroBits);
1223                 return r;
1224             }
1225         }
1226 
1227         return divideMagnitude(b, quotient, needRemainder);
1228     }
1229 
1230     /**
1231      * Computes {@code this/b} and {@code this%b} using the
1232      * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
1233      * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
1234      * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
1235      * multiples of 32 bits.<br/>
1236      * {@code this} and {@code b} must be nonnegative.
1237      * @param b the divisor
1238      * @param quotient output parameter for {@code this/b}
1239      * @return the remainder
1240      */
1241     MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
1242         int r = intLen;
1243         int s = b.intLen;
1244 
1245         // Clear the quotient
1246         quotient.offset = quotient.intLen = 0;
1247 
1248         if (r < s) {
1249             return this;
1250         } else {
1251             // Unlike Knuth division, we don't check for common powers of two here because
1252             // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1253             // additional benefit.
1254 
1255             // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1256             int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1257 
1258             int j = (s+m-1) / m;      // step 2a: j = ceil(s/m)
1259             int n = j * m;            // step 2b: block length in 32-bit units
1260             long n32 = 32L * n;         // block length in bits
1261             int sigma = (int) Math.max(0, n32 - b.bitLength());   // step 3: sigma = max{T | (2^T)*B < beta^n}
1262             MutableBigInteger bShifted = new MutableBigInteger(b);
1263             bShifted.safeLeftShift(sigma);   // step 4a: shift b so its length is a multiple of n
1264             MutableBigInteger aShifted = new MutableBigInteger (this);
1265             aShifted.safeLeftShift(sigma);     // step 4b: shift a by the same amount
1266 
1267             // step 5: t is the number of blocks needed to accommodate a plus one additional bit
1268             int t = (int) ((aShifted.bitLength()+n32) / n32);
1269             if (t < 2) {
1270                 t = 2;
1271             }
1272 
1273             // step 6: conceptually split a into blocks a[t-1], ..., a[0]
1274             MutableBigInteger a1 = aShifted.getBlock(t-1, t, n);   // the most significant block of a
1275 
1276             // step 7: z[t-2] = [a[t-1], a[t-2]]
1277             MutableBigInteger z = aShifted.getBlock(t-2, t, n);    // the second to most significant block
1278             z.addDisjoint(a1, n);   // z[t-2]
1279 
1280             // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
1281             MutableBigInteger qi = new MutableBigInteger();
1282             MutableBigInteger ri;
1283             for (int i=t-2; i > 0; i--) {
1284                 // step 8a: compute (qi,ri) such that z=b*qi+ri
1285                 ri = z.divide2n1n(bShifted, qi);
1286 
1287                 // step 8b: z = [ri, a[i-1]]
1288                 z = aShifted.getBlock(i-1, t, n);   // a[i-1]
1289                 z.addDisjoint(ri, n);
1290                 quotient.addShifted(qi, i*n);   // update q (part of step 9)
1291             }
1292             // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
1293             ri = z.divide2n1n(bShifted, qi);
1294             quotient.add(qi);
1295 
1296             ri.rightShift(sigma);   // step 9: a and b were shifted, so shift back
1297             return ri;
1298         }
1299     }
1300 
1301     /**
1302      * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
1303      * It divides a 2n-digit number by a n-digit number.<br/>
1304      * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
1305      * <br/>
1306      * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
1307      * @param b a positive number such that {@code b.bitLength()} is even
1308      * @param quotient output parameter for {@code this/b}
1309      * @return {@code this%b}
1310      */
1311     private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
1312         int n = b.intLen;
1313 
1314         // step 1: base case
1315         if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1316             return divideKnuth(b, quotient);
1317         }
1318 
1319         // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
1320         MutableBigInteger aUpper = new MutableBigInteger(this);
1321         aUpper.safeRightShift(32*(n/2));   // aUpper = [a1,a2,a3]
1322         keepLower(n/2);   // this = a4
1323 
1324         // step 3: q1=aUpper/b, r1=aUpper%b
1325         MutableBigInteger q1 = new MutableBigInteger();
1326         MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
1327 
1328         // step 4: quotient=[r1,this]/b, r2=[r1,this]%b
1329         addDisjoint(r1, n/2);   // this = [r1,this]
1330         MutableBigInteger r2 = divide3n2n(b, quotient);
1331 
1332         // step 5: let quotient=[q1,quotient] and return r2
1333         quotient.addDisjoint(q1, n/2);
1334         return r2;
1335     }
1336 
1337     /**
1338      * This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
1339      * It divides a 3n-digit number by a 2n-digit number.<br/>
1340      * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.<br/>
1341      * <br/>
1342      * {@code this} must be a nonnegative number such that {@code 2*this.bitLength() <= 3*b.bitLength()}
1343      * @param quotient output parameter for {@code this/b}
1344      * @return {@code this%b}
1345      */
1346     private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
1347         int n = b.intLen / 2;   // half the length of b in ints
1348 
1349         // step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
1350         MutableBigInteger a12 = new MutableBigInteger(this);
1351         a12.safeRightShift(32*n);
1352 
1353         // step 2: view b as [b1,b2] where each bi is n ints or less
1354         MutableBigInteger b1 = new MutableBigInteger(b);
1355         b1.safeRightShift(n * 32);
1356         BigInteger b2 = b.getLower(n);
1357 
1358         MutableBigInteger r;
1359         MutableBigInteger d;
1360         if (compareShifted(b, n) < 0) {
1361             // step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
1362             r = a12.divide2n1n(b1, quotient);
1363 
1364             // step 4: d=quotient*b2
1365             d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));
1366         } else {
1367             // step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
1368             quotient.ones(n);
1369             a12.add(b1);
1370             b1.leftShift(32*n);
1371             a12.subtract(b1);
1372             r = a12;
1373 
1374             // step 4: d=quotient*b2=(b2 << 32*n) - b2
1375             d = new MutableBigInteger(b2);
1376             d.leftShift(32 * n);
1377             d.subtract(new MutableBigInteger(b2));
1378         }
1379 
1380         // step 5: r = r*beta^n + a3 - d (paper says a4)
1381         // However, don't subtract d until after the while loop so r doesn't become negative
1382         r.leftShift(32 * n);
1383         r.addLower(this, n);
1384 
1385         // step 6: add b until r>=d
1386         while (r.compare(d) < 0) {
1387             r.add(b);
1388             quotient.subtract(MutableBigInteger.ONE);
1389         }
1390         r.subtract(d);
1391 
1392         return r;
1393     }
1394 
1395     /**
1396      * Returns a {@code MutableBigInteger} containing {@code blockLength} ints from
1397      * {@code this} number, starting at {@code index*blockLength}.<br/>
1398      * Used by Burnikel-Ziegler division.
1399      * @param index the block index
1400      * @param numBlocks the total number of blocks in {@code this} number
1401      * @param blockLength length of one block in units of 32 bits
1402      * @return
1403      */
1404     private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
1405         int blockStart = index * blockLength;
1406         if (blockStart >= intLen) {
1407             return new MutableBigInteger();
1408         }
1409 
1410         int blockEnd;
1411         if (index == numBlocks-1) {
1412             blockEnd = intLen;
1413         } else {
1414             blockEnd = (index+1) * blockLength;
1415         }
1416         if (blockEnd > intLen) {
1417             return new MutableBigInteger();
1418         }
1419 
1420         int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart);
1421         return new MutableBigInteger(newVal);
1422     }
1423 
1424     /** @see BigInteger#bitLength() */
1425     long bitLength() {
1426         if (intLen == 0)
1427             return 0;
1428         return intLen*32L - Integer.numberOfLeadingZeros(value[offset]);
1429     }
1430 
1431     /**
1432      * Internally used  to calculate the quotient of this div v and places the
1433      * quotient in the provided MutableBigInteger object and the remainder is
1434      * returned.
1435      *
1436      * @return the remainder of the division will be returned.
1437      */
1438     long divide(long v, MutableBigInteger quotient) {
1439         if (v == 0)
1440             throw new ArithmeticException("BigInteger divide by zero");
1441 
1442         // Dividend is zero
1443         if (intLen == 0) {
1444             quotient.intLen = quotient.offset = 0;
1445             return 0;
1446         }
1447         if (v < 0)
1448             v = -v;
1449 
1450         int d = (int)(v >>> 32);
1451         quotient.clear();
1452         // Special case on word divisor
1453         if (d == 0)
1454             return divideOneWord((int)v, quotient) & LONG_MASK;
1455         else {
1456             return divideLongMagnitude(v, quotient).toLong();
1457         }
1458     }
1459 
1460     private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
1461         int n2 = 32 - shift;
1462         int c=src[srcFrom];
1463         for (int i=0; i < srcLen-1; i++) {
1464             int b = c;
1465             c = src[++srcFrom];
1466             dst[dstFrom+i] = (b << shift) | (c >>> n2);
1467         }
1468         dst[dstFrom+srcLen-1] = c << shift;
1469     }
1470 
1471     /**
1472      * Divide this MutableBigInteger by the divisor.
1473      * The quotient will be placed into the provided quotient object &
1474      * the remainder object is returned.
1475      */
1476     private MutableBigInteger divideMagnitude(MutableBigInteger div,
1477                                               MutableBigInteger quotient,
1478                                               boolean needRemainder ) {
1479         // assert div.intLen > 1
1480         // D1 normalize the divisor
1481         int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
1482         // Copy divisor value to protect divisor
1483         final int dlen = div.intLen;
1484         int[] divisor;
1485         MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
1486         if (shift > 0) {
1487             divisor = new int[dlen];
1488             copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
1489             if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
1490                 int[] remarr = new int[intLen + 1];
1491                 rem = new MutableBigInteger(remarr);
1492                 rem.intLen = intLen;
1493                 rem.offset = 1;
1494                 copyAndShift(value,offset,intLen,remarr,1,shift);
1495             } else {
1496                 int[] remarr = new int[intLen + 2];
1497                 rem = new MutableBigInteger(remarr);
1498                 rem.intLen = intLen+1;
1499                 rem.offset = 1;
1500                 int rFrom = offset;
1501                 int c=0;
1502                 int n2 = 32 - shift;
1503                 for (int i=1; i < intLen+1; i++,rFrom++) {
1504                     int b = c;
1505                     c = value[rFrom];
1506                     remarr[i] = (b << shift) | (c >>> n2);
1507                 }
1508                 remarr[intLen+1] = c << shift;
1509             }
1510         } else {
1511             divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
1512             rem = new MutableBigInteger(new int[intLen + 1]);
1513             System.arraycopy(value, offset, rem.value, 1, intLen);
1514             rem.intLen = intLen;
1515             rem.offset = 1;
1516         }
1517 
1518         int nlen = rem.intLen;
1519 
1520         // Set the quotient size
1521         final int limit = nlen - dlen + 1;
1522         if (quotient.value.length < limit) {
1523             quotient.value = new int[limit];
1524             quotient.offset = 0;
1525         }
1526         quotient.intLen = limit;
1527         int[] q = quotient.value;
1528 
1529 
1530         // Must insert leading 0 in rem if its length did not change
1531         if (rem.intLen == nlen) {
1532             rem.offset = 0;
1533             rem.value[0] = 0;
1534             rem.intLen++;
1535         }
1536 
1537         int dh = divisor[0];
1538         long dhLong = dh & LONG_MASK;
1539         int dl = divisor[1];
1540 
1541         // D2 Initialize j
1542         for (int j=0; j < limit-1; j++) {
1543             // D3 Calculate qhat
1544             // estimate qhat
1545             int qhat = 0;
1546             int qrem = 0;
1547             boolean skipCorrection = false;
1548             int nh = rem.value[j+rem.offset];
1549             int nh2 = nh + 0x80000000;
1550             int nm = rem.value[j+1+rem.offset];
1551 
1552             if (nh == dh) {
1553                 qhat = ~0;
1554                 qrem = nh + nm;
1555                 skipCorrection = qrem + 0x80000000 < nh2;
1556             } else {
1557                 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1558                 if (nChunk >= 0) {
1559                     qhat = (int) (nChunk / dhLong);
1560                     qrem = (int) (nChunk - (qhat * dhLong));
1561                 } else {
1562                     long tmp = divWord(nChunk, dh);
1563                     qhat = (int) (tmp & LONG_MASK);
1564                     qrem = (int) (tmp >>> 32);
1565                 }
1566             }
1567 
1568             if (qhat == 0)
1569                 continue;
1570 
1571             if (!skipCorrection) { // Correct qhat
1572                 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1573                 long rs = ((qrem & LONG_MASK) << 32) | nl;
1574                 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1575 
1576                 if (unsignedLongCompare(estProduct, rs)) {
1577                     qhat--;
1578                     qrem = (int)((qrem & LONG_MASK) + dhLong);
1579                     if ((qrem & LONG_MASK) >=  dhLong) {
1580                         estProduct -= (dl & LONG_MASK);
1581                         rs = ((qrem & LONG_MASK) << 32) | nl;
1582                         if (unsignedLongCompare(estProduct, rs))
1583                             qhat--;
1584                     }
1585                 }
1586             }
1587 
1588             // D4 Multiply and subtract
1589             rem.value[j+rem.offset] = 0;
1590             int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1591 
1592             // D5 Test remainder
1593             if (borrow + 0x80000000 > nh2) {
1594                 // D6 Add back
1595                 divadd(divisor, rem.value, j+1+rem.offset);
1596                 qhat--;
1597             }
1598 
1599             // Store the quotient digit
1600             q[j] = qhat;
1601         } // D7 loop on j
1602         // D3 Calculate qhat
1603         // estimate qhat
1604         int qhat = 0;
1605         int qrem = 0;
1606         boolean skipCorrection = false;
1607         int nh = rem.value[limit - 1 + rem.offset];
1608         int nh2 = nh + 0x80000000;
1609         int nm = rem.value[limit + rem.offset];
1610 
1611         if (nh == dh) {
1612             qhat = ~0;
1613             qrem = nh + nm;
1614             skipCorrection = qrem + 0x80000000 < nh2;
1615         } else {
1616             long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1617             if (nChunk >= 0) {
1618                 qhat = (int) (nChunk / dhLong);
1619                 qrem = (int) (nChunk - (qhat * dhLong));
1620             } else {
1621                 long tmp = divWord(nChunk, dh);
1622                 qhat = (int) (tmp & LONG_MASK);
1623                 qrem = (int) (tmp >>> 32);
1624             }
1625         }
1626         if (qhat != 0) {
1627             if (!skipCorrection) { // Correct qhat
1628                 long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
1629                 long rs = ((qrem & LONG_MASK) << 32) | nl;
1630                 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1631 
1632                 if (unsignedLongCompare(estProduct, rs)) {
1633                     qhat--;
1634                     qrem = (int) ((qrem & LONG_MASK) + dhLong);
1635                     if ((qrem & LONG_MASK) >= dhLong) {
1636                         estProduct -= (dl & LONG_MASK);
1637                         rs = ((qrem & LONG_MASK) << 32) | nl;
1638                         if (unsignedLongCompare(estProduct, rs))
1639                             qhat--;
1640                     }
1641                 }
1642             }
1643 
1644 
1645             // D4 Multiply and subtract
1646             int borrow;
1647             rem.value[limit - 1 + rem.offset] = 0;
1648             if(needRemainder)
1649                 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1650             else
1651                 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1652 
1653             // D5 Test remainder
1654             if (borrow + 0x80000000 > nh2) {
1655                 // D6 Add back
1656                 if(needRemainder)
1657                     divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
1658                 qhat--;
1659             }
1660 
1661             // Store the quotient digit
1662             q[(limit - 1)] = qhat;
1663         }
1664 
1665 
1666         if (needRemainder) {
1667             // D8 Unnormalize
1668             if (shift > 0)
1669                 rem.rightShift(shift);
1670             rem.normalize();
1671         }
1672         quotient.normalize();
1673         return needRemainder ? rem : null;
1674     }
1675 
1676     /**
1677      * Divide this MutableBigInteger by the divisor represented by positive long
1678      * value. The quotient will be placed into the provided quotient object &
1679      * the remainder object is returned.
1680      */
1681     private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
1682         // Remainder starts as dividend with space for a leading zero
1683         MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1684         System.arraycopy(value, offset, rem.value, 1, intLen);
1685         rem.intLen = intLen;
1686         rem.offset = 1;
1687 
1688         int nlen = rem.intLen;
1689 
1690         int limit = nlen - 2 + 1;
1691         if (quotient.value.length < limit) {
1692             quotient.value = new int[limit];
1693             quotient.offset = 0;
1694         }
1695         quotient.intLen = limit;
1696         int[] q = quotient.value;
1697 
1698         // D1 normalize the divisor
1699         int shift = Long.numberOfLeadingZeros(ldivisor);
1700         if (shift > 0) {
1701             ldivisor<<=shift;
1702             rem.leftShift(shift);
1703         }
1704 
1705         // Must insert leading 0 in rem if its length did not change
1706         if (rem.intLen == nlen) {
1707             rem.offset = 0;
1708             rem.value[0] = 0;
1709             rem.intLen++;
1710         }
1711 
1712         int dh = (int)(ldivisor >>> 32);
1713         long dhLong = dh & LONG_MASK;
1714         int dl = (int)(ldivisor & LONG_MASK);
1715 
1716         // D2 Initialize j
1717         for (int j = 0; j < limit; j++) {
1718             // D3 Calculate qhat
1719             // estimate qhat
1720             int qhat = 0;
1721             int qrem = 0;
1722             boolean skipCorrection = false;
1723             int nh = rem.value[j + rem.offset];
1724             int nh2 = nh + 0x80000000;
1725             int nm = rem.value[j + 1 + rem.offset];
1726 
1727             if (nh == dh) {
1728                 qhat = ~0;
1729                 qrem = nh + nm;
1730                 skipCorrection = qrem + 0x80000000 < nh2;
1731             } else {
1732                 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1733                 if (nChunk >= 0) {
1734                     qhat = (int) (nChunk / dhLong);
1735                     qrem = (int) (nChunk - (qhat * dhLong));
1736                 } else {
1737                     long tmp = divWord(nChunk, dh);
1738                     qhat =(int)(tmp & LONG_MASK);
1739                     qrem = (int)(tmp>>>32);
1740                 }
1741             }
1742 
1743             if (qhat == 0)
1744                 continue;
1745 
1746             if (!skipCorrection) { // Correct qhat
1747                 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
1748                 long rs = ((qrem & LONG_MASK) << 32) | nl;
1749                 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1750 
1751                 if (unsignedLongCompare(estProduct, rs)) {
1752                     qhat--;
1753                     qrem = (int) ((qrem & LONG_MASK) + dhLong);
1754                     if ((qrem & LONG_MASK) >= dhLong) {
1755                         estProduct -= (dl & LONG_MASK);
1756                         rs = ((qrem & LONG_MASK) << 32) | nl;
1757                         if (unsignedLongCompare(estProduct, rs))
1758                             qhat--;
1759                     }
1760                 }
1761             }
1762 
1763             // D4 Multiply and subtract
1764             rem.value[j + rem.offset] = 0;
1765             int borrow = mulsubLong(rem.value, dh, dl, qhat,  j + rem.offset);
1766 
1767             // D5 Test remainder
1768             if (borrow + 0x80000000 > nh2) {
1769                 // D6 Add back
1770                 divaddLong(dh,dl, rem.value, j + 1 + rem.offset);
1771                 qhat--;
1772             }
1773 
1774             // Store the quotient digit
1775             q[j] = qhat;
1776         } // D7 loop on j
1777 
1778         // D8 Unnormalize
1779         if (shift > 0)
1780             rem.rightShift(shift);
1781 
1782         quotient.normalize();
1783         rem.normalize();
1784         return rem;
1785     }
1786 
1787     /**
1788      * A primitive used for division by long.
1789      * Specialized version of the method divadd.
1790      * dh is a high part of the divisor, dl is a low part
1791      */
1792     private int divaddLong(int dh, int dl, int[] result, int offset) {
1793         long carry = 0;
1794 
1795         long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK);
1796         result[1+offset] = (int)sum;
1797 
1798         sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
1799         result[offset] = (int)sum;
1800         carry = sum >>> 32;
1801         return (int)carry;
1802     }
1803 
1804     /**
1805      * This method is used for division by long.
1806      * Specialized version of the method sulsub.
1807      * dh is a high part of the divisor, dl is a low part
1808      */
1809     private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
1810         long xLong = x & LONG_MASK;
1811         offset += 2;
1812         long product = (dl & LONG_MASK) * xLong;
1813         long difference = q[offset] - product;
1814         q[offset--] = (int)difference;
1815         long carry = (product >>> 32)
1816                  + (((difference & LONG_MASK) >
1817                      (((~(int)product) & LONG_MASK))) ? 1:0);
1818         product = (dh & LONG_MASK) * xLong + carry;
1819         difference = q[offset] - product;
1820         q[offset--] = (int)difference;
1821         carry = (product >>> 32)
1822                  + (((difference & LONG_MASK) >
1823                      (((~(int)product) & LONG_MASK))) ? 1:0);
1824         return (int)carry;
1825     }
1826 
1827     /**
1828      * Compare two longs as if they were unsigned.
1829      * Returns true iff one is bigger than two.
1830      */
1831     private boolean unsignedLongCompare(long one, long two) {
1832         return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1833     }
1834 
1835     /**
1836      * This method divides a long quantity by an int to estimate
1837      * qhat for two multi precision numbers. It is used when
1838      * the signed value of n is less than zero.
1839      * Returns long value where high 32 bits contain remainder value and
1840      * low 32 bits contain quotient value.
1841      */
1842     static long divWord(long n, int d) {
1843         long dLong = d & LONG_MASK;
1844         long r;
1845         long q;
1846         if (dLong == 1) {
1847             q = (int)n;
1848             r = 0;
1849             return (r << 32) | (q & LONG_MASK);
1850         }
1851 
1852         // Approximate the quotient and remainder
1853         q = (n >>> 1) / (dLong >>> 1);
1854         r = n - q*dLong;
1855 
1856         // Correct the approximation
1857         while (r < 0) {
1858             r += dLong;
1859             q--;
1860         }
1861         while (r >= dLong) {
1862             r -= dLong;
1863             q++;
1864         }
1865         // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1866         return (r << 32) | (q & LONG_MASK);
1867     }
1868 
1869     /**
1870      * Calculate GCD of this and b. This and b are changed by the computation.
1871      */
1872     MutableBigInteger hybridGCD(MutableBigInteger b) {
1873         // Use Euclid's algorithm until the numbers are approximately the
1874         // same length, then use the binary GCD algorithm to find the GCD.
1875         MutableBigInteger a = this;
1876         MutableBigInteger q = new MutableBigInteger();
1877 
1878         while (b.intLen != 0) {
1879             if (Math.abs(a.intLen - b.intLen) < 2)
1880                 return a.binaryGCD(b);
1881 
1882             MutableBigInteger r = a.divide(b, q);
1883             a = b;
1884             b = r;
1885         }
1886         return a;
1887     }
1888 
1889     /**
1890      * Calculate GCD of this and v.
1891      * Assumes that this and v are not zero.
1892      */
1893     private MutableBigInteger binaryGCD(MutableBigInteger v) {
1894         // Algorithm B from Knuth section 4.5.2
1895         MutableBigInteger u = this;
1896         MutableBigInteger r = new MutableBigInteger();
1897 
1898         // step B1
1899         int s1 = u.getLowestSetBit();
1900         int s2 = v.getLowestSetBit();
1901         int k = (s1 < s2) ? s1 : s2;
1902         if (k != 0) {
1903             u.rightShift(k);
1904             v.rightShift(k);
1905         }
1906 
1907         // step B2
1908         boolean uOdd = (k == s1);
1909         MutableBigInteger t = uOdd ? v: u;
1910         int tsign = uOdd ? -1 : 1;
1911 
1912         int lb;
1913         while ((lb = t.getLowestSetBit()) >= 0) {
1914             // steps B3 and B4
1915             t.rightShift(lb);
1916             // step B5
1917             if (tsign > 0)
1918                 u = t;
1919             else
1920                 v = t;
1921 
1922             // Special case one word numbers
1923             if (u.intLen < 2 && v.intLen < 2) {
1924                 int x = u.value[u.offset];
1925                 int y = v.value[v.offset];
1926                 x  = binaryGcd(x, y);
1927                 r.value[0] = x;
1928                 r.intLen = 1;
1929                 r.offset = 0;
1930                 if (k > 0)
1931                     r.leftShift(k);
1932                 return r;
1933             }
1934 
1935             // step B6
1936             if ((tsign = u.difference(v)) == 0)
1937                 break;
1938             t = (tsign >= 0) ? u : v;
1939         }
1940 
1941         if (k > 0)
1942             u.leftShift(k);
1943         return u;
1944     }
1945 
1946     /**
1947      * Calculate GCD of a and b interpreted as unsigned integers.
1948      */
1949     static int binaryGcd(int a, int b) {
1950         if (b == 0)
1951             return a;
1952         if (a == 0)
1953             return b;
1954 
1955         // Right shift a & b till their last bits equal to 1.
1956         int aZeros = Integer.numberOfTrailingZeros(a);
1957         int bZeros = Integer.numberOfTrailingZeros(b);
1958         a >>>= aZeros;
1959         b >>>= bZeros;
1960 
1961         int t = (aZeros < bZeros ? aZeros : bZeros);
1962 
1963         while (a != b) {
1964             if ((a+0x80000000) > (b+0x80000000)) {  // a > b as unsigned
1965                 a -= b;
1966                 a >>>= Integer.numberOfTrailingZeros(a);
1967             } else {
1968                 b -= a;
1969                 b >>>= Integer.numberOfTrailingZeros(b);
1970             }
1971         }
1972         return a<<t;
1973     }
1974 
1975     /**
1976      * Returns the modInverse of this mod p. This and p are not affected by
1977      * the operation.
1978      */
1979     MutableBigInteger mutableModInverse(MutableBigInteger p) {
1980         // Modulus is odd, use Schroeppel's algorithm
1981         if (p.isOdd())
1982             return modInverse(p);
1983 
1984         // Base and modulus are even, throw exception
1985         if (isEven())
1986             throw new ArithmeticException("BigInteger not invertible.");
1987 
1988         // Get even part of modulus expressed as a power of 2
1989         int powersOf2 = p.getLowestSetBit();
1990 
1991         // Construct odd part of modulus
1992         MutableBigInteger oddMod = new MutableBigInteger(p);
1993         oddMod.rightShift(powersOf2);
1994 
1995         if (oddMod.isOne())
1996             return modInverseMP2(powersOf2);
1997 
1998         // Calculate 1/a mod oddMod
1999         MutableBigInteger oddPart = modInverse(oddMod);
2000 
2001         // Calculate 1/a mod evenMod
2002         MutableBigInteger evenPart = modInverseMP2(powersOf2);
2003 
2004         // Combine the results using Chinese Remainder Theorem
2005         MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
2006         MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
2007 
2008         MutableBigInteger temp1 = new MutableBigInteger();
2009         MutableBigInteger temp2 = new MutableBigInteger();
2010         MutableBigInteger result = new MutableBigInteger();
2011 
2012         oddPart.leftShift(powersOf2);
2013         oddPart.multiply(y1, result);
2014 
2015         evenPart.multiply(oddMod, temp1);
2016         temp1.multiply(y2, temp2);
2017 
2018         result.add(temp2);
2019         return result.divide(p, temp1);
2020     }
2021 
2022     /*
2023      * Calculate the multiplicative inverse of this mod 2^k.
2024      */
2025     MutableBigInteger modInverseMP2(int k) {
2026         if (isEven())
2027             throw new ArithmeticException("Non-invertible. (GCD != 1)");
2028 
2029         if (k > 64)
2030             return euclidModInverse(k);
2031 
2032         int t = inverseMod32(value[offset+intLen-1]);
2033 
2034         if (k < 33) {
2035             t = (k == 32 ? t : t & ((1 << k) - 1));
2036             return new MutableBigInteger(t);
2037         }
2038 
2039         long pLong = (value[offset+intLen-1] & LONG_MASK);
2040         if (intLen > 1)
2041             pLong |=  ((long)value[offset+intLen-2] << 32);
2042         long tLong = t & LONG_MASK;
2043         tLong = tLong * (2 - pLong * tLong);  // 1 more Newton iter step
2044         tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
2045 
2046         MutableBigInteger result = new MutableBigInteger(new int[2]);
2047         result.value[0] = (int)(tLong >>> 32);
2048         result.value[1] = (int)tLong;
2049         result.intLen = 2;
2050         result.normalize();
2051         return result;
2052     }
2053 
2054     /**
2055      * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.
2056      */
2057     static int inverseMod32(int val) {
2058         // Newton's iteration!
2059         int t = val;
2060         t *= 2 - val*t;
2061         t *= 2 - val*t;
2062         t *= 2 - val*t;
2063         t *= 2 - val*t;
2064         return t;
2065     }
2066 
2067     /**
2068      * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
2069      */
2070     static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
2071         // Copy the mod to protect original
2072         return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
2073     }
2074 
2075     /**
2076      * Calculate the multiplicative inverse of this mod mod, where mod is odd.
2077      * This and mod are not changed by the calculation.
2078      *
2079      * This method implements an algorithm due to Richard Schroeppel, that uses
2080      * the same intermediate representation as Montgomery Reduction
2081      * ("Montgomery Form").  The algorithm is described in an unpublished
2082      * manuscript entitled "Fast Modular Reciprocals."
2083      */
2084     private MutableBigInteger modInverse(MutableBigInteger mod) {
2085         MutableBigInteger p = new MutableBigInteger(mod);
2086         MutableBigInteger f = new MutableBigInteger(this);
2087         MutableBigInteger g = new MutableBigInteger(p);
2088         SignedMutableBigInteger c = new SignedMutableBigInteger(1);
2089         SignedMutableBigInteger d = new SignedMutableBigInteger();
2090         MutableBigInteger temp = null;
2091         SignedMutableBigInteger sTemp = null;
2092 
2093         int k = 0;
2094         // Right shift f k times until odd, left shift d k times
2095         if (f.isEven()) {
2096             int trailingZeros = f.getLowestSetBit();
2097             f.rightShift(trailingZeros);
2098             d.leftShift(trailingZeros);
2099             k = trailingZeros;
2100         }
2101 
2102         // The Almost Inverse Algorithm
2103         while (!f.isOne()) {
2104             // If gcd(f, g) != 1, number is not invertible modulo mod
2105             if (f.isZero())
2106                 throw new ArithmeticException("BigInteger not invertible.");
2107 
2108             // If f < g exchange f, g and c, d
2109             if (f.compare(g) < 0) {
2110                 temp = f; f = g; g = temp;
2111                 sTemp = d; d = c; c = sTemp;
2112             }
2113 
2114             // If f == g (mod 4)
2115             if (((f.value[f.offset + f.intLen - 1] ^
2116                  g.value[g.offset + g.intLen - 1]) & 3) == 0) {
2117                 f.subtract(g);
2118                 c.signedSubtract(d);
2119             } else { // If f != g (mod 4)
2120                 f.add(g);
2121                 c.signedAdd(d);
2122             }
2123 
2124             // Right shift f k times until odd, left shift d k times
2125             int trailingZeros = f.getLowestSetBit();
2126             f.rightShift(trailingZeros);
2127             d.leftShift(trailingZeros);
2128             k += trailingZeros;
2129         }
2130 
2131         while (c.sign < 0)
2132            c.signedAdd(p);
2133 
2134         return fixup(c, p, k);
2135     }
2136 
2137     /**
2138      * The Fixup Algorithm
2139      * Calculates X such that X = C * 2^(-k) (mod P)
2140      * Assumes C<P and P is odd.
2141      */
2142     static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
2143                                                                       int k) {
2144         MutableBigInteger temp = new MutableBigInteger();
2145         // Set r to the multiplicative inverse of p mod 2^32
2146         int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
2147 
2148         for (int i=0, numWords = k >> 5; i < numWords; i++) {
2149             // V = R * c (mod 2^j)
2150             int  v = r * c.value[c.offset + c.intLen-1];
2151             // c = c + (v * p)
2152             p.mul(v, temp);
2153             c.add(temp);
2154             // c = c / 2^j
2155             c.intLen--;
2156         }
2157         int numBits = k & 0x1f;
2158         if (numBits != 0) {
2159             // V = R * c (mod 2^j)
2160             int v = r * c.value[c.offset + c.intLen-1];
2161             v &= ((1<<numBits) - 1);
2162             // c = c + (v * p)
2163             p.mul(v, temp);
2164             c.add(temp);
2165             // c = c / 2^j
2166             c.rightShift(numBits);
2167         }
2168 
2169         // In theory, c may be greater than p at this point (Very rare!)
2170         while (c.compare(p) >= 0)
2171             c.subtract(p);
2172 
2173         return c;
2174     }
2175 
2176     /**
2177      * Uses the extended Euclidean algorithm to compute the modInverse of base
2178      * mod a modulus that is a power of 2. The modulus is 2^k.
2179      */
2180     MutableBigInteger euclidModInverse(int k) {
2181         MutableBigInteger b = new MutableBigInteger(1);
2182         b.leftShift(k);
2183         MutableBigInteger mod = new MutableBigInteger(b);
2184 
2185         MutableBigInteger a = new MutableBigInteger(this);
2186         MutableBigInteger q = new MutableBigInteger();
2187         MutableBigInteger r = b.divide(a, q);
2188 
2189         MutableBigInteger swapper = b;
2190         // swap b & r
2191         b = r;
2192         r = swapper;
2193 
2194         MutableBigInteger t1 = new MutableBigInteger(q);
2195         MutableBigInteger t0 = new MutableBigInteger(1);
2196         MutableBigInteger temp = new MutableBigInteger();
2197 
2198         while (!b.isOne()) {
2199             r = a.divide(b, q);
2200 
2201             if (r.intLen == 0)
2202                 throw new ArithmeticException("BigInteger not invertible.");
2203 
2204             swapper = r;
2205             a = swapper;
2206 
2207             if (q.intLen == 1)
2208                 t1.mul(q.value[q.offset], temp);
2209             else
2210                 q.multiply(t1, temp);
2211             swapper = q;
2212             q = temp;
2213             temp = swapper;
2214             t0.add(q);
2215 
2216             if (a.isOne())
2217                 return t0;
2218 
2219             r = b.divide(a, q);
2220 
2221             if (r.intLen == 0)
2222                 throw new ArithmeticException("BigInteger not invertible.");
2223 
2224             swapper = b;
2225             b =  r;
2226 
2227             if (q.intLen == 1)
2228                 t0.mul(q.value[q.offset], temp);
2229             else
2230                 q.multiply(t0, temp);
2231             swapper = q; q = temp; temp = swapper;
2232 
2233             t1.add(q);
2234         }
2235         mod.subtract(t1);
2236         return mod;
2237     }
2238 }