Java > Open Source Codes > java > math > BigInteger


1 /*
2  * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
3  * SUN PROPRIETARY/CONFIDENTIAL. Use is subject to license terms.
4  */
5
6 /*
7  * @(#)BigInteger.java 1.70 05/08/09
8  */
9
10 package java.math;
11
12 import java.util.Random ;
13 import java.io.*;
14
15 /**
16  * Immutable arbitrary-precision integers. All operations behave as if
17  * BigIntegers were represented in two's-complement notation (like Java's
18  * primitive integer types). BigInteger provides analogues to all of Java's
19  * primitive integer operators, and all relevant methods from java.lang.Math.
20  * Additionally, BigInteger provides operations for modular arithmetic, GCD
21  * calculation, primality testing, prime generation, bit manipulation,
22  * and a few other miscellaneous operations.
23  * <p>
24  * Semantics of arithmetic operations exactly mimic those of Java's integer
25  * arithmetic operators, as defined in <i>The Java Language Specification</i>.
26  * For example, division by zero throws an <tt>ArithmeticException</tt>, and
27  * division of a negative by a positive yields a negative (or zero) remainder.
28  * All of the details in the Spec concerning overflow are ignored, as
29  * BigIntegers are made as large as necessary to accommodate the results of an
30  * operation.
31  * <p>
32  * Semantics of shift operations extend those of Java's shift operators
33  * to allow for negative shift distances. A right-shift with a negative
34  * shift distance results in a left shift, and vice-versa. The unsigned
35  * right shift operator (&gt;&gt;&gt;) is omitted, as this operation makes
36  * little sense in combination with the "infinite word size" abstraction
37  * provided by this class.
38  * <p>
39  * Semantics of bitwise logical operations exactly mimic those of Java's
40  * bitwise integer operators. The binary operators (<tt>and</tt>,
41  * <tt>or</tt>, <tt>xor</tt>) implicitly perform sign extension on the shorter
42  * of the two operands prior to performing the operation.
43  * <p>
44  * Comparison operations perform signed integer comparisons, analogous to
45  * those performed by Java's relational and equality operators.
46  * <p>
47  * Modular arithmetic operations are provided to compute residues, perform
48  * exponentiation, and compute multiplicative inverses. These methods always
49  * return a non-negative result, between <tt>0</tt> and <tt>(modulus - 1)</tt>,
50  * inclusive.
51  * <p>
52  * Bit operations operate on a single bit of the two's-complement
53  * representation of their operand. If necessary, the operand is sign-
54  * extended so that it contains the designated bit. None of the single-bit
55  * operations can produce a BigInteger with a different sign from the
56  * BigInteger being operated on, as they affect only a single bit, and the
57  * "infinite word size" abstraction provided by this class ensures that there
58  * are infinitely many "virtual sign bits" preceding each BigInteger.
59  * <p>
60  * For the sake of brevity and clarity, pseudo-code is used throughout the
61  * descriptions of BigInteger methods. The pseudo-code expression
62  * <tt>(i + j)</tt> is shorthand for "a BigInteger whose value is
63  * that of the BigInteger <tt>i</tt> plus that of the BigInteger <tt>j</tt>."
64  * The pseudo-code expression <tt>(i == j)</tt> is shorthand for
65  * "<tt>true</tt> if and only if the BigInteger <tt>i</tt> represents the same
66  * value as the BigInteger <tt>j</tt>." Other pseudo-code expressions are
67  * interpreted similarly.
68  * <p>
69  * All methods and constructors in this class throw
70  * <CODE>NullPointerException</CODE> when passed
71  * a null object reference for any input parameter.
72  *
73  * @see BigDecimal
74  * @version 1.70, 08/09/05
75  * @author Josh Bloch
76  * @author Michael McCloskey
77  * @since JDK1.1
78  */
79
80 public class BigInteger extends Number implements Comparable <BigInteger > {
81     /**
82      * The signum of this BigInteger: -1 for negative, 0 for zero, or
83      * 1 for positive. Note that the BigInteger zero <i>must</i> have
84      * a signum of 0. This is necessary to ensures that there is exactly one
85      * representation for each BigInteger value.
86      *
87      * @serial
88      */
89     int signum;
90
91     /**
92      * The magnitude of this BigInteger, in <i>big-endian</i> order: the
93      * zeroth element of this array is the most-significant int of the
94      * magnitude. The magnitude must be "minimal" in that the most-significant
95      * int (<tt>mag[0]</tt>) must be non-zero. This is necessary to
96      * ensure that there is exactly one representation for each BigInteger
97      * value. Note that this implies that the BigInteger zero has a
98      * zero-length mag array.
99      */
100     int[] mag;
101
102     // These "redundant fields" are initialized with recognizable nonsense
103 // values, and cached the first time they are needed (or never, if they
104 // aren't needed).
105
106     /**
107      * The bitCount of this BigInteger, as returned by bitCount(), or -1
108      * (either value is acceptable).
109      *
110      * @serial
111      * @see #bitCount
112      */
113     private int bitCount = -1;
114
115     /**
116      * The bitLength of this BigInteger, as returned by bitLength(), or -1
117      * (either value is acceptable).
118      *
119      * @serial
120      * @see #bitLength
121      */
122     private int bitLength = -1;
123
124     /**
125      * The lowest set bit of this BigInteger, as returned by getLowestSetBit(),
126      * or -2 (either value is acceptable).
127      *
128      * @serial
129      * @see #getLowestSetBit
130      */
131     private int lowestSetBit = -2;
132
133     /**
134      * The index of the lowest-order byte in the magnitude of this BigInteger
135      * that contains a nonzero byte, or -2 (either value is acceptable). The
136      * least significant byte has int-number 0, the next byte in order of
137      * increasing significance has byte-number 1, and so forth.
138      *
139      * @serial
140      */
141     private int firstNonzeroByteNum = -2;
142
143     /**
144      * The index of the lowest-order int in the magnitude of this BigInteger
145      * that contains a nonzero int, or -2 (either value is acceptable). The
146      * least significant int has int-number 0, the next int in order of
147      * increasing significance has int-number 1, and so forth.
148      */
149     private int firstNonzeroIntNum = -2;
150
151     /**
152      * This mask is used to obtain the value of an int as if it were unsigned.
153      */
154     private final static long LONG_MASK = 0xffffffffL;
155
156     //Constructors
157
158     /**
159      * Translates a byte array containing the two's-complement binary
160      * representation of a BigInteger into a BigInteger. The input array is
161      * assumed to be in <i>big-endian</i> byte-order: the most significant
162      * byte is in the zeroth element.
163      *
164      * @param val big-endian two's-complement binary representation of
165      * BigInteger.
166      * @throws NumberFormatException <tt>val</tt> is zero bytes long.
167      */
168     public BigInteger(byte[] val) {
169     if (val.length == 0)
170         throw new NumberFormatException ("Zero length BigInteger");
171
172     if (val[0] < 0) {
173             mag = makePositive(val);
174         signum = -1;
175     } else {
176         mag = stripLeadingZeroBytes(val);
177         signum = (mag.length == 0 ? 0 : 1);
178     }
179     }
180
181     /**
182      * This private constructor translates an int array containing the
183      * two's-complement binary representation of a BigInteger into a
184      * BigInteger. The input array is assumed to be in <i>big-endian</i>
185      * int-order: the most significant int is in the zeroth element.
186      */
187     private BigInteger(int[] val) {
188     if (val.length == 0)
189         throw new NumberFormatException ("Zero length BigInteger");
190
191     if (val[0] < 0) {
192             mag = makePositive(val);
193         signum = -1;
194     } else {
195         mag = trustedStripLeadingZeroInts(val);
196         signum = (mag.length == 0 ? 0 : 1);
197     }
198     }
199
200     /**
201      * Translates the sign-magnitude representation of a BigInteger into a
202      * BigInteger. The sign is represented as an integer signum value: -1 for
203      * negative, 0 for zero, or 1 for positive. The magnitude is a byte array
204      * in <i>big-endian</i> byte-order: the most significant byte is in the
205      * zeroth element. A zero-length magnitude array is permissible, and will
206      * result inin a BigInteger value of 0, whether signum is -1, 0 or 1.
207      *
208      * @param signum signum of the number (-1 for negative, 0 for zero, 1
209      * for positive).
210      * @param magnitude big-endian binary representation of the magnitude of
211      * the number.
212      * @throws NumberFormatException <tt>signum</tt> is not one of the three
213      * legal values (-1, 0, and 1), or <tt>signum</tt> is 0 and
214      * <tt>magnitude</tt> contains one or more non-zero bytes.
215      */
216     public BigInteger(int signum, byte[] magnitude) {
217     this.mag = stripLeadingZeroBytes(magnitude);
218
219     if (signum < -1 || signum > 1)
220         throw(new NumberFormatException ("Invalid signum value"));
221
222     if (this.mag.length==0) {
223         this.signum = 0;
224     } else {
225         if (signum == 0)
226         throw(new NumberFormatException ("signum-magnitude mismatch"));
227         this.signum = signum;
228     }
229     }
230
231     /**
232      * A constructor for internal use that translates the sign-magnitude
233      * representation of a BigInteger into a BigInteger. It checks the
234      * arguments and copies the magnitude so this constructor would be
235      * safe for external use.
236      */
237     private BigInteger(int signum, int[] magnitude) {
238     this.mag = stripLeadingZeroInts(magnitude);
239
240     if (signum < -1 || signum > 1)
241         throw(new NumberFormatException ("Invalid signum value"));
242
243     if (this.mag.length==0) {
244         this.signum = 0;
245     } else {
246         if (signum == 0)
247         throw(new NumberFormatException ("signum-magnitude mismatch"));
248         this.signum = signum;
249     }
250     }
251
252     /**
253      * Translates the String representation of a BigInteger in the specified
254      * radix into a BigInteger. The String representation consists of an
255      * optional minus sign followed by a sequence of one or more digits in the
256      * specified radix. The character-to-digit mapping is provided by
257      * <tt>Character.digit</tt>. The String may not contain any extraneous
258      * characters (whitespace, for example).
259      *
260      * @param val String representation of BigInteger.
261      * @param radix radix to be used in interpreting <tt>val</tt>.
262      * @throws NumberFormatException <tt>val</tt> is not a valid representation
263      * of a BigInteger in the specified radix, or <tt>radix</tt> is
264      * outside the range from {@link Character#MIN_RADIX} to
265      * {@link Character#MAX_RADIX}, inclusive.
266      * @see Character#digit
267      */
268     public BigInteger(String val, int radix) {
269     int cursor = 0, numDigits;
270         int len = val.length();
271
272     if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
273         throw new NumberFormatException ("Radix out of range");
274     if (val.length() == 0)
275         throw new NumberFormatException ("Zero length BigInteger");
276
277     // Check for minus sign
278 signum = 1;
279         int index = val.lastIndexOf("-");
280         if (index != -1) {
281             if (index == 0) {
282                 if (val.length() == 1)
283                     throw new NumberFormatException ("Zero length BigInteger");
284                 signum = -1;
285                 cursor = 1;
286             } else {
287                 throw new NumberFormatException ("Illegal embedded minus sign");
288             }
289         }
290
291         // Skip leading zeros and compute number of digits in magnitude
292 while (cursor < len &&
293                Character.digit(val.charAt(cursor),radix) == 0)
294         cursor++;
295     if (cursor == len) {
296         signum = 0;
297         mag = ZERO.mag;
298         return;
299     } else {
300         numDigits = len - cursor;
301     }
302
303         // Pre-allocate array of expected size. May be too large but can
304 // never be too small. Typically exact.
305 int numBits = (int)(((numDigits * bitsPerDigit[radix]) >>> 10) + 1);
306         int numWords = (numBits + 31) /32;
307         mag = new int[numWords];
308
309     // Process first (potentially short) digit group
310 int firstGroupLen = numDigits % digitsPerInt[radix];
311     if (firstGroupLen == 0)
312         firstGroupLen = digitsPerInt[radix];
313     String group = val.substring(cursor, cursor += firstGroupLen);
314         mag[mag.length - 1] = Integer.parseInt(group, radix);
315     if (mag[mag.length - 1] < 0)
316         throw new NumberFormatException ("Illegal digit");
317         
318     // Process remaining digit groups
319 int superRadix = intRadix[radix];
320         int groupVal = 0;
321     while (cursor < val.length()) {
322         group = val.substring(cursor, cursor += digitsPerInt[radix]);
323         groupVal = Integer.parseInt(group, radix);
324         if (groupVal < 0)
325         throw new NumberFormatException ("Illegal digit");
326             destructiveMulAdd(mag, superRadix, groupVal);
327     }
328         // Required for cases where the array was overallocated.
329 mag = trustedStripLeadingZeroInts(mag);
330     }
331
332     // Constructs a new BigInteger using a char array with radix=10
333 BigInteger(char[] val) {
334         int cursor = 0, numDigits;
335         int len = val.length;
336
337     // Check for leading minus sign
338 signum = 1;
339     if (val[0] == '-') {
340         if (len == 1)
341         throw new NumberFormatException ("Zero length BigInteger");
342         signum = -1;
343         cursor = 1;
344     }
345
346         // Skip leading zeros and compute number of digits in magnitude
347 while (cursor < len && Character.digit(val[cursor], 10) == 0)
348         cursor++;
349     if (cursor == len) {
350         signum = 0;
351         mag = ZERO.mag;
352         return;
353     } else {
354         numDigits = len - cursor;
355     }
356
357         // Pre-allocate array of expected size
358 int numWords;
359         if (len < 10) {
360             numWords = 1;
361         } else {
362             int numBits = (int)(((numDigits * bitsPerDigit[10]) >>> 10) + 1);
363             numWords = (numBits + 31) /32;
364         }
365         mag = new int[numWords];
366  
367     // Process first (potentially short) digit group
368 int firstGroupLen = numDigits % digitsPerInt[10];
369     if (firstGroupLen == 0)
370         firstGroupLen = digitsPerInt[10];
371         mag[mag.length-1] = parseInt(val, cursor, cursor += firstGroupLen);
372         
373     // Process remaining digit groups
374 while (cursor < len) {
375         int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
376             destructiveMulAdd(mag, intRadix[10], groupVal);
377     }
378         mag = trustedStripLeadingZeroInts(mag);
379     }
380
381     // Create an integer with the digits between the two indexes
382 // Assumes start < end. The result may be negative, but it
383 // is to be treated as an unsigned value.
384 private int parseInt(char[] source, int start, int end) {
385         int result = Character.digit(source[start++], 10);
386         if (result == -1)
387             throw new NumberFormatException (new String (source));
388
389         for (int index = start; index<end; index++) {
390             int nextVal = Character.digit(source[index], 10);
391             if (nextVal == -1)
392                 throw new NumberFormatException (new String (source));
393             result = 10*result + nextVal;
394         }
395
396         return result;
397     }
398
399     // bitsPerDigit in the given radix times 1024
400 // Rounded up to avoid underallocation.
401 private static long bitsPerDigit[] = { 0, 0,
402         1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
403         3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
404         4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
405                                            5253, 5295};
406
407     // Multiply x array times word y in place, and add word z
408 private static void destructiveMulAdd(int[] x, int y, int z) {
409         // Perform the multiplication word by word
410 long ylong = y & LONG_MASK;
411         long zlong = z & LONG_MASK;
412         int len = x.length;
413
414         long product = 0;
415         long carry = 0;
416         for (int i = len-1; i >= 0; i--) {
417             product = ylong * (x[i] & LONG_MASK) + carry;
418             x[i] = (int)product;
419             carry = product >>> 32;
420         }
421
422         // Perform the addition
423 long sum = (x[len-1] & LONG_MASK) + zlong;
424         x[len-1] = (int)sum;
425         carry = sum >>> 32;
426         for (int i = len-2; i >= 0; i--) {
427             sum = (x[i] & LONG_MASK) + carry;
428             x[i] = (int)sum;
429             carry = sum >>> 32;
430         }
431     }
432
433     /**
434      * Translates the decimal String representation of a BigInteger into a
435      * BigInteger. The String representation consists of an optional minus
436      * sign followed by a sequence of one or more decimal digits. The
437      * character-to-digit mapping is provided by <tt>Character.digit</tt>.
438      * The String may not contain any extraneous characters (whitespace, for
439      * example).
440      *
441      * @param val decimal String representation of BigInteger.
442      * @throws NumberFormatException <tt>val</tt> is not a valid representation
443      * of a BigInteger.
444      * @see Character#digit
445      */
446     public BigInteger(String val) {
447     this(val, 10);
448     }
449
450     /**
451      * Constructs a randomly generated BigInteger, uniformly distributed over
452      * the range <tt>0</tt> to <tt>(2<sup>numBits</sup> - 1)</tt>, inclusive.
453      * The uniformity of the distribution assumes that a fair source of random
454      * bits is provided in <tt>rnd</tt>. Note that this constructor always
455      * constructs a non-negative BigInteger.
456      *
457      * @param numBits maximum bitLength of the new BigInteger.
458      * @param rnd source of randomness to be used in computing the new
459      * BigInteger.
460      * @throws IllegalArgumentException <tt>numBits</tt> is negative.
461      * @see #bitLength
462      */
463     public BigInteger(int numBits, Random rnd) {
464     this(1, randomBits(numBits, rnd));
465     }
466
467     private static byte[] randomBits(int numBits, Random rnd) {
468     if (numBits < 0)
469         throw new IllegalArgumentException ("numBits must be non-negative");
470     int numBytes = (numBits+7)/8;
471     byte[] randomBits = new byte[numBytes];
472
473     // Generate random bytes and mask out any excess bits
474 if (numBytes > 0) {
475         rnd.nextBytes(randomBits);
476         int excessBits = 8*numBytes - numBits;
477         randomBits[0] &= (1 << (8-excessBits)) - 1;
478     }
479     return randomBits;
480     }
481
482     /**
483      * Constructs a randomly generated positive BigInteger that is probably
484      * prime, with the specified bitLength.<p>
485      *
486      * It is recommended that the {@link #probablePrime probablePrime}
487      * method be used in preference to this constructor unless there
488      * is a compelling need to specify a certainty.
489      *
490      * @param bitLength bitLength of the returned BigInteger.
491      * @param certainty a measure of the uncertainty that the caller is
492      * willing to tolerate. The probability that the new BigInteger
493      * represents a prime number will exceed
494      * <tt>(1 - 1/2<sup>certainty</sup></tt>). The execution time of
495      * this constructor is proportional to the value of this parameter.
496      * @param rnd source of random bits used to select candidates to be
497      * tested for primality.
498      * @throws ArithmeticException <tt>bitLength &lt; 2</tt>.
499      * @see #bitLength
500      */
501     public BigInteger(int bitLength, int certainty, Random rnd) {
502         BigInteger prime;
503
504     if (bitLength < 2)
505         throw new ArithmeticException ("bitLength < 2");
506         // The cutoff of 95 was chosen empirically for best performance
507 prime = (bitLength < 95 ? smallPrime(bitLength, certainty, rnd)
508                                 : largePrime(bitLength, certainty, rnd));
509     signum = 1;
510     mag = prime.mag;
511     }
512
513     // Minimum size in bits that the requested prime number has
514 // before we use the large prime number generating algorithms
515 private static final int SMALL_PRIME_THRESHOLD = 95;
516
517     // Certainty required to meet the spec of probablePrime
518 private static final int DEFAULT_PRIME_CERTAINTY = 100;
519
520     /**
521      * Returns a positive BigInteger that is probably prime, with the
522      * specified bitLength. The probability that a BigInteger returned
523      * by this method is composite does not exceed 2<sup>-100</sup>.
524      *
525      * @param bitLength bitLength of the returned BigInteger.
526      * @param rnd source of random bits used to select candidates to be
527      * tested for primality.
528      * @return a BigInteger of <tt>bitLength</tt> bits that is probably prime
529      * @throws ArithmeticException <tt>bitLength &lt; 2</tt>.
530      * @see #bitLength
531      */
532     public static BigInteger probablePrime(int bitLength, Random rnd) {
533     if (bitLength < 2)
534         throw new ArithmeticException ("bitLength < 2");
535
536         // The cutoff of 95 was chosen empirically for best performance
537 return (bitLength < SMALL_PRIME_THRESHOLD ?
538                 smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
539                 largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
540     }
541
542     /**
543      * Find a random number of the specified bitLength that is probably prime.
544      * This method is used for smaller primes, its performance degrades on
545      * larger bitlengths.
546      *
547      * This method assumes bitLength > 1.
548      */
549     private static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
550         int magLen = (bitLength + 31) >>> 5;
551         int temp[] = new int[magLen];
552         int highBit = 1 << ((bitLength+31) & 0x1f); // High bit of high int
553 int highMask = (highBit << 1) - 1; // Bits to keep in high int
554
555         while(true) {
556             // Construct a candidate
557 for (int i=0; i<magLen; i++)
558                 temp[i] = rnd.nextInt();
559             temp[0] = (temp[0] & highMask) | highBit; // Ensure exact length
560 if (bitLength > 2)
561                 temp[magLen-1] |= 1; // Make odd if bitlen > 2
562
563             BigInteger p = new BigInteger (temp, 1);
564
565             // Do cheap "pre-test" if applicable
566 if (bitLength > 6) {
567                 long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
568                 if ((r%3==0) || (r%5==0) || (r%7==0) || (r%11==0) ||
569                     (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
570                     (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
571                     continue; // Candidate is composite; try another
572 }
573             
574             // All candidates of bitLength 2 and 3 are prime by this point
575 if (bitLength < 4)
576                 return p;
577
578             // Do expensive test if we survive pre-test (or it's inapplicable)
579 if (p.primeToCertainty(certainty))
580                 return p;
581         }
582     }
583
584     private static final BigInteger SMALL_PRIME_PRODUCT
585                        = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);
586
587     /**
588      * Find a random number of the specified bitLength that is probably prime.
589      * This method is more appropriate for larger bitlengths since it uses
590      * a sieve to eliminate most composites before using a more expensive
591      * test.
592      */
593     private static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
594         BigInteger p;
595         p = new BigInteger (bitLength, rnd).setBit(bitLength-1);
596         p.mag[p.mag.length-1] &= 0xfffffffe;
597
598         // Use a sieve length likely to contain the next prime number
599 int searchLen = (bitLength / 20) * 64;
600         BitSieve searchSieve = new BitSieve (p, searchLen);
601         BigInteger candidate = searchSieve.retrieve(p, certainty);
602
603         while ((candidate == null) || (candidate.bitLength() != bitLength)) {
604             p = p.add(BigInteger.valueOf(2*searchLen));
605             if (p.bitLength() != bitLength)
606                 p = new BigInteger (bitLength, rnd).setBit(bitLength-1);
607             p.mag[p.mag.length-1] &= 0xfffffffe;
608             searchSieve = new BitSieve (p, searchLen);
609             candidate = searchSieve.retrieve(p, certainty);
610         }
611         return candidate;
612     }
613
614    /**
615     * Returns the first integer greater than this <code>BigInteger</code> that
616     * is probably prime. The probability that the number returned by this
617     * method is composite does not exceed 2<sup>-100</sup>. This method will
618     * never skip over a prime when searching: if it returns <tt>p</tt>, there
619     * is no prime <tt>q</tt> such that <tt>this &lt; q &lt; p</tt>.
620     *
621     * @return the first integer greater than this <code>BigInteger</code> that
622     * is probably prime.
623     * @throws ArithmeticException <tt>this &lt; 0</tt>.
624     * @since 1.5
625     */
626     public BigInteger nextProbablePrime() {
627         if (this.signum < 0)
628             throw new ArithmeticException ("start < 0: " + this);
629         
630         // Handle trivial cases
631 if ((this.signum == 0) || this.equals(ONE))
632             return TWO;
633
634         BigInteger result = this.add(ONE);
635
636         // Fastpath for small numbers
637 if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
638  
639             // Ensure an odd number
640 if (!result.testBit(0))
641                 result = result.add(ONE);
642
643             while(true) {
644                 // Do cheap "pre-test" if applicable
645 if (result.bitLength() > 6) {
646                     long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
647                     if ((r%3==0) || (r%5==0) || (r%7==0) || (r%11==0) ||
648                         (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
649                         (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
650                         result = result.add(TWO);
651                         continue; // Candidate is composite; try another
652 }
653                 }
654             
655                 // All candidates of bitLength 2 and 3 are prime by this point
656 if (result.bitLength() < 4)
657                     return result;
658
659                 // The expensive test
660 if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY))
661                     return result;
662
663                 result = result.add(TWO);
664             }
665         }
666
667         // Start at previous even number
668 if (result.testBit(0))
669             result = result.subtract(ONE);
670
671         // Looking for the next large prime
672 int searchLen = (result.bitLength() / 20) * 64;
673
674         while(true) {
675            BitSieve searchSieve = new BitSieve (result, searchLen);
676            BigInteger candidate = searchSieve.retrieve(result,
677                                                      DEFAULT_PRIME_CERTAINTY);
678            if (candidate != null)
679                return candidate;
680            result = result.add(BigInteger.valueOf(2 * searchLen));
681         }
682     }
683
684     /**
685      * Returns <tt>true</tt> if this BigInteger is probably prime,
686      * <tt>false</tt> if it's definitely composite.
687      *
688      * This method assumes bitLength > 2.
689      *
690      * @param certainty a measure of the uncertainty that the caller is
691      * willing to tolerate: if the call returns <tt>true</tt>
692      * the probability that this BigInteger is prime exceeds
693      * <tt>(1 - 1/2<sup>certainty</sup>)</tt>. The execution time of
694      * this method is proportional to the value of this parameter.
695      * @return <tt>true</tt> if this BigInteger is probably prime,
696      * <tt>false</tt> if it's definitely composite.
697      */
698     boolean primeToCertainty(int certainty) {
699         int rounds = 0;
700         int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;
701
702         // The relationship between the certainty and the number of rounds
703 // we perform is given in the draft standard ANSI X9.80, "PRIME
704 // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
705 int sizeInBits = this.bitLength();
706         if (sizeInBits < 100) {
707             rounds = 50;
708             rounds = n < rounds ? n : rounds;
709             return passesMillerRabin(rounds);
710         }
711
712         if (sizeInBits < 256) {
713             rounds = 27;
714         } else if (sizeInBits < 512) {
715             rounds = 15;
716         } else if (sizeInBits < 768) {
717             rounds = 8;
718         } else if (sizeInBits < 1024) {
719             rounds = 4;
720         } else {
721             rounds = 2;
722         }
723         rounds = n < rounds ? n : rounds;
724
725         return passesMillerRabin(rounds) && passesLucasLehmer();
726     }
727
728     /**
729      * Returns true iff this BigInteger is a Lucas-Lehmer probable prime.
730      *
731      * The following assumptions are made:
732      * This BigInteger is a positive, odd number.
733      */
734     private boolean passesLucasLehmer() {
735         BigInteger thisPlusOne = this.add(ONE);
736
737         // Step 1
738 int d = 5;
739         while (jacobiSymbol(d, this) != -1) {
740             // 5, -7, 9, -11, ...
741 d = (d<0) ? Math.abs(d)+2 : -(d+2);
742         }
743         
744         // Step 2
745 BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);
746
747         // Step 3
748 return u.mod(this).equals(ZERO);
749     }
750
751     /**
752      * Computes Jacobi(p,n).
753      * Assumes n positive, odd, n>=3.
754      */
755     private static int jacobiSymbol(int p, BigInteger n) {
756         if (p == 0)
757             return 0;
758
759         // Algorithm and comments adapted from Colin Plumb's C library.
760 int j = 1;
761     int u = n.mag[n.mag.length-1];
762
763         // Make p positive
764 if (p < 0) {
765             p = -p;
766             int n8 = u & 7;
767             if ((n8 == 3) || (n8 == 7))
768                 j = -j; // 3 (011) or 7 (111) mod 8
769 }
770
771     // Get rid of factors of 2 in p
772 while ((p & 3) == 0)
773             p >>= 2;
774     if ((p & 1) == 0) {
775             p >>= 1;
776             if (((u ^ (u>>1)) & 2) != 0)
777                 j = -j; // 3 (011) or 5 (101) mod 8
778 }
779     if (p == 1)
780         return j;
781     // Then, apply quadratic reciprocity
782 if ((p & u & 2) != 0) // p = u = 3 (mod 4)?
783 j = -j;
784     // And reduce u mod p
785 u = n.mod(BigInteger.valueOf(p)).intValue();
786
787     // Now compute Jacobi(u,p), u < p
788 while (u != 0) {
789             while ((u & 3) == 0)
790                 u >>= 2;
791             if ((u & 1) == 0) {
792                 u >>= 1;
793                 if (((p ^ (p>>1)) & 2) != 0)
794                     j = -j; // 3 (011) or 5 (101) mod 8
795 }
796             if (u == 1)
797                 return j;
798             // Now both u and p are odd, so use quadratic reciprocity
799 assert (u < p);
800             int t = u; u = p; p = t;
801             if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
802 j = -j;
803             // Now u >= p, so it can be reduced
804 u %= p;
805     }
806     return 0;
807     }
808
809     private static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
810         BigInteger d = BigInteger.valueOf(z);
811         BigInteger u = ONE; BigInteger u2;
812         BigInteger v = ONE; BigInteger v2;
813
814         for (int i=k.bitLength()-2; i>=0; i--) {
815             u2 = u.multiply(v).mod(n);
816
817             v2 = v.square().add(d.multiply(u.square())).mod(n);
818             if (v2.testBit(0)) {
819                 v2 = n.subtract(v2);
820                 v2.signum = - v2.signum;
821             }
822             v2 = v2.shiftRight(1);
823
824             u = u2; v = v2;
825             if (k.testBit(i)) {
826                 u2 = u.add(v).mod(n);
827                 if (u2.testBit(0)) {
828                     u2 = n.subtract(u2);
829                     u2.signum = - u2.signum;
830                 }
831                 u2 = u2.shiftRight(1);
832                 
833                 v2 = v.add(d.multiply(u)).mod(n);
834                 if (v2.testBit(0)) {
835                     v2 = n.subtract(v2);
836                     v2.signum = - v2.signum;
837                 }
838                 v2 = v2.shiftRight(1);
839
840                 u = u2; v = v2;
841             }
842         }
843         return u;
844     }
845
846     /**
847      * Returns true iff this BigInteger passes the specified number of
848      * Miller-Rabin tests. This test is taken from the DSA spec (NIST FIPS
849      * 186-2).
850  &nb