rem.value[j+rem.offset] = 0;            int borrow = mulsub(rem.value, d, qhat, dlen, j+rem.offset);            // D5 Test remainder            if (borrow + 0x80000000 > nh2) {                // D6 Add back                divadd(d, rem.value, j+1+rem.offset);                qhat--;            }            // Store the quotient digit            q[j] = qhat;        } // D7 loop on j        // D8 Unnormalize        if (shift > 0)            rem.rightShift(shift);        rem.normalize();        quotient.normalize();    }    /**     * Compare two longs as if they were unsigned.     * Returns true iff one is bigger than two.     */    private boolean unsignedLongCompare(long one, long two) {        return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);    }    /**     * This method divides a long quantity by an int to estimate     * qhat for two multi precision numbers. It is used when      * the signed value of n is less than zero.     */    private void divWord(int[] result, long n, int d) {        long dLong = d & LONG_MASK;        if (dLong == 1) {            result[0] = (int)n;            result[1] = 0;            return;        }                        // Approximate the quotient and remainder        long q = (n >>> 1) / (dLong >>> 1);        long r = n - q*dLong;        // Correct the approximation        while (r < 0) {            r += dLong;            q--;        }        while (r >= dLong) {            r -= dLong;            q++;        }        // n - q*dlong == r && 0 <= r <dLong, hence we're done.        result[0] = (int)q;        result[1] = (int)r;    }    /**     * Calculate GCD of this and b. This and b are changed by the computation.     */    MutableBigInteger hybridGCD(MutableBigInteger b) {        // Use Euclid's algorithm until the numbers are approximately the        // same length, then use the binary GCD algorithm to find the GCD.        MutableBigInteger a = this;        MutableBigInteger q = new MutableBigInteger(),                          r = new MutableBigInteger();        while (b.intLen != 0) {            if (Math.abs(a.intLen - b.intLen) < 2)                return a.binaryGCD(b);            a.divide(b, q, r);            MutableBigInteger swapper = a;            a = b; b = r; r = swapper;        }        return a;    }    /**     * Calculate GCD of this and v.     * Assumes that this and v are not zero.     */    private MutableBigInteger binaryGCD(MutableBigInteger v) {        // Algorithm B from Knuth section 4.5.2        MutableBigInteger u = this;        MutableBigInteger q = new MutableBigInteger(),            r = new MutableBigInteger();        // step B1        int s1 = u.getLowestSetBit();        int s2 = v.getLowestSetBit();        int k = (s1 < s2) ? s1 : s2;        if (k != 0) {            u.rightShift(k);            v.rightShift(k);        }        // step B2        boolean uOdd = (k==s1);        MutableBigInteger t = uOdd ? v: u;        int tsign = uOdd ? -1 : 1;        int lb;        while ((lb = t.getLowestSetBit()) >= 0) {            // steps B3 and B4            t.rightShift(lb);            // step B5            if (tsign > 0)                u = t;            else                v = t;            // Special case one word numbers            if (u.intLen < 2 && v.intLen < 2) {                int x = u.value[u.offset];                int y = v.value[v.offset];                x  = binaryGcd(x, y);                r.value[0] = x;                r.intLen = 1;                r.offset = 0;                if (k > 0)                    r.leftShift(k);                return r;            }                            // step B6            if ((tsign = u.difference(v)) == 0)                break;            t = (tsign >= 0) ? u : v;        }        if (k > 0)            u.leftShift(k);        return u;    }    /**     * Calculate GCD of a and b interpreted as unsigned integers.     */    static int binaryGcd(int a, int b) {        if (b==0)            return a;        if (a==0)            return b;        int x;        int aZeros = 0;        while ((x = (int)a & 0xff) == 0) {            a >>>= 8;            aZeros += 8;        }        int y = BigInteger.trailingZeroTable[x];        aZeros += y;        a >>>= y;        int bZeros = 0;        while ((x = (int)b & 0xff) == 0) {            b >>>= 8;            bZeros += 8;        }        y = BigInteger.trailingZeroTable[x];        bZeros += y;        b >>>= y;        int t = (aZeros < bZeros ? aZeros : bZeros);        while (a != b) {            if ((a+0x80000000) > (b+0x80000000)) {  // a > b as unsigned                a -= b;                while ((x = (int)a & 0xff) == 0)                    a >>>= 8;                a >>>= BigInteger.trailingZeroTable[x];            } else {                b -= a;                while ((x = (int)b & 0xff) == 0)                    b >>>= 8;                b >>>= BigInteger.trailingZeroTable[x];            }        }        return a<<t;    }    /**     * Returns the modInverse of this mod p. This and p are not affected by     * the operation.     */    MutableBigInteger mutableModInverse(MutableBigInteger p) {        // Modulus is odd, use Schroeppel's algorithm        if (p.isOdd())            return modInverse(p);        // Base and modulus are even, throw exception        if (isEven())            throw new ArithmeticException("BigInteger not invertible.");        // Get even part of modulus expressed as a power of 2        int powersOf2 = p.getLowestSetBit();        // Construct odd part of modulus        MutableBigInteger oddMod = new MutableBigInteger(p);        oddMod.rightShift(powersOf2);        if (oddMod.isOne())            return modInverseMP2(powersOf2);        // Calculate 1/a mod oddMod        MutableBigInteger oddPart = modInverse(oddMod);        // Calculate 1/a mod evenMod        MutableBigInteger evenPart = modInverseMP2(powersOf2);        // Combine the results using Chinese Remainder Theorem        MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);        MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);        MutableBigInteger temp1 = new MutableBigInteger();        MutableBigInteger temp2 = new MutableBigInteger();        MutableBigInteger result = new MutableBigInteger();        oddPart.leftShift(powersOf2);        oddPart.multiply(y1, result);        evenPart.multiply(oddMod, temp1);        temp1.multiply(y2, temp2);        result.add(temp2);        result.divide(p, temp1, temp2);        return temp2;    }    /*     * Calculate the multiplicative inverse of this mod 2^k.     */    MutableBigInteger modInverseMP2(int k) {        if (isEven())            throw new ArithmeticException("Non-invertible. (GCD != 1)");        if (k > 64)            return euclidModInverse(k);        int t = inverseMod32(value[offset+intLen-1]);        if (k < 33) {            t = (k == 32 ? t : t & ((1 << k) - 1));            return new MutableBigInteger(t);        }                   long pLong = (value[offset+intLen-1] & LONG_MASK);        if (intLen > 1)            pLong |=  ((long)value[offset+intLen-2] << 32);        long tLong = t & LONG_MASK;        tLong = tLong * (2 - pLong * tLong);  // 1 more Newton iter step        tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));        MutableBigInteger result = new MutableBigInteger(new int[2]);        result.value[0] = (int)(tLong >>> 32);        result.value[1] = (int)tLong;        result.intLen = 2;        result.normalize();        return result;       }    /*     * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.     */    static int inverseMod32(int val) {        // Newton's iteration!        int t = val;        t *= 2 - val*t;        t *= 2 - val*t;        t *= 2 - val*t;        t *= 2 - val*t;        return t;    }    /*     * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.     */    static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {        // Copy the mod to protect original        return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);    }    /**     * Calculate the multiplicative inverse of this mod mod, where mod is odd.     * This and mod are not changed by the calculation.     *     * This method implements an algorithm due to Richard Schroeppel, that uses     * the same intermediate representation as Montgomery Reduction     * ("Montgomery Form").  The algorithm is described in an unpublished     * manuscript entitled "Fast Modular Reciprocals."     */    private MutableBigInteger modInverse(MutableBigInteger mod) {        MutableBigInteger p = new MutableBigInteger(mod);        MutableBigInteger f = new MutableBigInteger(this);        MutableBigInteger g = new MutableBigInteger(p);        SignedMutableBigInteger c = new SignedMutableBigInteger(1);        SignedMutableBigInteger d = new SignedMutableBigInteger();        MutableBigInteger temp = null;        SignedMutableBigInteger sTemp = null;        int k = 0;        // Right shift f k times until odd, left shift d k times        if (f.isEven()) {            int trailingZeros = f.getLowestSetBit();            f.rightShift(trailingZeros);            d.leftShift(trailingZeros);            k = trailingZeros;        }                // The Almost Inverse Algorithm        while(!f.isOne()) {             // If gcd(f, g) != 1, number is not invertible modulo mod            if (f.isZero())                throw new ArithmeticException("BigInteger not invertible.");            // If f < g exchange f, g and c, d            if (f.compare(g) < 0) {                temp = f; f = g; g = temp;                sTemp = d; d = c; c = sTemp;            }            // If f == g (mod 4)             if (((f.value[f.offset + f.intLen - 1] ^                 g.value[g.offset + g.intLen - 1]) & 3) == 0) {                f.subtract(g);                c.signedSubtract(d);            } else { // If f != g (mod 4)                f.add(g);                c.signedAdd(d);            }            // Right shift f k times until odd, left shift d k times            int trailingZeros = f.getLowestSetBit();            f.rightShift(trailingZeros);            d.leftShift(trailingZeros);            k += trailingZeros;        }        while (c.sign < 0)           c.signedAdd(p);        return fixup(c, p, k);    }    /*     * The Fixup Algorithm     * Calculates X such that X = C * 2^(-k) (mod P)     * Assumes C<P and P is odd.     */    static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,                                                                      int k) {        MutableBigInteger temp = new MutableBigInteger();        // Set r to the multiplicative inverse of p mod 2^32        int r = -inverseMod32(p.value[p.offset+p.intLen-1]);        for(int i=0, numWords = k >> 5; i<numWords; i++) {            // V = R * c (mod 2^j)            int  v = r * c.value[c.offset + c.intLen-1];            // c = c + (v * p)            p.mul(v, temp);            c.add(temp);            // c = c / 2^j            c.intLen--;        }        int numBits = k & 0x1f;        if (numBits != 0) {            // V = R * c (mod 2^j)            int v = r * c.value[c.offset + c.intLen-1];                     v &= ((1<<numBits) - 1);                     // c = c + (v * p)            p.mul(v, temp);            c.add(temp);              // c = c / 2^j            c.rightShift(numBits);        }        // In theory, c may be greater than p at this point (Very rare!)        while (c.compare(p) >= 0)            c.subtract(p);        return c;    }    /**     * Uses the extended Euclidean algorithm to compute the modInverse of base     * mod a modulus that is a power of 2. The modulus is 2^k.     */    MutableBigInteger euclidModInverse(int k) {        MutableBigInteger b = new MutableBigInteger(1);        b.leftShift(k);        MutableBigInteger mod = new MutableBigInteger(b);        MutableBigInteger a = new MutableBigInteger(this);        MutableBigInteger q = new MutableBigInteger();        MutableBigInteger r = new MutableBigInteger();               b.divide(a, q, r);        MutableBigInteger swapper = b; b = r; r = swapper;        MutableBigInteger t1 = new MutableBigInteger(q);        MutableBigInteger t0 = new MutableBigInteger(1);        MutableBigInteger temp = new MutableBigInteger();        while (!b.isOne()) {            a.divide(b, q, r);            if (r.intLen == 0)                throw new ArithmeticException("BigInteger not invertible.");                        swapper = r; r = a; a = swapper;            if (q.intLen == 1)                t1.mul(q.value[q.offset], temp);            else                q.multiply(t1, temp);            swapper = q; q = temp; temp = swapper;            t0.add(q);                       if (a.isOne())                return t0;                       b.divide(a, q, r);            if (r.intLen == 0)                throw new ArithmeticException("BigInteger not invertible.");            swapper = b; b = r; r = swapper;            if (q.intLen == 1)                t0.mul(q.value[q.offset], temp);            else                q.multiply(t0, temp);            swapper = q; q = temp; temp = swapper;            t1.add(q);        }        mod.subtract(t1);        return mod;    }}