biginteger.java

java源代码 请看看啊 提点宝贵的意见

         *
         * is a copy of everything below the main diagonal:
         *                       de
         *                 cd ce
         *           bc bd be
         *     ab ac ad ae
         *
         * Thus, the sum is 2 * (off the diagonal) + diagonal.
         *
         * This is accumulated beginning with the diagonal (which
         * consist of the squares of the digits of the input), which is then
         * divided by two, the off-diagonal added, and multiplied by two
         * again.  The low bit is simply a copy of the low bit of the
         * input, so it doesn't need special care.
         */
        int zlen = len << 1;
        if (z == null || z.length < zlen)
            z = new int[zlen];
        
        // Store the squares, right shifted one bit (i.e., divided by 2)
        int lastProductLowWord = 0;
        for (int j=0, i=0; j<len; j++) {
            long piece = (x[j] & LONG_MASK);
            long product = piece * piece;
            z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
            z[i++] = (int)(product >>> 1);
            lastProductLowWord = (int)product;
        }

        // Add in off-diagonal sums
        for (int i=len, offset=1; i>0; i--, offset+=2) {
            int t = x[i-1];
            t = mulAdd(z, x, offset, i-1, t);
            addOne(z, offset-1, i, t);
        }

        // Shift back up and set low bit
        primitiveLeftShift(z, zlen, 1);
        z[zlen-1] |= x[len-1] & 1;

        return z;
    }

    /**
     * Returns a BigInteger whose value is <tt>(this / val)</tt>.
     *
     * @param  val value by which this BigInteger is to be divided.
     * @return <tt>this / val</tt>
     * @throws ArithmeticException <tt>val==0</tt>
     */
    public BigInteger divide(BigInteger val) {
        MutableBigInteger q = new MutableBigInteger(),
                          r = new MutableBigInteger(),
                          a = new MutableBigInteger(this.mag),
                          b = new MutableBigInteger(val.mag);

        a.divide(b, q, r);
        return new BigInteger(q, this.signum * val.signum);
    }

    /**
     * Returns an array of two BigIntegers containing <tt>(this / val)</tt>
     * followed by <tt>(this % val)</tt>.
     *
     * @param  val value by which this BigInteger is to be divided, and the
     *	       remainder computed.
     * @return an array of two BigIntegers: the quotient <tt>(this / val)</tt>
     *	       is the initial element, and the remainder <tt>(this % val)</tt>
     *	       is the final element.
     * @throws ArithmeticException <tt>val==0</tt>
     */
    public BigInteger[] divideAndRemainder(BigInteger val) {
        BigInteger[] result = new BigInteger[2];
        MutableBigInteger q = new MutableBigInteger(),
                          r = new MutableBigInteger(),
                          a = new MutableBigInteger(this.mag),
                          b = new MutableBigInteger(val.mag);
        a.divide(b, q, r);
        result[0] = new BigInteger(q, this.signum * val.signum);
        result[1] = new BigInteger(r, this.signum);
        return result;
    }

    /**
     * Returns a BigInteger whose value is <tt>(this % val)</tt>.
     *
     * @param  val value by which this BigInteger is to be divided, and the
     *	       remainder computed.
     * @return <tt>this % val</tt>
     * @throws ArithmeticException <tt>val==0</tt>
     */
    public BigInteger remainder(BigInteger val) {
        MutableBigInteger q = new MutableBigInteger(),
                          r = new MutableBigInteger(),
                          a = new MutableBigInteger(this.mag),
                          b = new MutableBigInteger(val.mag);

        a.divide(b, q, r);
        return new BigInteger(r, this.signum);
    }

    /**
     * Returns a BigInteger whose value is <tt>(this<sup>exponent</sup>)</tt>.
     * Note that <tt>exponent</tt> is an integer rather than a BigInteger.
     *
     * @param  exponent exponent to which this BigInteger is to be raised.
     * @return <tt>this<sup>exponent</sup></tt>
     * @throws ArithmeticException <tt>exponent</tt> is negative.  (This would
     *	       cause the operation to yield a non-integer value.)
     */
    public BigInteger pow(int exponent) {
	if (exponent < 0)
	    throw new ArithmeticException("Negative exponent");
	if (signum==0)
	    return (exponent==0 ? ONE : this);

	// Perform exponentiation using repeated squaring trick
        int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
	int[] baseToPow2 = this.mag;
        int[] result = {1};

	while (exponent != 0) {
	    if ((exponent & 1)==1) {
		result = multiplyToLen(result, result.length, 
                                       baseToPow2, baseToPow2.length, null);
		result = trustedStripLeadingZeroInts(result);
	    }
	    if ((exponent >>>= 1) != 0) {
                baseToPow2 = squareToLen(baseToPow2, baseToPow2.length, null);
		baseToPow2 = trustedStripLeadingZeroInts(baseToPow2);
	    }
	}
	return new BigInteger(result, newSign);
    }

    /**
     * Returns a BigInteger whose value is the greatest common divisor of
     * <tt>abs(this)</tt> and <tt>abs(val)</tt>.  Returns 0 if
     * <tt>this==0 &amp;&amp; val==0</tt>.
     *
     * @param  val value with with the GCD is to be computed.
     * @return <tt>GCD(abs(this), abs(val))</tt>
     */
    public BigInteger gcd(BigInteger val) {
        if (val.signum == 0)
	    return this.abs();
	else if (this.signum == 0)
	    return val.abs();

        MutableBigInteger a = new MutableBigInteger(this);
        MutableBigInteger b = new MutableBigInteger(val);

        MutableBigInteger result = a.hybridGCD(b);

        return new BigInteger(result, 1);
    }

    /**
     * Left shift int array a up to len by n bits. Returns the array that
     * results from the shift since space may have to be reallocated.
     */
    private static int[] leftShift(int[] a, int len, int n) {
        int nInts = n >>> 5;
        int nBits = n&0x1F;
        int bitsInHighWord = bitLen(a[0]);
        
        // If shift can be done without recopy, do so
        if (n <= (32-bitsInHighWord)) {
            primitiveLeftShift(a, len, nBits);
            return a;
        } else { // Array must be resized
            if (nBits <= (32-bitsInHighWord)) {
                int result[] = new int[nInts+len];
                for (int i=0; i<len; i++)
                    result[i] = a[i];
                primitiveLeftShift(result, result.length, nBits);
                return result;
            } else {
                int result[] = new int[nInts+len+1];
                for (int i=0; i<len; i++)
                    result[i] = a[i];
                primitiveRightShift(result, result.length, 32 - nBits);
                return result;
            }
        }
    }

    // shifts a up to len right n bits assumes no leading zeros, 0<n<32
    static void primitiveRightShift(int[] a, int len, int n) {
        int n2 = 32 - n;
        for (int i=len-1, c=a[i]; i>0; i--) {
            int b = c;
            c = a[i-1];
            a[i] = (c << n2) | (b >>> n);
        }
        a[0] >>>= n;
    }

    // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
    static void primitiveLeftShift(int[] a, int len, int n) {
        if (len == 0 || n == 0)
            return;

        int n2 = 32 - n;
        for (int i=0, c=a[i], m=i+len-1; i<m; i++) {
            int b = c;
            c = a[i+1];
            a[i] = (b << n) | (c >>> n2);
        }
        a[len-1] <<= n;
    }

    /**
     * Calculate bitlength of contents of the first len elements an int array,
     * assuming there are no leading zero ints.
     */
    private static int bitLength(int[] val, int len) {
        if (len==0)
            return 0;
        return ((len-1)<<5) + bitLen(val[0]);
    }

    /**
     * Returns a BigInteger whose value is the absolute value of this
     * BigInteger. 
     *
     * @return <tt>abs(this)</tt>
     */
    public BigInteger abs() {
	return (signum >= 0 ? this : this.negate());
    }

    /**
     * Returns a BigInteger whose value is <tt>(-this)</tt>.
     *
     * @return <tt>-this</tt>
     */
    public BigInteger negate() {
	return new BigInteger(this.mag, -this.signum);
    }

    /**
     * Returns the signum function of this BigInteger.
     *
     * @return -1, 0 or 1 as the value of this BigInteger is negative, zero or
     *	       positive.
     */
    public int signum() {
	return this.signum;
    }

    // Modular Arithmetic Operations

    /**
     * Returns a BigInteger whose value is <tt>(this mod m</tt>).  This method
     * differs from <tt>remainder</tt> in that it always returns a
     * <i>non-negative</i> BigInteger.
     *
     * @param  m the modulus.
     * @return <tt>this mod m</tt>
     * @throws ArithmeticException <tt>m &lt;= 0</tt>
     * @see    #remainder
     */
    public BigInteger mod(BigInteger m) {
	if (m.signum <= 0)
	    throw new ArithmeticException("BigInteger: modulus not positive");

	BigInteger result = this.remainder(m);
	return (result.signum >= 0 ? result : result.add(m));
    }

    /**
     * Returns a BigInteger whose value is
     * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike <tt>pow</tt>, this
     * method permits negative exponents.)
     *
     * @param  exponent the exponent.
     * @param  m the modulus.
     * @return <tt>this<sup>exponent</sup> mod m</tt>
     * @throws ArithmeticException <tt>m &lt;= 0</tt>
     * @see    #modInverse
     */
    public BigInteger modPow(BigInteger exponent, BigInteger m) {
	if (m.signum <= 0)
	    throw new ArithmeticException("BigInteger: modulus not positive");

	// Trivial cases
	if (exponent.signum == 0)
            return (m.equals(ONE) ? ZERO : ONE);

        if (this.equals(ONE))
            return (m.equals(ONE) ? ZERO : ONE);

        if (this.equals(ZERO) && exponent.signum >= 0)
            return ZERO;

        if (this.equals(negConst[1]) && (!exponent.testBit(0)))
            return (m.equals(ONE) ? ZERO : ONE);
            
	boolean invertResult;
	if ((invertResult = (exponent.signum < 0)))
	    exponent = exponent.negate();

	BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
			   ? this.mod(m) : this);
	BigInteger result;
	if (m.testBit(0)) { // odd modulus
	    result = base.oddModPow(exponent, m);
	} else {
	    /*
	     * Even modulus.  Tear it into an "odd part" (m1) and power of two
             * (m2), exponentiate mod m1, manually exponentiate mod m2, and
             * use Chinese Remainder Theorem to combine results.
	     */

	    // Tear m apart into odd part (m1) and power of 2 (m2)
	    int p = m.getLowestSetBit();   // Max pow of 2 that divides m

	    BigInteger m1 = m.shiftRight(p);  // m/2**p
	    BigInteger m2 = ONE.shiftLeft(p); // 2**p

            // Calculate new base from m1
            BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
                                ? this.mod(m1) : this);

            // Caculate (base ** exponent) mod m1.
            BigInteger a1 = (m1.equals(ONE) ? ZERO :
                             base2.oddModPow(exponent, m1));

	    // Calculate (this ** exponent) mod m2
	    BigInteger a2 = base.modPow2(exponent, p);

	    // Combine results using Chinese Remainder Theorem
	    BigInteger y1 = m2.modInverse(m1);
	    BigInteger y2 = m1.modInverse(m2);

	    result = a1.multiply(m2).multiply(y1).add
		     (a2.multiply(m1).multiply(y2)).mod(m);
	}

	return (invertResult ? result.modInverse(m) : result);
    }

    static int[] bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
                                                Integer.MAX_VALUE}; // Sentinel

    /**
     * Returns a BigInteger whose value is x to the power of y mod z.
     * Assumes: z is odd && x < z.
     */
    private BigInteger oddModPow(BigInteger y, BigInteger z) {
    /*
     * The algorithm is adapted from Colin Plumb's C library.
     *
     * The window algorithm:
     * The idea is to keep a running product of b1 = n^(high-order bits of exp)
     * and then keep appending exponent bits to it.  The following patterns
     * apply to a 3-bit window (k = 3):
     * To append   0: square
     * To append   1: square, multiply by n^1
     * To append  10: square, multiply by n^1, square
     * To append  11: square, square, multiply by n^3
     * To append 100: square, multiply by n^1, square, square
     * To append 101: square, square, square, multiply by n^5
     * To append 110: square, square, multiply by n^3, square
     * To append 111: square, square, square, multiply by n^7
     *
     * Since each pattern involves only one multiply, the longer the pattern
     * the better, except that a 0 (no multiplies) can be appended directly.
     * We precompute a table of odd powers of n, up to 2^k, and can then
     * multiply k bits of exponent at a time.  Actually, assuming random
     * exponents, there is on average one zero bit between needs to
     * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
     * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
     * you have to do one multiply per k+1 bits of exponent.
     *
     * The loop walks down the exponent, squaring the result buffer as
     * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
     * filled with the upcoming exponent bits.  (What is read after the
     * end of the exponent is unimportant, but it is filled with zero here.)
     * When the most-significant bit of this buffer becomes set, i.e.
     * (buf & tblmask) != 0, we have to decide what pattern to multiply
     * by, and when to do it.  We decide, remember to do it in future