biginteger.java

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

    private int parseInt(char[] source, int start, int end) {
        int result = Character.digit(source[start++], 10);
        if (result == -1)
            throw new NumberFormatException(new String(source));

        for (int index = start; index<end; index++) {
            int nextVal = Character.digit(source[index], 10);
            if (nextVal == -1)
                throw new NumberFormatException(new String(source));
            result = 10*result + nextVal;
        }

        return result;
    }

    // bitsPerDigit in the given radix times 1024
    // Rounded up to avoid underallocation.
    private static long bitsPerDigit[] = { 0, 0,
        1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
        3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
        4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
                                           5253, 5295};

    // Multiply x array times word y in place, and add word z
    private static void destructiveMulAdd(int[] x, int y, int z) {
        // Perform the multiplication word by word
        long ylong = y & LONG_MASK;
        long zlong = z & LONG_MASK;
        int len = x.length;

        long product = 0;
        long carry = 0;
        for (int i = len-1; i >= 0; i--) {
            product = ylong * (x[i] & LONG_MASK) + carry;
            x[i] = (int)product;
            carry = product >>> 32;
        }

        // Perform the addition
        long sum = (x[len-1] & LONG_MASK) + zlong;
        x[len-1] = (int)sum;
        carry = sum >>> 32;
        for (int i = len-2; i >= 0; i--) {
            sum = (x[i] & LONG_MASK) + carry;
            x[i] = (int)sum;
            carry = sum >>> 32;
        }
    }

    /**
     * Translates the decimal String representation of a BigInteger into a
     * BigInteger.  The String representation consists of an optional minus
     * sign followed by a sequence of one or more decimal digits.  The
     * character-to-digit mapping is provided by <tt>Character.digit</tt>.
     * The String may not contain any extraneous characters (whitespace, for
     * example).
     *
     * @param val decimal String representation of BigInteger.
     * @throws NumberFormatException <tt>val</tt> is not a valid representation
     *	       of a BigInteger.
     * @see    Character#digit
     */
    public BigInteger(String val) {
	this(val, 10);
    }

    /**
     * Constructs a randomly generated BigInteger, uniformly distributed over
     * the range <tt>0</tt> to <tt>(2<sup>numBits</sup> - 1)</tt>, inclusive.
     * The uniformity of the distribution assumes that a fair source of random
     * bits is provided in <tt>rnd</tt>.  Note that this constructor always
     * constructs a non-negative BigInteger.
     *
     * @param  numBits maximum bitLength of the new BigInteger.
     * @param  rnd source of randomness to be used in computing the new
     *	       BigInteger.
     * @throws IllegalArgumentException <tt>numBits</tt> is negative.
     * @see #bitLength
     */
    public BigInteger(int numBits, Random rnd) {
	this(1, randomBits(numBits, rnd));
    }

    private static byte[] randomBits(int numBits, Random rnd) {
	if (numBits < 0)
	    throw new IllegalArgumentException("numBits must be non-negative");
	int numBytes = (numBits+7)/8;
	byte[] randomBits = new byte[numBytes];

	// Generate random bytes and mask out any excess bits
	if (numBytes > 0) {
	    rnd.nextBytes(randomBits);
	    int excessBits = 8*numBytes - numBits;
	    randomBits[0] &= (1 << (8-excessBits)) - 1;
	}
	return randomBits;
    }

    /**
     * Constructs a randomly generated positive BigInteger that is probably
     * prime, with the specified bitLength.<p>
     *
     * It is recommended that the {@link #probablePrime probablePrime}
     * method be used in preference to this constructor unless there
     * is a compelling need to specify a certainty.
     *
     * @param  bitLength bitLength of the returned BigInteger.
     * @param  certainty a measure of the uncertainty that the caller is
     *         willing to tolerate.  The probability that the new BigInteger
     *	       represents a prime number will exceed
     *	       <tt>(1 - 1/2<sup>certainty</sup></tt>).  The execution time of
     *	       this constructor is proportional to the value of this parameter.
     * @param  rnd source of random bits used to select candidates to be
     *	       tested for primality.
     * @throws ArithmeticException <tt>bitLength &lt; 2</tt>.
     * @see    #bitLength
     */
    public BigInteger(int bitLength, int certainty, Random rnd) {
        BigInteger prime;

	if (bitLength < 2)
	    throw new ArithmeticException("bitLength < 2");
        // The cutoff of 95 was chosen empirically for best performance
        prime = (bitLength < 95 ? smallPrime(bitLength, certainty, rnd)
                                : largePrime(bitLength, certainty, rnd));
	signum = 1;
	mag = prime.mag;
    }

    // Minimum size in bits that the requested prime number has
    // before we use the large prime number generating algorithms
    private static final int SMALL_PRIME_THRESHOLD = 95;

    /**
     * Returns a positive BigInteger that is probably prime, with the
     * specified bitLength. The probability that a BigInteger returned
     * by this method is composite does not exceed 2<sup>-100</sup>.
     *
     * @param  bitLength bitLength of the returned BigInteger.
     * @param  rnd source of random bits used to select candidates to be
     *	       tested for primality.
     * @return a BigInteger of <tt>bitLength</tt> bits that is probably prime
     * @throws ArithmeticException <tt>bitLength &lt; 2</tt>.
     * @see    #bitLength
     */
    public static BigInteger probablePrime(int bitLength, Random rnd) {
	if (bitLength < 2)
	    throw new ArithmeticException("bitLength < 2");

        // The cutoff of 95 was chosen empirically for best performance
        return (bitLength < SMALL_PRIME_THRESHOLD ?
                smallPrime(bitLength, 100, rnd) :
                largePrime(bitLength, 100, rnd));
    }

    /**
     * Find a random number of the specified bitLength that is probably prime.
     * This method is used for smaller primes, its performance degrades on
     * larger bitlengths.
     *
     * This method assumes bitLength > 1.
     */
    private static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
        int magLen = (bitLength + 31) >>> 5;
        int temp[] = new int[magLen];
        int highBit = 1 << ((bitLength+31) & 0x1f);  // High bit of high int
        int highMask = (highBit << 1) - 1;  // Bits to keep in high int

        while(true) {
            // Construct a candidate
            for (int i=0; i<magLen; i++)
                temp[i] = rnd.nextInt();
            temp[0] = (temp[0] & highMask) | highBit;  // Ensure exact length
            if (bitLength > 2)
                temp[magLen-1] |= 1;  // Make odd if bitlen > 2

            BigInteger p = new BigInteger(temp, 1);

            // Do cheap "pre-test" if applicable
            if (bitLength > 6) {
                long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
                if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) || 
                    (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) || 
                    (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
                    continue; // Candidate is composite; try another
            }
            
            // All candidates of bitLength 2 and 3 are prime by this point
            if (bitLength < 4)
                return p;

            // Do expensive test if we survive pre-test (or it's inapplicable)
            if (p.primeToCertainty(certainty))
                return p;
        }
    }

    private static final BigInteger SMALL_PRIME_PRODUCT
                       = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);

    /**
     * Find a random number of the specified bitLength that is probably prime.
     * This method is more appropriate for larger bitlengths since it uses
     * a sieve to eliminate most composites before using a more expensive
     * test.
     */
    private static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
        BigInteger p;
        p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
        p.mag[p.mag.length-1] &= 0xfffffffe;

        // Use a sieve length likely to contain the next prime number
        int searchLen = (bitLength / 20) * 64;
        BitSieve searchSieve = new BitSieve(p, searchLen);
        BigInteger candidate = searchSieve.retrieve(p, certainty);

        while ((candidate == null) || (candidate.bitLength() != bitLength)) {
            p = p.add(BigInteger.valueOf(2*searchLen));
            if (p.bitLength() != bitLength)
                p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
            p.mag[p.mag.length-1] &= 0xfffffffe;
            searchSieve = new BitSieve(p, searchLen);
            candidate = searchSieve.retrieve(p, certainty);
        }
        return candidate;
    }

    /**
     * Returns <tt>true</tt> if this BigInteger is probably prime,
     * <tt>false</tt> if it's definitely composite.
     *
     * This method assumes bitLength > 2.
     *
     * @param  certainty a measure of the uncertainty that the caller is
     *	       willing to tolerate: if the call returns <tt>true</tt>
     *	       the probability that this BigInteger is prime exceeds
     *	       <tt>(1 - 1/2<sup>certainty</sup>)</tt>.  The execution time of
     * 	       this method is proportional to the value of this parameter.
     * @return <tt>true</tt> if this BigInteger is probably prime,
     * 	       <tt>false</tt> if it's definitely composite.
     */
    boolean primeToCertainty(int certainty) {
        int rounds = 0;
        int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;

        // The relationship between the certainty and the number of rounds
        // we perform is given in the draft standard ANSI X9.80, "PRIME
        // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
        int sizeInBits = this.bitLength();
        if (sizeInBits < 100) {
            rounds = 50;
            rounds = n < rounds ? n : rounds;
            return passesMillerRabin(rounds);
        }

        if (sizeInBits < 256) {
            rounds = 27;
        } else if (sizeInBits < 512) {
            rounds = 15;
        } else if (sizeInBits < 768) {
            rounds = 8;
        } else if (sizeInBits < 1024) {
            rounds = 4;
        } else {
            rounds = 2;
        }
        rounds = n < rounds ? n : rounds;

        return passesMillerRabin(rounds) && passesLucasLehmer();
    }

    /**
     * Returns true iff this BigInteger is a Lucas-Lehmer probable prime.
     *
     * The following assumptions are made:
     * This BigInteger is a positive, odd number.
     */
    private boolean passesLucasLehmer() {
        BigInteger thisPlusOne = this.add(ONE);

        // Step 1
        int d = 5;
        while (jacobiSymbol(d, this) != -1) {
            // 5, -7, 9, -11, ...
            d = (d<0) ? Math.abs(d)+2 : -(d+2);
        }
        
        // Step 2
        BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);

        // Step 3
        return u.mod(this).equals(ZERO);
    }

    /**
     * Computes Jacobi(p,n).
     * Assumes n is positive, odd.
     */
    int jacobiSymbol(int p, BigInteger n) {
        if (p == 0)
            return 0;

        // Algorithm and comments adapted from Colin Plumb's C library.
	int j = 1;
	int u = n.mag[n.mag.length-1];

        // Make p positive
        if (p < 0) {
            p = -p;
            int n8 = u & 7;
            if ((n8 == 3) || (n8 == 7))
                j = -j; // 3 (011) or 7 (111) mod 8
        }

	// Get rid of factors of 2 in p
	while ((p & 3) == 0)
            p >>= 2;
	if ((p & 1) == 0) {
            p >>= 1;
            if (((u ^ u>>1) & 2) != 0)
                j = -j;	// 3 (011) or 5 (101) mod 8
	}
	if (p == 1)
	    return j;
	// Then, apply quadratic reciprocity
	if ((p & u & 2) != 0)	// p = u = 3 (mod 4)?
	    j = -j;
	// And reduce u mod p
	u = n.mod(BigInteger.valueOf(p)).intValue();

	// Now compute Jacobi(u,p), u < p
	while (u != 0) {
            while ((u & 3) == 0)
                u >>= 2;
            if ((u & 1) == 0) {
                u >>= 1;
                if (((p ^ p>>1) & 2) != 0)
                    j = -j;	// 3 (011) or 5 (101) mod 8
            }
            if (u == 1)
                return j;
            // Now both u and p are odd, so use quadratic reciprocity
            if (u < p) {
                int t = u; u = p; p = t;
                if ((u & p & 2)	!= 0)// u = p = 3 (mod 4)?
                    j = -j;
            }
            // Now u >= p, so it can be reduced
            u %= p;
	}
	return 0;
    }

    private static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
        BigInteger d = BigInteger.valueOf(z);
        BigInteger u = ONE; BigInteger u2;
        BigInteger v = ONE; BigInteger v2;

        for (int i=k.bitLength()-2; i>=0; i--) {
            u2 = u.multiply(v).mod(n);

            v2 = v.square().add(d.multiply(u.square())).mod(n);
            if (v2.testBit(0)) {
                v2 = n.subtract(v2);
                v2.signum = - v2.signum;
            }
            v2 = v2.shiftRight(1);

            u = u2; v = v2;
            if (k.testBit(i)) {
                u2 = u.add(v).mod(n);
                if (u2.testBit(0)) {
                    u2 = n.subtract(u2);
                    u2.signum = - u2.signum;
                }
                u2 = u2.shiftRight(1);
                
                v2 = v.add(d.multiply(u)).mod(n);
                if (v2.testBit(0)) {
                    v2 = n.subtract(v2);
                    v2.signum = - v2.signum;
                }
                v2 = v2.shiftRight(1);