OSDN Git Service

2001-08-17 Mark J Roberts <mjr@anarcast.net>
[pf3gnuchains/gcc-fork.git] / libjava / java / math / BigInteger.java
1 // BigInteger.java -- an arbitrary-precision integer
2
3 /* Copyright (C) 1999, 2000, 2001  Free Software Foundation
4
5    This file is part of libgcj.
6
7 This software is copyrighted work licensed under the terms of the
8 Libgcj License.  Please consult the file "LIBGCJ_LICENSE" for
9 details.  */
10
11 package java.math;
12 import gnu.gcj.math.*;
13 import java.util.Random;
14 import java.io.ObjectInputStream;
15 import java.io.ObjectOutputStream;
16 import java.io.IOException;
17
18 /**
19  * @author Warren Levy <warrenl@cygnus.com>
20  * @date December 20, 1999.
21  */
22
23 /**
24  * Written using on-line Java Platform 1.2 API Specification, as well
25  * as "The Java Class Libraries", 2nd edition (Addison-Wesley, 1998) and
26  * "Applied Cryptography, Second Edition" by Bruce Schneier (Wiley, 1996).
27
28  * 
29  * Based primarily on IntNum.java BitOps.java by Per Bothner <per@bothner.com>
30  * (found in Kawa 1.6.62).
31  *
32  * Status:  Believed complete and correct.
33  */
34
35 public class BigInteger extends Number implements Comparable
36 {
37   /** All integers are stored in 2's-complement form.
38    * If words == null, the ival is the value of this BigInteger.
39    * Otherwise, the first ival elements of words make the value
40    * of this BigInteger, stored in little-endian order, 2's-complement form. */
41   transient private int ival;
42   transient private int[] words;
43
44   // Serialization fields.
45   private int bitCount = -1;
46   private int bitLength = -1;
47   private int firstNonzeroByteNum = -2;
48   private int lowestSetBit = -2;
49   private byte[] magnitude;
50   private int signum;
51   private static final long serialVersionUID = -8287574255936472291L;
52
53
54   /** We pre-allocate integers in the range minFixNum..maxFixNum. */
55   private static final int minFixNum = -100;
56   private static final int maxFixNum = 1024;
57   private static final int numFixNum = maxFixNum-minFixNum+1;
58   private static final BigInteger[] smallFixNums = new BigInteger[numFixNum];
59
60   static {
61     for (int i = numFixNum;  --i >= 0; )
62       smallFixNums[i] = new BigInteger(i + minFixNum);
63   }
64
65   // JDK1.2
66   public static final BigInteger ZERO = smallFixNums[-minFixNum];
67
68   // JDK1.2
69   public static final BigInteger ONE = smallFixNums[1 - minFixNum];
70
71   /* Rounding modes: */
72   private static final int FLOOR = 1;
73   private static final int CEILING = 2;
74   private static final int TRUNCATE = 3;
75   private static final int ROUND = 4;
76
77   /** When checking the probability of primes, it is most efficient to
78    * first check the factoring of small primes, so we'll use this array.
79    */
80   private static final int[] primes =
81     {   2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,
82        47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101, 103, 107,
83       109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
84       191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251 };
85
86   private BigInteger()
87   {
88   }
89
90   /* Create a new (non-shared) BigInteger, and initialize to an int. */
91   private BigInteger(int value)
92   {
93     ival = value;
94   }
95
96   public BigInteger(String val, int radix)
97   {
98     BigInteger result = valueOf(val, radix);
99     this.ival = result.ival;
100     this.words = result.words;
101   }
102
103   public BigInteger(String val)
104   {
105     this(val, 10);
106   }
107
108   /* Create a new (non-shared) BigInteger, and initialize from a byte array. */
109   public BigInteger(byte[] val)
110   {
111     if (val == null || val.length < 1)
112       throw new NumberFormatException();
113
114     words = byteArrayToIntArray(val, val[0] < 0 ? -1 : 0);
115     BigInteger result = make(words, words.length);
116     this.ival = result.ival;
117     this.words = result.words;
118   }
119
120   public BigInteger(int signum, byte[] magnitude)
121   {
122     if (magnitude == null || signum > 1 || signum < -1)
123       throw new NumberFormatException();
124
125     if (signum == 0)
126       {
127         int i;
128         for (i = magnitude.length - 1; i >= 0 && magnitude[i] == 0; --i)
129           ;
130         if (i >= 0)
131           throw new NumberFormatException();
132         return;
133       }
134
135     // Magnitude is always positive, so don't ever pass a sign of -1.
136     words = byteArrayToIntArray(magnitude, 0);
137     BigInteger result = make(words, words.length);
138     this.ival = result.ival;
139     this.words = result.words;
140
141     if (signum < 0)
142       setNegative();
143   }
144
145   public BigInteger(int numBits, Random rnd)
146   {
147     this(1, randBytes(numBits, rnd));
148   }
149
150   private static byte[] randBytes(int numBits, Random rnd)
151   {
152     if (numBits < 0)
153       throw new IllegalArgumentException();
154
155     int extra = numBits % 8;
156     byte[] b = new byte[numBits / 8 + (extra > 0 ? 1 : 0)];
157     rnd.nextBytes(b);
158     if (extra > 0)
159       b[0] &= ~((~0) << (8 - extra));
160     return b;
161   }
162
163   public BigInteger(int bitLength, int certainty, Random rnd)
164   {
165     this(bitLength, rnd);
166
167     // Keep going until we find a probable prime.
168     while (true)
169       {
170         if (isProbablePrime(certainty))
171           return;
172
173         BigInteger next = new BigInteger(bitLength, rnd);
174         this.ival = next.ival;
175         this.words = next.words;
176       }
177   }
178
179   /** Return a (possibly-shared) BigInteger with a given long value. */
180   private static BigInteger make(long value)
181   {
182     if (value >= minFixNum && value <= maxFixNum)
183       return smallFixNums[(int)value - minFixNum];
184     int i = (int) value;
185     if ((long)i == value)
186       return new BigInteger(i);
187     BigInteger result = alloc(2);
188     result.ival = 2;
189     result.words[0] = i;
190     result.words[1] = (int) (value >> 32);
191     return result;
192   }
193
194   // FIXME: Could simply rename 'make' method above as valueOf while
195   // changing all instances of 'make'.  Don't do this until this class
196   // is done as the Kawa class this is based on has 'make' methods
197   // with other parameters; wait to see if they are used in BigInteger.
198   public static BigInteger valueOf(long val)
199   {
200     return make(val);
201   }
202
203   /** Make a canonicalized BigInteger from an array of words.
204    * The array may be reused (without copying). */
205   private static BigInteger make(int[] words, int len)
206   {
207     if (words == null)
208       return make(len);
209     len = BigInteger.wordsNeeded(words, len);
210     if (len <= 1)
211       return len == 0 ? ZERO : make(words[0]);
212     BigInteger num = new BigInteger();
213     num.words = words;
214     num.ival = len;
215     return num;
216   }
217
218   /** Convert a big-endian byte array to a little-endian array of words. */
219   private static int[] byteArrayToIntArray(byte[] bytes, int sign)
220   {
221     // Determine number of words needed.
222     int[] words = new int[bytes.length/4 + 1];
223     int nwords = words.length;
224
225     // Create a int out of modulo 4 high order bytes.
226     int bptr = 0;
227     int word = sign;
228     for (int i = bytes.length % 4; i > 0; --i, bptr++)
229       word = (word << 8) | (((int) bytes[bptr]) & 0xff);
230     words[--nwords] = word;
231
232     // Elements remaining in byte[] are a multiple of 4.
233     while (nwords > 0)
234       words[--nwords] = bytes[bptr++] << 24 |
235                         (((int) bytes[bptr++]) & 0xff) << 16 |
236                         (((int) bytes[bptr++]) & 0xff) << 8 |
237                         (((int) bytes[bptr++]) & 0xff);
238     return words;
239   }
240
241   /** Allocate a new non-shared BigInteger.
242    * @param nwords number of words to allocate
243    */
244   private static BigInteger alloc(int nwords)
245   {
246     if (nwords <= 1)
247       return new BigInteger();
248     BigInteger result = new BigInteger();
249     result.words = new int[nwords];
250     return result;
251   }
252
253   /** Change words.length to nwords.
254    * We allow words.length to be upto nwords+2 without reallocating.
255    */
256   private void realloc(int nwords)
257   {
258     if (nwords == 0)
259       {
260         if (words != null)
261           {
262             if (ival > 0)
263               ival = words[0];
264             words = null;
265           }
266       }
267     else if (words == null
268              || words.length < nwords
269              || words.length > nwords + 2)
270       {
271         int[] new_words = new int [nwords];
272         if (words == null)
273           {
274             new_words[0] = ival;
275             ival = 1;
276           }
277         else
278           {
279             if (nwords < ival)
280               ival = nwords;
281             System.arraycopy(words, 0, new_words, 0, ival);
282           }
283         words = new_words;
284       }
285   }
286
287   private final boolean isNegative()
288   {
289     return (words == null ? ival : words[ival - 1]) < 0;
290   }
291
292   public int signum()
293   {
294     int top = words == null ? ival : words[ival-1];
295     if (top == 0 && words == null)
296       return 0;
297     return top < 0 ? -1 : 1;
298   }
299
300   private static int compareTo(BigInteger x, BigInteger y)
301   {
302     if (x.words == null && y.words == null)
303       return x.ival < y.ival ? -1 : x.ival > y.ival ? 1 : 0;
304     boolean x_negative = x.isNegative();
305     boolean y_negative = y.isNegative();
306     if (x_negative != y_negative)
307       return x_negative ? -1 : 1;
308     int x_len = x.words == null ? 1 : x.ival;
309     int y_len = y.words == null ? 1 : y.ival;
310     if (x_len != y_len)
311       return (x_len > y_len) != x_negative ? 1 : -1;
312     return MPN.cmp(x.words, y.words, x_len);
313   }
314
315   // JDK1.2
316   public int compareTo(Object obj)
317   {
318     if (obj instanceof BigInteger)
319       return compareTo(this, (BigInteger) obj);
320     throw new ClassCastException();
321   }
322
323   public int compareTo(BigInteger val)
324   {
325     return compareTo(this, val);
326   }
327
328   public BigInteger min(BigInteger val)
329   {
330     return compareTo(this, val) < 0 ? this : val;
331   }
332
333   public BigInteger max(BigInteger val)
334   {
335     return compareTo(this, val) > 0 ? this : val;
336   }
337
338   private final boolean isOdd()
339   {
340     int low = words == null ? ival : words[0];
341     return (low & 1) != 0;
342   }
343
344   private final boolean isZero()
345   {
346     return words == null && ival == 0;
347   }
348
349   private final boolean isOne()
350   {
351     return words == null && ival == 1;
352   }
353
354   private final boolean isMinusOne()
355   {
356     return words == null && ival == -1;
357   }
358
359   /** Calculate how many words are significant in words[0:len-1].
360    * Returns the least value x such that x>0 && words[0:x-1]==words[0:len-1],
361    * when words is viewed as a 2's complement integer.
362    */
363   private static int wordsNeeded(int[] words, int len)
364   {
365     int i = len;
366     if (i > 0)
367       {
368         int word = words[--i];
369         if (word == -1)
370           {
371             while (i > 0 && (word = words[i - 1]) < 0)
372               {
373                 i--;
374                 if (word != -1) break;
375               }
376           }
377         else
378           {
379             while (word == 0 && i > 0 && (word = words[i - 1]) >= 0)  i--;
380           }
381       }
382     return i + 1;
383   }
384
385   private BigInteger canonicalize()
386   {
387     if (words != null
388         && (ival = BigInteger.wordsNeeded(words, ival)) <= 1)
389       {
390         if (ival == 1)
391           ival = words[0];
392         words = null;
393       }
394     if (words == null && ival >= minFixNum && ival <= maxFixNum)
395       return smallFixNums[(int) ival - minFixNum];
396     return this;
397   }
398
399   /** Add two ints, yielding a BigInteger. */
400   private static final BigInteger add(int x, int y)
401   {
402     return BigInteger.make((long) x + (long) y);
403   }
404
405   /** Add a BigInteger and an int, yielding a new BigInteger. */
406   private static BigInteger add(BigInteger x, int y)
407   {
408     if (x.words == null)
409       return BigInteger.add(x.ival, y);
410     BigInteger result = new BigInteger(0);
411     result.setAdd(x, y);
412     return result.canonicalize();
413   }
414
415   /** Set this to the sum of x and y.
416    * OK if x==this. */
417   private void setAdd(BigInteger x, int y)
418   {
419     if (x.words == null)
420       {
421         set((long) x.ival + (long) y);
422         return;
423       }
424     int len = x.ival;
425     realloc(len + 1);
426     long carry = y;
427     for (int i = 0;  i < len;  i++)
428       {
429         carry += ((long) x.words[i] & 0xffffffffL);
430         words[i] = (int) carry;
431         carry >>= 32;
432       }
433     if (x.words[len - 1] < 0)
434       carry--;
435     words[len] = (int) carry;
436     ival = wordsNeeded(words, len + 1);
437   }
438
439   /** Destructively add an int to this. */
440   private final void setAdd(int y)
441   {
442     setAdd(this, y);
443   }
444
445   /** Destructively set the value of this to a long. */
446   private final void set(long y)
447   {
448     int i = (int) y;
449     if ((long) i == y)
450       {
451         ival = i;
452         words = null;
453       }
454     else
455       {
456         realloc(2);
457         words[0] = i;
458         words[1] = (int) (y >> 32);
459         ival = 2;
460       }
461   }
462
463   /** Destructively set the value of this to the given words.
464   * The words array is reused, not copied. */
465   private final void set(int[] words, int length)
466   {
467     this.ival = length;
468     this.words = words;
469   }
470
471   /** Destructively set the value of this to that of y. */
472   private final void set(BigInteger y)
473   {
474     if (y.words == null)
475       set(y.ival);
476     else if (this != y)
477       {
478         realloc(y.ival);
479         System.arraycopy(y.words, 0, words, 0, y.ival);
480         ival = y.ival;
481       }
482   }
483
484   /** Add two BigIntegers, yielding their sum as another BigInteger. */
485   private static BigInteger add(BigInteger x, BigInteger y, int k)
486   {
487     if (x.words == null && y.words == null)
488       return BigInteger.make((long) k * (long) y.ival + (long) x.ival);
489     if (k != 1)
490       {
491         if (k == -1)
492           y = BigInteger.neg(y);
493         else
494           y = BigInteger.times(y, BigInteger.make(k));
495       }
496     if (x.words == null)
497       return BigInteger.add(y, x.ival);
498     if (y.words == null)
499       return BigInteger.add(x, y.ival);
500     // Both are big
501     int len;
502     if (y.ival > x.ival)
503       { // Swap so x is longer then y.
504         BigInteger tmp = x;  x = y;  y = tmp;
505       }
506     BigInteger result = alloc(x.ival + 1);
507     int i = y.ival;
508     long carry = MPN.add_n(result.words, x.words, y.words, i);
509     long y_ext = y.words[i - 1] < 0 ? 0xffffffffL : 0;
510     for (; i < x.ival;  i++)
511       {
512         carry += ((long) x.words[i] & 0xffffffffL) + y_ext;;
513         result.words[i] = (int) carry;
514         carry >>>= 32;
515       }
516     if (x.words[i - 1] < 0)
517       y_ext--;
518     result.words[i] = (int) (carry + y_ext);
519     result.ival = i+1;
520     return result.canonicalize();
521   }
522
523   public BigInteger add(BigInteger val)
524   {
525     return add(this, val, 1);
526   }
527
528   public BigInteger subtract(BigInteger val)
529   {
530     return add(this, val, -1);
531   }
532
533   private static final BigInteger times(BigInteger x, int y)
534   {
535     if (y == 0)
536       return ZERO;
537     if (y == 1)
538       return x;
539     int[] xwords = x.words;
540     int xlen = x.ival;
541     if (xwords == null)
542       return BigInteger.make((long) xlen * (long) y);
543     boolean negative;
544     BigInteger result = BigInteger.alloc(xlen + 1);
545     if (xwords[xlen - 1] < 0)
546       {
547         negative = true;
548         negate(result.words, xwords, xlen);
549         xwords = result.words;
550       }
551     else
552       negative = false;
553     if (y < 0)
554       {
555         negative = !negative;
556         y = -y;
557       }
558     result.words[xlen] = MPN.mul_1(result.words, xwords, xlen, y);
559     result.ival = xlen + 1;
560     if (negative)
561       result.setNegative();
562     return result.canonicalize();
563   }
564
565   private static final BigInteger times(BigInteger x, BigInteger y)
566   {
567     if (y.words == null)
568       return times(x, y.ival);
569     if (x.words == null)
570       return times(y, x.ival);
571     boolean negative = false;
572     int[] xwords;
573     int[] ywords;
574     int xlen = x.ival;
575     int ylen = y.ival;
576     if (x.isNegative())
577       {
578         negative = true;
579         xwords = new int[xlen];
580         negate(xwords, x.words, xlen);
581       }
582     else
583       {
584         negative = false;
585         xwords = x.words;
586       }
587     if (y.isNegative())
588       {
589         negative = !negative;
590         ywords = new int[ylen];
591         negate(ywords, y.words, ylen);
592       }
593     else
594       ywords = y.words;
595     // Swap if x is shorter then y.
596     if (xlen < ylen)
597       {
598         int[] twords = xwords;  xwords = ywords;  ywords = twords;
599         int tlen = xlen;  xlen = ylen;  ylen = tlen;
600       }
601     BigInteger result = BigInteger.alloc(xlen+ylen);
602     MPN.mul(result.words, xwords, xlen, ywords, ylen);
603     result.ival = xlen+ylen;
604     if (negative)
605       result.setNegative();
606     return result.canonicalize();
607   }
608
609   public BigInteger multiply(BigInteger y)
610   {
611     return times(this, y);
612   }
613
614   private static void divide(long x, long y,
615                              BigInteger quotient, BigInteger remainder,
616                              int rounding_mode)
617   {
618     boolean xNegative, yNegative;
619     if (x < 0)
620       {
621         xNegative = true;
622         if (x == Long.MIN_VALUE)
623           {
624             divide(BigInteger.make(x), BigInteger.make(y),
625                    quotient, remainder, rounding_mode);
626             return;
627           }
628         x = -x;
629       }
630     else
631       xNegative = false;
632
633     if (y < 0)
634       {
635         yNegative = true;
636         if (y == Long.MIN_VALUE)
637           {
638             if (rounding_mode == TRUNCATE)
639               { // x != Long.Min_VALUE implies abs(x) < abs(y)
640                 if (quotient != null)
641                   quotient.set(0);
642                 if (remainder != null)
643                   remainder.set(x);
644               }
645             else
646               divide(BigInteger.make(x), BigInteger.make(y),
647                       quotient, remainder, rounding_mode);
648             return;
649           }
650         y = -y;
651       }
652     else
653       yNegative = false;
654
655     long q = x / y;
656     long r = x % y;
657     boolean qNegative = xNegative ^ yNegative;
658
659     boolean add_one = false;
660     if (r != 0)
661       {
662         switch (rounding_mode)
663           {
664           case TRUNCATE:
665             break;
666           case CEILING:
667           case FLOOR:
668             if (qNegative == (rounding_mode == FLOOR))
669               add_one = true;
670             break;
671           case ROUND:
672             add_one = r > ((y - (q & 1)) >> 1);
673             break;
674           }
675       }
676     if (quotient != null)
677       {
678         if (add_one)
679           q++;
680         if (qNegative)
681           q = -q;
682         quotient.set(q);
683       }
684     if (remainder != null)
685       {
686         // The remainder is by definition: X-Q*Y
687         if (add_one)
688           {
689             // Subtract the remainder from Y.
690             r = y - r;
691             // In this case, abs(Q*Y) > abs(X).
692             // So sign(remainder) = -sign(X).
693             xNegative = ! xNegative;
694           }
695         else
696           {
697             // If !add_one, then: abs(Q*Y) <= abs(X).
698             // So sign(remainder) = sign(X).
699           }
700         if (xNegative)
701           r = -r;
702         remainder.set(r);
703       }
704   }
705
706   /** Divide two integers, yielding quotient and remainder.
707    * @param x the numerator in the division
708    * @param y the denominator in the division
709    * @param quotient is set to the quotient of the result (iff quotient!=null)
710    * @param remainder is set to the remainder of the result
711    *  (iff remainder!=null)
712    * @param rounding_mode one of FLOOR, CEILING, TRUNCATE, or ROUND.
713    */
714   private static void divide(BigInteger x, BigInteger y,
715                              BigInteger quotient, BigInteger remainder,
716                              int rounding_mode)
717   {
718     if ((x.words == null || x.ival <= 2)
719         && (y.words == null || y.ival <= 2))
720       {
721         long x_l = x.longValue();
722         long y_l = y.longValue();
723         if (x_l != Long.MIN_VALUE && y_l != Long.MIN_VALUE)
724           {
725             divide(x_l, y_l, quotient, remainder, rounding_mode);
726             return;
727           }
728       }
729
730     boolean xNegative = x.isNegative();
731     boolean yNegative = y.isNegative();
732     boolean qNegative = xNegative ^ yNegative;
733
734     int ylen = y.words == null ? 1 : y.ival;
735     int[] ywords = new int[ylen];
736     y.getAbsolute(ywords);
737     while (ylen > 1 && ywords[ylen - 1] == 0)  ylen--;
738
739     int xlen = x.words == null ? 1 : x.ival;
740     int[] xwords = new int[xlen+2];
741     x.getAbsolute(xwords);
742     while (xlen > 1 && xwords[xlen-1] == 0)  xlen--;
743
744     int qlen, rlen;
745
746     int cmpval = MPN.cmp(xwords, xlen, ywords, ylen);
747     if (cmpval < 0)  // abs(x) < abs(y)
748       { // quotient = 0;  remainder = num.
749         int[] rwords = xwords;  xwords = ywords;  ywords = rwords;
750         rlen = xlen;  qlen = 1;  xwords[0] = 0;
751       }
752     else if (cmpval == 0)  // abs(x) == abs(y)
753       {
754         xwords[0] = 1;  qlen = 1;  // quotient = 1
755         ywords[0] = 0;  rlen = 1;  // remainder = 0;
756       }
757     else if (ylen == 1)
758       {
759         qlen = xlen;
760         // Need to leave room for a word of leading zeros if dividing by 1
761         // and the dividend has the high bit set.  It might be safe to
762         // increment qlen in all cases, but it certainly is only necessary
763         // in the following case.
764         if (ywords[0] == 1 && xwords[xlen-1] < 0)
765           qlen++;
766         rlen = 1;
767         ywords[0] = MPN.divmod_1(xwords, xwords, xlen, ywords[0]);
768       }
769     else  // abs(x) > abs(y)
770       {
771         // Normalize the denominator, i.e. make its most significant bit set by
772         // shifting it normalization_steps bits to the left.  Also shift the
773         // numerator the same number of steps (to keep the quotient the same!).
774
775         int nshift = MPN.count_leading_zeros(ywords[ylen - 1]);
776         if (nshift != 0)
777           {
778             // Shift up the denominator setting the most significant bit of
779             // the most significant word.
780             MPN.lshift(ywords, 0, ywords, ylen, nshift);
781
782             // Shift up the numerator, possibly introducing a new most
783             // significant word.
784             int x_high = MPN.lshift(xwords, 0, xwords, xlen, nshift);
785             xwords[xlen++] = x_high;
786           }
787
788         if (xlen == ylen)
789           xwords[xlen++] = 0;
790         MPN.divide(xwords, xlen, ywords, ylen);
791         rlen = ylen;
792         MPN.rshift0 (ywords, xwords, 0, rlen, nshift);
793
794         qlen = xlen + 1 - ylen;
795         if (quotient != null)
796           {
797             for (int i = 0;  i < qlen;  i++)
798               xwords[i] = xwords[i+ylen];
799           }
800       }
801
802     if (ywords[rlen-1] < 0)
803       {
804         ywords[rlen] = 0;
805         rlen++;
806       }
807
808     // Now the quotient is in xwords, and the remainder is in ywords.
809
810     boolean add_one = false;
811     if (rlen > 1 || ywords[0] != 0)
812       { // Non-zero remainder i.e. in-exact quotient.
813         switch (rounding_mode)
814           {
815           case TRUNCATE:
816             break;
817           case CEILING:
818           case FLOOR:
819             if (qNegative == (rounding_mode == FLOOR))
820               add_one = true;
821             break;
822           case ROUND:
823             // int cmp = compareTo(remainder<<1, abs(y));
824             BigInteger tmp = remainder == null ? new BigInteger() : remainder;
825             tmp.set(ywords, rlen);
826             tmp = shift(tmp, 1);
827             if (yNegative)
828               tmp.setNegative();
829             int cmp = compareTo(tmp, y);
830             // Now cmp == compareTo(sign(y)*(remainder<<1), y)
831             if (yNegative)
832               cmp = -cmp;
833             add_one = (cmp == 1) || (cmp == 0 && (xwords[0]&1) != 0);
834           }
835       }
836     if (quotient != null)
837       {
838         quotient.set(xwords, qlen);
839         if (qNegative)
840           {
841             if (add_one)  // -(quotient + 1) == ~(quotient)
842               quotient.setInvert();
843             else
844               quotient.setNegative();
845           }
846         else if (add_one)
847           quotient.setAdd(1);
848       }
849     if (remainder != null)
850       {
851         // The remainder is by definition: X-Q*Y
852         remainder.set(ywords, rlen);
853         if (add_one)
854           {
855             // Subtract the remainder from Y:
856             // abs(R) = abs(Y) - abs(orig_rem) = -(abs(orig_rem) - abs(Y)).
857             BigInteger tmp;
858             if (y.words == null)
859               {
860                 tmp = remainder;
861                 tmp.set(yNegative ? ywords[0] + y.ival : ywords[0] - y.ival);
862               }
863             else
864               tmp = BigInteger.add(remainder, y, yNegative ? 1 : -1);
865             // Now tmp <= 0.
866             // In this case, abs(Q) = 1 + floor(abs(X)/abs(Y)).
867             // Hence, abs(Q*Y) > abs(X).
868             // So sign(remainder) = -sign(X).
869             if (xNegative)
870               remainder.setNegative(tmp);
871             else
872               remainder.set(tmp);
873           }
874         else
875           {
876             // If !add_one, then: abs(Q*Y) <= abs(X).
877             // So sign(remainder) = sign(X).
878             if (xNegative)
879               remainder.setNegative();
880           }
881       }
882   }
883
884   public BigInteger divide(BigInteger val)
885   {
886     if (val.isZero())
887       throw new ArithmeticException("divisor is zero");
888
889     BigInteger quot = new BigInteger();
890     divide(this, val, quot, null, TRUNCATE);
891     return quot.canonicalize();
892   }
893
894   public BigInteger remainder(BigInteger val)
895   {
896     if (val.isZero())
897       throw new ArithmeticException("divisor is zero");
898
899     BigInteger rem = new BigInteger();
900     divide(this, val, null, rem, TRUNCATE);
901     return rem.canonicalize();
902   }
903
904   public BigInteger[] divideAndRemainder(BigInteger val)
905   {
906     if (val.isZero())
907       throw new ArithmeticException("divisor is zero");
908
909     BigInteger[] result = new BigInteger[2];
910     result[0] = new BigInteger();
911     result[1] = new BigInteger();
912     divide(this, val, result[0], result[1], TRUNCATE);
913     result[0].canonicalize();
914     result[1].canonicalize();
915     return result;
916   }
917
918   public BigInteger mod(BigInteger m)
919   {
920     if (m.isNegative() || m.isZero())
921       throw new ArithmeticException("non-positive modulus");
922
923     BigInteger rem = new BigInteger();
924     divide(this, m, null, rem, FLOOR);
925     return rem.canonicalize();
926   }
927
928   /** Calculate power for BigInteger exponents.
929    * @param y exponent assumed to be non-negative. */
930   private BigInteger pow(BigInteger y)
931   {
932     if (isOne())
933       return this;
934     if (isMinusOne())
935       return y.isOdd () ? this : ONE;
936     if (y.words == null && y.ival >= 0)
937       return pow(y.ival);
938
939     // Assume exponent is non-negative.
940     if (isZero())
941       return this;
942
943     // Implemented by repeated squaring and multiplication.
944     BigInteger pow2 = this;
945     BigInteger r = null;
946     for (;;)  // for (i = 0;  ; i++)
947       {
948         // pow2 == x**(2**i)
949         // prod = x**(sum(j=0..i-1, (y>>j)&1))
950         if (y.isOdd())
951           r = r == null ? pow2 : times(r, pow2);  // r *= pow2
952         y = BigInteger.shift(y, -1);
953         if (y.isZero())
954           break;
955         // pow2 *= pow2;
956         pow2 = times(pow2, pow2);
957       }
958     return r == null ? ONE : r;
959   }
960
961   /** Calculate the integral power of a BigInteger.
962    * @param exponent the exponent (must be non-negative)
963    */
964   public BigInteger pow(int exponent)
965   {
966     if (exponent <= 0)
967       {
968         if (exponent == 0)
969           return ONE;
970         else
971           throw new ArithmeticException("negative exponent");
972       }
973     if (isZero())
974       return this;
975     int plen = words == null ? 1 : ival;  // Length of pow2.
976     int blen = ((bitLength() * exponent) >> 5) + 2 * plen;
977     boolean negative = isNegative() && (exponent & 1) != 0;
978     int[] pow2 = new int [blen];
979     int[] rwords = new int [blen];
980     int[] work = new int [blen];
981     getAbsolute(pow2);  // pow2 = abs(this);
982     int rlen = 1;
983     rwords[0] = 1; // rwords = 1;
984     for (;;)  // for (i = 0;  ; i++)
985       {
986         // pow2 == this**(2**i)
987         // prod = this**(sum(j=0..i-1, (exponent>>j)&1))
988         if ((exponent & 1) != 0)
989           { // r *= pow2
990             MPN.mul(work, pow2, plen, rwords, rlen);
991             int[] temp = work;  work = rwords;  rwords = temp;
992             rlen += plen;
993             while (rwords[rlen - 1] == 0)  rlen--;
994           }
995         exponent >>= 1;
996         if (exponent == 0)
997           break;
998         // pow2 *= pow2;
999         MPN.mul(work, pow2, plen, pow2, plen);
1000         int[] temp = work;  work = pow2;  pow2 = temp;  // swap to avoid a copy
1001         plen *= 2;
1002         while (pow2[plen - 1] == 0)  plen--;
1003       }
1004     if (rwords[rlen - 1] < 0)
1005       rlen++;
1006     if (negative)
1007       negate(rwords, rwords, rlen);
1008     return BigInteger.make(rwords, rlen);
1009   }
1010
1011   private static final int[] euclidInv(int a, int b, int prevDiv)
1012   {
1013     // Storage for return values, plus one slot for a temp int (see below).
1014     int[] xy;
1015
1016     if (b == 0)
1017       throw new ArithmeticException("not invertible");
1018     else if (b == 1)
1019       {
1020         // Success:  values are indeed invertible!
1021         // Bottom of the recursion reached; start unwinding.
1022         xy = new int[3];
1023         xy[0] = -prevDiv;
1024         xy[1] = 1;
1025         return xy;
1026       }
1027
1028     xy = euclidInv(b, a % b, a / b);    // Recursion happens here.
1029
1030     // xy[2] is just temp storage for intermediate results in the following
1031     // calculation.  This saves us a bit of space over having an int
1032     // allocated at every level of this recursive method.
1033     xy[2] = xy[0];
1034     xy[0] = xy[2] * -prevDiv + xy[1];
1035     xy[1] = xy[2];
1036     return xy;
1037   }
1038
1039   private static final BigInteger[]
1040     euclidInv(BigInteger a, BigInteger b, BigInteger prevDiv)
1041   {
1042     // FIXME: This method could be more efficient memory-wise and should be
1043     // modified as such since it is recursive.
1044
1045     // Storage for return values, plus one slot for a temp int (see below).
1046     BigInteger[] xy;
1047
1048     if (b.isZero())
1049       throw new ArithmeticException("not invertible");
1050     else if (b.isOne())
1051       {
1052         // Success:  values are indeed invertible!
1053         // Bottom of the recursion reached; start unwinding.
1054         xy = new BigInteger[3];
1055         xy[0] = neg(prevDiv);
1056         xy[1] = ONE;
1057         return xy;
1058       }
1059
1060     // Recursion happens in the following conditional!
1061
1062     // If a just contains an int, then use integer math for the rest.
1063     if (a.words == null)
1064       {
1065         int[] xyInt = euclidInv(b.ival, a.ival % b.ival, a.ival / b.ival);
1066         xy = new BigInteger[3];
1067         xy[0] = new BigInteger(xyInt[0]);
1068         xy[1] = new BigInteger(xyInt[1]);
1069       }
1070     else
1071       {
1072         BigInteger rem = new BigInteger();
1073         BigInteger quot = new BigInteger();
1074         divide(a, b, quot, rem, FLOOR);
1075         xy = euclidInv(b, rem, quot);
1076       }
1077
1078     // xy[2] is just temp storage for intermediate results in the following
1079     // calculation.  This saves us a bit of space over having a BigInteger
1080     // allocated at every level of this recursive method.
1081     xy[2] = xy[0];
1082     xy[0] = add(xy[1], times(xy[2], prevDiv), -1);
1083     xy[1] = xy[2];
1084     return xy;
1085   }
1086
1087   public BigInteger modInverse(BigInteger y)
1088   {
1089     if (y.isNegative() || y.isZero())
1090       throw new ArithmeticException("non-positive modulo");
1091
1092     // Degenerate cases.
1093     if (y.isOne())
1094       return ZERO;
1095     else if (isOne())
1096       return ONE;
1097
1098     // Use Euclid's algorithm as in gcd() but do this recursively
1099     // rather than in a loop so we can use the intermediate results as we
1100     // unwind from the recursion.
1101     // Used http://www.math.nmsu.edu/~crypto/EuclideanAlgo.html as reference.
1102     BigInteger result = new BigInteger();
1103     int xval = ival;
1104     int yval = y.ival;
1105     boolean swapped = false;
1106
1107     if (y.words == null)
1108       {
1109         // The result is guaranteed to be less than the modulus, y (which is
1110         // an int), so simplify this by working with the int result of this
1111         // modulo y.  Also, if this is negative, make it positive via modulo
1112         // math.  Note that BigInteger.mod() must be used even if this is
1113         // already an int as the % operator would provide a negative result if
1114         // this is negative, BigInteger.mod() never returns negative values.
1115         if (words != null || isNegative())
1116           xval = mod(y).ival;
1117
1118         // Swap values so x > y.
1119         if (yval > xval)
1120           {
1121             int tmp = xval; xval = yval; yval = tmp;
1122             swapped = true;
1123           }
1124         // Normally, the result is in the 2nd element of the array, but
1125         // if originally x < y, then x and y were swapped and the result
1126         // is in the 1st element of the array.
1127         result.ival =
1128           euclidInv(yval, xval % yval, xval / yval)[swapped ? 0 : 1];
1129
1130         // Result can't be negative, so make it positive by adding the
1131         // original modulus, y.ival (not the possibly "swapped" yval).
1132         if (result.ival < 0)
1133           result.ival += y.ival;
1134       }
1135     else
1136       {
1137         BigInteger x = this;
1138
1139         // As above, force this to be a positive value via modulo math.
1140         if (isNegative())
1141           x = mod(y);
1142
1143         // Swap values so x > y.
1144         if (x.compareTo(y) < 0)
1145           {
1146             BigInteger tmp = x; x = y; y = tmp;
1147             swapped = true;
1148           }
1149         // As above (for ints), result will be in the 2nd element unless
1150         // the original x and y were swapped.
1151         BigInteger rem = new BigInteger();
1152         BigInteger quot = new BigInteger();
1153         divide(x, y, quot, rem, FLOOR);
1154         result = euclidInv(y, rem, quot)[swapped ? 0 : 1];
1155
1156         // Result can't be negative, so make it positive by adding the
1157         // original modulus, y (which is now x if they were swapped).
1158         if (result.isNegative())
1159           result = add(result, swapped ? x : y, 1);
1160       }
1161     
1162     return result;
1163   }
1164
1165   public BigInteger modPow(BigInteger exponent, BigInteger m)
1166   {
1167     if (m.isNegative() || m.isZero())
1168       throw new ArithmeticException("non-positive modulo");
1169
1170     if (exponent.isNegative())
1171       return modInverse(m);
1172     if (exponent.isOne())
1173       return mod(m);
1174
1175     // To do this naively by first raising this to the power of exponent
1176     // and then performing modulo m would be extremely expensive, especially
1177     // for very large numbers.  The solution is found in Number Theory
1178     // where a combination of partial powers and modulos can be done easily.
1179     //
1180     // We'll use the algorithm for Additive Chaining which can be found on
1181     // p. 244 of "Applied Cryptography, Second Edition" by Bruce Schneier.
1182     BigInteger s, t, u;
1183     int i;
1184
1185     s = ONE;
1186     t = this;
1187     u = exponent;
1188
1189     while (!u.isZero())
1190       {
1191         if (u.and(ONE).isOne())
1192           s = times(s, t).mod(m);
1193         u = u.shiftRight(1);
1194         t = times(t, t).mod(m);
1195       }
1196
1197     return s;
1198   }
1199
1200   /** Calculate Greatest Common Divisor for non-negative ints. */
1201   private static final int gcd(int a, int b)
1202   {
1203     // Euclid's algorithm, copied from libg++.
1204     if (b > a)
1205       {
1206         int tmp = a; a = b; b = tmp;
1207       }
1208     for(;;)
1209       {
1210         if (b == 0)
1211           return a;
1212         else if (b == 1)
1213           return b;
1214         else
1215           {
1216             int tmp = b;
1217             b = a % b;
1218             a = tmp;
1219           }
1220       }
1221   }
1222
1223   public BigInteger gcd(BigInteger y)
1224   {
1225     int xval = ival;
1226     int yval = y.ival;
1227     if (words == null)
1228       {
1229         if (xval == 0)
1230           return BigInteger.abs(y);
1231         if (y.words == null
1232             && xval != Integer.MIN_VALUE && yval != Integer.MIN_VALUE)
1233           {
1234             if (xval < 0)
1235               xval = -xval;
1236             if (yval < 0)
1237               yval = -yval;
1238             return BigInteger.make(BigInteger.gcd(xval, yval));
1239           }
1240         xval = 1;
1241       }
1242     if (y.words == null)
1243       {
1244         if (yval == 0)
1245           return BigInteger.abs(this);
1246         yval = 1;
1247       }
1248     int len = (xval > yval ? xval : yval) + 1;
1249     int[] xwords = new int[len];
1250     int[] ywords = new int[len];
1251     getAbsolute(xwords);
1252     y.getAbsolute(ywords);
1253     len = MPN.gcd(xwords, ywords, len);
1254     BigInteger result = new BigInteger(0);
1255     result.ival = len;
1256     result.words = xwords;
1257     return result.canonicalize();
1258   }
1259
1260   public boolean isProbablePrime(int certainty)
1261   {
1262     /** We'll use the Rabin-Miller algorithm for doing a probabilistic
1263      * primality test.  It is fast, easy and has faster decreasing odds of a
1264      * composite passing than with other tests.  This means that this
1265      * method will actually have a probability much greater than the
1266      * 1 - .5^certainty specified in the JCL (p. 117), but I don't think
1267      * anyone will complain about better performance with greater certainty.
1268      *
1269      * The Rabin-Miller algorithm can be found on pp. 259-261 of "Applied
1270      * Cryptography, Second Edition" by Bruce Schneier.
1271      */
1272
1273     // First rule out small prime factors and assure the number is odd.
1274     for (int i = 0; i < primes.length; i++)
1275       {
1276         if (words == null && ival == primes[i])
1277           return true;
1278         if (remainder(make(primes[i])).isZero())
1279           return false;
1280       }
1281
1282     // Now perform the Rabin-Miller test.
1283     // NB: I know that this can be simplified programatically, but
1284     // I have tried to keep it as close as possible to the algorithm
1285     // as written in the Schneier book for reference purposes.
1286
1287     // Set b to the number of times 2 evenly divides (this - 1).
1288     // I.e. 2^b is the largest power of 2 that divides (this - 1).
1289     BigInteger pMinus1 = add(this, -1);
1290     int b = pMinus1.getLowestSetBit();
1291
1292     // Set m such that this = 1 + 2^b * m.
1293     BigInteger m = pMinus1.divide(make(2L << b - 1));
1294
1295     Random rand = new Random();
1296     while (certainty-- > 0)
1297       {
1298         // Pick a random number greater than 1 and less than this.
1299         // The algorithm says to pick a small number to make the calculations
1300         // go faster, but it doesn't say how small; we'll use 2 to 1024.
1301         int a = rand.nextInt();
1302         a = (a < 0 ? -a : a) % 1023 + 2;
1303
1304         BigInteger z = make(a).modPow(m, this);
1305         if (z.isOne() || z.equals(pMinus1))
1306           continue;                     // Passes the test; may be prime.
1307
1308         int i;
1309         for (i = 0; i < b; )
1310           {
1311             if (z.isOne())
1312               return false;
1313             i++;
1314             if (z.equals(pMinus1))
1315               break;                    // Passes the test; may be prime.
1316
1317             z = z.modPow(make(2), this);
1318           }
1319
1320         if (i == b && !z.equals(pMinus1))
1321           return false;
1322       }
1323     return true;
1324   }
1325
1326   private void setInvert()
1327   {
1328     if (words == null)
1329       ival = ~ival;
1330     else
1331       {
1332         for (int i = ival;  --i >= 0; )
1333           words[i] = ~words[i];
1334       }
1335   }
1336
1337   private void setShiftLeft(BigInteger x, int count)
1338   {
1339     int[] xwords;
1340     int xlen;
1341     if (x.words == null)
1342       {
1343         if (count < 32)
1344           {
1345             set((long) x.ival << count);
1346             return;
1347           }
1348         xwords = new int[1];
1349         xwords[0] = x.ival;
1350         xlen = 1;
1351       }
1352     else
1353       {
1354         xwords = x.words;
1355         xlen = x.ival;
1356       }
1357     int word_count = count >> 5;
1358     count &= 31;
1359     int new_len = xlen + word_count;
1360     if (count == 0)
1361       {
1362         realloc(new_len);
1363         for (int i = xlen;  --i >= 0; )
1364           words[i+word_count] = xwords[i];
1365       }
1366     else
1367       {
1368         new_len++;
1369         realloc(new_len);
1370         int shift_out = MPN.lshift(words, word_count, xwords, xlen, count);
1371         count = 32 - count;
1372         words[new_len-1] = (shift_out << count) >> count;  // sign-extend.
1373       }
1374     ival = new_len;
1375     for (int i = word_count;  --i >= 0; )
1376       words[i] = 0;
1377   }
1378
1379   private void setShiftRight(BigInteger x, int count)
1380   {
1381     if (x.words == null)
1382       set(count < 32 ? x.ival >> count : x.ival < 0 ? -1 : 0);
1383     else if (count == 0)
1384       set(x);
1385     else
1386       {
1387         boolean neg = x.isNegative();
1388         int word_count = count >> 5;
1389         count &= 31;
1390         int d_len = x.ival - word_count;
1391         if (d_len <= 0)
1392           set(neg ? -1 : 0);
1393         else
1394           {
1395             if (words == null || words.length < d_len)
1396               realloc(d_len);
1397             MPN.rshift0 (words, x.words, word_count, d_len, count);
1398             ival = d_len;
1399             if (neg)
1400               words[d_len-1] |= -2 << (31 - count);
1401           }
1402       }
1403   }
1404
1405   private void setShift(BigInteger x, int count)
1406   {
1407     if (count > 0)
1408       setShiftLeft(x, count);
1409     else
1410       setShiftRight(x, -count);
1411   }
1412
1413   private static BigInteger shift(BigInteger x, int count)
1414   {
1415     if (x.words == null)
1416       {
1417         if (count <= 0)
1418           return make(count > -32 ? x.ival >> (-count) : x.ival < 0 ? -1 : 0);
1419         if (count < 32)
1420           return make((long) x.ival << count);
1421       }
1422     if (count == 0)
1423       return x;
1424     BigInteger result = new BigInteger(0);
1425     result.setShift(x, count);
1426     return result.canonicalize();
1427   }
1428
1429   public BigInteger shiftLeft(int n)
1430   {
1431     return shift(this, n);
1432   }
1433
1434   public BigInteger shiftRight(int n)
1435   {
1436     return shift(this, -n);
1437   }
1438
1439   private void format(int radix, StringBuffer buffer)
1440   {
1441     if (words == null)
1442       buffer.append(Integer.toString(ival, radix));
1443     else if (ival <= 2)
1444       buffer.append(Long.toString(longValue(), radix));
1445     else
1446       {
1447         boolean neg = isNegative();
1448         int[] work;
1449         if (neg || radix != 16)
1450           {
1451             work = new int[ival];
1452             getAbsolute(work);
1453           }
1454         else
1455           work = words;
1456         int len = ival;
1457
1458         int buf_size = len * (MPN.chars_per_word(radix) + 1);
1459         if (radix == 16)
1460           {
1461             if (neg)
1462               buffer.append('-');
1463             int buf_start = buffer.length();
1464             for (int i = len;  --i >= 0; )
1465               {
1466                 int word = work[i];
1467                 for (int j = 8;  --j >= 0; )
1468                   {
1469                     int hex_digit = (word >> (4 * j)) & 0xF;
1470                     // Suppress leading zeros:
1471                     if (hex_digit > 0 || buffer.length() > buf_start)
1472                       buffer.append(Character.forDigit(hex_digit, 16));
1473                   }
1474               }
1475           }
1476         else
1477           {
1478             int i = buffer.length();
1479             for (;;)
1480               {
1481                 int digit = MPN.divmod_1(work, work, len, radix);
1482                 buffer.append(Character.forDigit(digit, radix));
1483                 while (len > 0 && work[len-1] == 0) len--;
1484                 if (len == 0)
1485                   break;
1486               }
1487             if (neg)
1488               buffer.append('-');
1489             /* Reverse buffer. */
1490             int j = buffer.length() - 1;
1491             while (i < j)
1492               {
1493                 char tmp = buffer.charAt(i);
1494                 buffer.setCharAt(i, buffer.charAt(j));
1495                 buffer.setCharAt(j, tmp);
1496                 i++;  j--;
1497               }
1498           }
1499       }
1500   }
1501
1502   public String toString()
1503   {
1504     return toString(10);
1505   }
1506
1507   public String toString(int radix)
1508   {
1509     if (words == null)
1510       return Integer.toString(ival, radix);
1511     else if (ival <= 2)
1512       return Long.toString(longValue(), radix);
1513     int buf_size = ival * (MPN.chars_per_word(radix) + 1);
1514     StringBuffer buffer = new StringBuffer(buf_size);
1515     format(radix, buffer);
1516     return buffer.toString();
1517   }
1518
1519   public int intValue()
1520   {
1521     if (words == null)
1522       return ival;
1523     return words[0];
1524   }
1525
1526   public long longValue()
1527   {
1528     if (words == null)
1529       return ival;
1530     if (ival == 1)
1531       return words[0];
1532     return ((long)words[1] << 32) + ((long)words[0] & 0xffffffffL);
1533   }
1534
1535   public int hashCode()
1536   {
1537     // FIXME: May not match hashcode of JDK.
1538     return words == null ? ival : (words[0] + words[ival - 1]);
1539   }
1540
1541   /* Assumes x and y are both canonicalized. */
1542   private static boolean equals(BigInteger x, BigInteger y)
1543   {
1544     if (x.words == null && y.words == null)
1545       return x.ival == y.ival;
1546     if (x.words == null || y.words == null || x.ival != y.ival)
1547       return false;
1548     for (int i = x.ival; --i >= 0; )
1549       {
1550         if (x.words[i] != y.words[i])
1551           return false;
1552       }
1553     return true;
1554   }
1555
1556   /* Assumes this and obj are both canonicalized. */
1557   public boolean equals(Object obj)
1558   {
1559     if (obj == null || ! (obj instanceof BigInteger))
1560       return false;
1561     return BigInteger.equals(this, (BigInteger) obj);
1562   }
1563
1564   private static BigInteger valueOf(String s, int radix)
1565        throws NumberFormatException
1566   {
1567     int len = s.length();
1568     // Testing (len < MPN.chars_per_word(radix)) would be more accurate,
1569     // but slightly more expensive, for little practical gain.
1570     if (len <= 15 && radix <= 16)
1571       return BigInteger.make(Long.parseLong(s, radix));
1572     
1573     int byte_len = 0;
1574     byte[] bytes = new byte[len];
1575     boolean negative = false;
1576     for (int i = 0;  i < len;  i++)
1577       {
1578         char ch = s.charAt(i);
1579         if (ch == '-')
1580           negative = true;
1581         else if (ch == '_' || (byte_len == 0 && (ch == ' ' || ch == '\t')))
1582           continue;
1583         else
1584           {
1585             int digit = Character.digit(ch, radix);
1586             if (digit < 0)
1587               break;
1588             bytes[byte_len++] = (byte) digit;
1589           }
1590       }
1591     return valueOf(bytes, byte_len, negative, radix);
1592   }
1593
1594   private static BigInteger valueOf(byte[] digits, int byte_len,
1595                                     boolean negative, int radix)
1596   {
1597     int chars_per_word = MPN.chars_per_word(radix);
1598     int[] words = new int[byte_len / chars_per_word + 1];
1599     int size = MPN.set_str(words, digits, byte_len, radix);
1600     if (size == 0)
1601       return ZERO;
1602     if (words[size-1] < 0)
1603       words[size++] = 0;
1604     if (negative)
1605       negate(words, words, size);
1606     return make(words, size);
1607   }
1608
1609   public double doubleValue()
1610   {
1611     if (words == null)
1612       return (double) ival;
1613     if (ival <= 2)
1614       return (double) longValue();
1615     if (isNegative())
1616       return BigInteger.neg(this).roundToDouble(0, true, false);
1617     else
1618       return roundToDouble(0, false, false);
1619   }
1620
1621   public float floatValue()
1622   {
1623     return (float) doubleValue();
1624   }
1625
1626   /** Return true if any of the lowest n bits are one.
1627    * (false if n is negative).  */
1628   private boolean checkBits(int n)
1629   {
1630     if (n <= 0)
1631       return false;
1632     if (words == null)
1633       return n > 31 || ((ival & ((1 << n) - 1)) != 0);
1634     int i;
1635     for (i = 0; i < (n >> 5) ; i++)
1636       if (words[i] != 0)
1637         return true;
1638     return (n & 31) != 0 && (words[i] & ((1 << (n & 31)) - 1)) != 0;
1639   }
1640
1641   /** Convert a semi-processed BigInteger to double.
1642    * Number must be non-negative.  Multiplies by a power of two, applies sign,
1643    * and converts to double, with the usual java rounding.
1644    * @param exp power of two, positive or negative, by which to multiply
1645    * @param neg true if negative
1646    * @param remainder true if the BigInteger is the result of a truncating
1647    * division that had non-zero remainder.  To ensure proper rounding in
1648    * this case, the BigInteger must have at least 54 bits.  */
1649   private double roundToDouble(int exp, boolean neg, boolean remainder)
1650   {
1651     // Compute length.
1652     int il = bitLength();
1653
1654     // Exponent when normalized to have decimal point directly after
1655     // leading one.  This is stored excess 1023 in the exponent bit field.
1656     exp += il - 1;
1657
1658     // Gross underflow.  If exp == -1075, we let the rounding
1659     // computation determine whether it is minval or 0 (which are just
1660     // 0x0000 0000 0000 0001 and 0x0000 0000 0000 0000 as bit
1661     // patterns).
1662     if (exp < -1075)
1663       return neg ? -0.0 : 0.0;
1664
1665     // gross overflow
1666     if (exp > 1023)
1667       return neg ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1668
1669     // number of bits in mantissa, including the leading one.
1670     // 53 unless it's denormalized
1671     int ml = (exp >= -1022 ? 53 : 53 + exp + 1022);
1672
1673     // Get top ml + 1 bits.  The extra one is for rounding.
1674     long m;
1675     int excess_bits = il - (ml + 1);
1676     if (excess_bits > 0)
1677       m = ((words == null) ? ival >> excess_bits
1678            : MPN.rshift_long(words, ival, excess_bits));
1679     else
1680       m = longValue() << (- excess_bits);
1681
1682     // Special rounding for maxval.  If the number exceeds maxval by
1683     // any amount, even if it's less than half a step, it overflows.
1684     if (exp == 1023 && ((m >> 1) == (1L << 53) - 1))
1685       {
1686         if (remainder || checkBits(il - ml))
1687           return neg ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1688         else
1689           return neg ? - Double.MAX_VALUE : Double.MAX_VALUE;
1690       }
1691
1692     // Normal round-to-even rule: round up if the bit dropped is a one, and
1693     // the bit above it or any of the bits below it is a one.
1694     if ((m & 1) == 1
1695         && ((m & 2) == 2 || remainder || checkBits(excess_bits)))
1696       {
1697         m += 2;
1698         // Check if we overflowed the mantissa
1699         if ((m & (1L << 54)) != 0)
1700           {
1701             exp++;
1702             // renormalize
1703             m >>= 1;
1704           }
1705         // Check if a denormalized mantissa was just rounded up to a
1706         // normalized one.
1707         else if (ml == 52 && (m & (1L << 53)) != 0)
1708           exp++;
1709       }
1710         
1711     // Discard the rounding bit
1712     m >>= 1;
1713
1714     long bits_sign = neg ? (1L << 63) : 0;
1715     exp += 1023;
1716     long bits_exp = (exp <= 0) ? 0 : ((long)exp) << 52;
1717     long bits_mant = m & ~(1L << 52);
1718     return Double.longBitsToDouble(bits_sign | bits_exp | bits_mant);
1719   }
1720
1721   /** Copy the abolute value of this into an array of words.
1722    * Assumes words.length >= (this.words == null ? 1 : this.ival).
1723    * Result is zero-extended, but need not be a valid 2's complement number.
1724    */
1725     
1726   private void getAbsolute(int[] words)
1727   {
1728     int len;
1729     if (this.words == null)
1730       {
1731         len = 1;
1732         words[0] = this.ival;
1733       }
1734     else
1735       {
1736         len = this.ival;
1737         for (int i = len;  --i >= 0; )
1738           words[i] = this.words[i];
1739       }
1740     if (words[len - 1] < 0)
1741       negate(words, words, len);
1742     for (int i = words.length;  --i > len; )
1743       words[i] = 0;
1744   }
1745
1746   /** Set dest[0:len-1] to the negation of src[0:len-1].
1747    * Return true if overflow (i.e. if src is -2**(32*len-1)).
1748    * Ok for src==dest. */
1749   private static boolean negate(int[] dest, int[] src, int len)
1750   {
1751     long carry = 1;
1752     boolean negative = src[len-1] < 0;
1753     for (int i = 0;  i < len;  i++)
1754       {
1755         carry += ((long) (~src[i]) & 0xffffffffL);
1756         dest[i] = (int) carry;
1757         carry >>= 32;
1758       }
1759     return (negative && dest[len-1] < 0);
1760   }
1761
1762   /** Destructively set this to the negative of x.
1763    * It is OK if x==this.*/
1764   private void setNegative(BigInteger x)
1765   {
1766     int len = x.ival;
1767     if (x.words == null)
1768       {
1769         if (len == Integer.MIN_VALUE)
1770           set(- (long) len);
1771         else
1772           set(-len);
1773         return;
1774       }
1775     realloc(len + 1);
1776     if (BigInteger.negate(words, x.words, len))
1777       words[len++] = 0;
1778     ival = len;
1779   }
1780
1781   /** Destructively negate this. */
1782   private final void setNegative()
1783   {
1784     setNegative(this);
1785   }
1786
1787   private static BigInteger abs(BigInteger x)
1788   {
1789     return x.isNegative() ? neg(x) : x;
1790   }
1791
1792   public BigInteger abs()
1793   {
1794     return abs(this);
1795   }
1796
1797   private static BigInteger neg(BigInteger x)
1798   {
1799     if (x.words == null && x.ival != Integer.MIN_VALUE)
1800       return make(- x.ival);
1801     BigInteger result = new BigInteger(0);
1802     result.setNegative(x);
1803     return result.canonicalize();
1804   }
1805
1806   public BigInteger negate()
1807   {
1808     return BigInteger.neg(this);
1809   }
1810
1811   /** Calculates ceiling(log2(this < 0 ? -this : this+1))
1812    * See Common Lisp: the Language, 2nd ed, p. 361.
1813    */
1814   public int bitLength()
1815   {
1816     if (words == null)
1817       return MPN.intLength(ival);
1818     else
1819       return MPN.intLength(words, ival);
1820   }
1821
1822   public byte[] toByteArray()
1823   {
1824     // Determine number of bytes needed.  The method bitlength returns
1825     // the size without the sign bit, so add one bit for that and then
1826     // add 7 more to emulate the ceil function using integer math.
1827     byte[] bytes = new byte[(bitLength() + 1 + 7) / 8];
1828     int nbytes = bytes.length;
1829
1830     int wptr = 0;
1831     int word;
1832
1833     // Deal with words array until one word or less is left to process.
1834     // If BigInteger is an int, then it is in ival and nbytes will be <= 4.
1835     while (nbytes > 4)
1836       {
1837         word = words[wptr++];
1838         for (int i = 4; i > 0; --i, word >>= 8)
1839           bytes[--nbytes] = (byte) word;
1840       }
1841
1842     // Deal with the last few bytes.  If BigInteger is an int, use ival.
1843     word = (words == null) ? ival : words[wptr];
1844     for ( ; nbytes > 0; word >>= 8)
1845       bytes[--nbytes] = (byte) word;
1846
1847     return bytes;
1848   }
1849
1850   /** Return the boolean opcode (for bitOp) for swapped operands.
1851    * I.e. bitOp(swappedOp(op), x, y) == bitOp(op, y, x).
1852    */
1853   private static int swappedOp(int op)
1854   {
1855     return
1856     "\000\001\004\005\002\003\006\007\010\011\014\015\012\013\016\017"
1857     .charAt(op);
1858   }
1859
1860   /** Do one the the 16 possible bit-wise operations of two BigIntegers. */
1861   private static BigInteger bitOp(int op, BigInteger x, BigInteger y)
1862   {
1863     switch (op)
1864       {
1865         case 0:  return ZERO;
1866         case 1:  return x.and(y);
1867         case 3:  return x;
1868         case 5:  return y;
1869         case 15: return make(-1);
1870       }
1871     BigInteger result = new BigInteger();
1872     setBitOp(result, op, x, y);
1873     return result.canonicalize();
1874   }
1875
1876   /** Do one the the 16 possible bit-wise operations of two BigIntegers. */
1877   private static void setBitOp(BigInteger result, int op,
1878                                BigInteger x, BigInteger y)
1879   {
1880     if (y.words == null) ;
1881     else if (x.words == null || x.ival < y.ival)
1882       {
1883         BigInteger temp = x;  x = y;  y = temp;
1884         op = swappedOp(op);
1885       }
1886     int xi;
1887     int yi;
1888     int xlen, ylen;
1889     if (y.words == null)
1890       {
1891         yi = y.ival;
1892         ylen = 1;
1893       }
1894     else
1895       {
1896         yi = y.words[0];
1897         ylen = y.ival;
1898       }
1899     if (x.words == null)
1900       {
1901         xi = x.ival;
1902         xlen = 1;
1903       }
1904     else
1905       {
1906         xi = x.words[0];
1907         xlen = x.ival;
1908       }
1909     if (xlen > 1)
1910       result.realloc(xlen);
1911     int[] w = result.words;
1912     int i = 0;
1913     // Code for how to handle the remainder of x.
1914     // 0:  Truncate to length of y.
1915     // 1:  Copy rest of x.
1916     // 2:  Invert rest of x.
1917     int finish = 0;
1918     int ni;
1919     switch (op)
1920       {
1921       case 0:  // clr
1922         ni = 0;
1923         break;
1924       case 1: // and
1925         for (;;)
1926           {
1927             ni = xi & yi;
1928             if (i+1 >= ylen) break;
1929             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1930           }
1931         if (yi < 0) finish = 1;
1932         break;
1933       case 2: // andc2
1934         for (;;)
1935           {
1936             ni = xi & ~yi;
1937             if (i+1 >= ylen) break;
1938             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1939           }
1940         if (yi >= 0) finish = 1;
1941         break;
1942       case 3:  // copy x
1943         ni = xi;
1944         finish = 1;  // Copy rest
1945         break;
1946       case 4: // andc1
1947         for (;;)
1948           {
1949             ni = ~xi & yi;
1950             if (i+1 >= ylen) break;
1951             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1952           }
1953         if (yi < 0) finish = 2;
1954         break;
1955       case 5: // copy y
1956         for (;;)
1957           {
1958             ni = yi;
1959             if (i+1 >= ylen) break;
1960             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1961           }
1962         break;
1963       case 6:  // xor
1964         for (;;)
1965           {
1966             ni = xi ^ yi;
1967             if (i+1 >= ylen) break;
1968             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1969           }
1970         finish = yi < 0 ? 2 : 1;
1971         break;
1972       case 7:  // ior
1973         for (;;)
1974           {
1975             ni = xi | yi;
1976             if (i+1 >= ylen) break;
1977             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1978           }
1979         if (yi >= 0) finish = 1;
1980         break;
1981       case 8:  // nor
1982         for (;;)
1983           {
1984             ni = ~(xi | yi);
1985             if (i+1 >= ylen) break;
1986             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1987           }
1988         if (yi >= 0)  finish = 2;
1989         break;
1990       case 9:  // eqv [exclusive nor]
1991         for (;;)
1992           {
1993             ni = ~(xi ^ yi);
1994             if (i+1 >= ylen) break;
1995             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
1996           }
1997         finish = yi >= 0 ? 2 : 1;
1998         break;
1999       case 10:  // c2
2000         for (;;)
2001           {
2002             ni = ~yi;
2003             if (i+1 >= ylen) break;
2004             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2005           }
2006         break;
2007       case 11:  // orc2
2008         for (;;)
2009           {
2010             ni = xi | ~yi;
2011             if (i+1 >= ylen) break;
2012             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2013           }
2014         if (yi < 0)  finish = 1;
2015         break;
2016       case 12:  // c1
2017         ni = ~xi;
2018         finish = 2;
2019         break;
2020       case 13:  // orc1
2021         for (;;)
2022           {
2023             ni = ~xi | yi;
2024             if (i+1 >= ylen) break;
2025             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2026           }
2027         if (yi >= 0) finish = 2;
2028         break;
2029       case 14:  // nand
2030         for (;;)
2031           {
2032             ni = ~(xi & yi);
2033             if (i+1 >= ylen) break;
2034             w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2035           }
2036         if (yi < 0) finish = 2;
2037         break;
2038       default:
2039       case 15:  // set
2040         ni = -1;
2041         break;
2042       }
2043     // Here i==ylen-1; w[0]..w[i-1] have the correct result;
2044     // and ni contains the correct result for w[i+1].
2045     if (i+1 == xlen)
2046       finish = 0;
2047     switch (finish)
2048       {
2049       case 0:
2050         if (i == 0 && w == null)
2051           {
2052             result.ival = ni;
2053             return;
2054           }
2055         w[i++] = ni;
2056         break;
2057       case 1:  w[i] = ni;  while (++i < xlen)  w[i] = x.words[i];  break;
2058       case 2:  w[i] = ni;  while (++i < xlen)  w[i] = ~x.words[i];  break;
2059       }
2060     result.ival = i;
2061   }
2062
2063   /** Return the logical (bit-wise) "and" of a BigInteger and an int. */
2064   private static BigInteger and(BigInteger x, int y)
2065   {
2066     if (x.words == null)
2067       return BigInteger.make(x.ival & y);
2068     if (y >= 0)
2069       return BigInteger.make(x.words[0] & y);
2070     int len = x.ival;
2071     int[] words = new int[len];
2072     words[0] = x.words[0] & y;
2073     while (--len > 0)
2074       words[len] = x.words[len];
2075     return BigInteger.make(words, x.ival);
2076   }
2077
2078   /** Return the logical (bit-wise) "and" of two BigIntegers. */
2079   public BigInteger and(BigInteger y)
2080   {
2081     if (y.words == null)
2082       return and(this, y.ival);
2083     else if (words == null)
2084       return and(y, ival);
2085
2086     BigInteger x = this;
2087     if (ival < y.ival)
2088       {
2089         BigInteger temp = this;  x = y;  y = temp;
2090       }
2091     int i;
2092     int len = y.isNegative() ? x.ival : y.ival;
2093     int[] words = new int[len];
2094     for (i = 0;  i < y.ival;  i++)
2095       words[i] = x.words[i] & y.words[i];
2096     for ( ; i < len;  i++)
2097       words[i] = x.words[i];
2098     return BigInteger.make(words, len);
2099   }
2100
2101   /** Return the logical (bit-wise) "(inclusive) or" of two BigIntegers. */
2102   public BigInteger or(BigInteger y)
2103   {
2104     return bitOp(7, this, y);
2105   }
2106
2107   /** Return the logical (bit-wise) "exclusive or" of two BigIntegers. */
2108   public BigInteger xor(BigInteger y)
2109   {
2110     return bitOp(6, this, y);
2111   }
2112
2113   /** Return the logical (bit-wise) negation of a BigInteger. */
2114   public BigInteger not()
2115   {
2116     return bitOp(12, this, ZERO);
2117   }
2118
2119   public BigInteger andNot(BigInteger val)
2120   {
2121     return and(val.not());
2122   }
2123
2124   public BigInteger clearBit(int n)
2125   {
2126     if (n < 0)
2127       throw new ArithmeticException();
2128
2129     return and(ONE.shiftLeft(n).not());
2130   }
2131
2132   public BigInteger setBit(int n)
2133   {
2134     if (n < 0)
2135       throw new ArithmeticException();
2136
2137     return or(ONE.shiftLeft(n));
2138   }
2139
2140   public boolean testBit(int n)
2141   {
2142     if (n < 0)
2143       throw new ArithmeticException();
2144
2145     return !and(ONE.shiftLeft(n)).isZero();
2146   }
2147
2148   public BigInteger flipBit(int n)
2149   {
2150     if (n < 0)
2151       throw new ArithmeticException();
2152
2153     return xor(ONE.shiftLeft(n));
2154   }
2155
2156   public int getLowestSetBit()
2157   {
2158     if (isZero())
2159       return -1;
2160
2161     if (words == null)
2162       return MPN.findLowestBit(ival);
2163     else
2164       return MPN.findLowestBit(words);
2165   }
2166
2167   // bit4count[I] is number of '1' bits in I.
2168   private static final byte[] bit4_count = { 0, 1, 1, 2,  1, 2, 2, 3,
2169                                              1, 2, 2, 3,  2, 3, 3, 4};
2170
2171   private static int bitCount(int i)
2172   {
2173     int count = 0;
2174     while (i != 0)
2175       {
2176         count += bit4_count[i & 15];
2177         i >>>= 4;
2178       }
2179     return count;
2180   }
2181
2182   private static int bitCount(int[] x, int len)
2183   {
2184     int count = 0;
2185     while (--len >= 0)
2186       count += bitCount(x[len]);
2187     return count;
2188   }
2189
2190   /** Count one bits in a BigInteger.
2191    * If argument is negative, count zero bits instead. */
2192   public int bitCount()
2193   {
2194     int i, x_len;
2195     int[] x_words = words;
2196     if (x_words == null)
2197       {
2198         x_len = 1;
2199         i = bitCount(ival);
2200       }
2201     else
2202       {
2203         x_len = ival;
2204         i = bitCount(x_words, x_len);
2205       }
2206     return isNegative() ? x_len * 32 - i : i;
2207   }
2208
2209   private void readObject(ObjectInputStream s)
2210     throws IOException, ClassNotFoundException
2211   {
2212     s.defaultReadObject();
2213     words = byteArrayToIntArray(magnitude, signum < 0 ? -1 : 0);
2214     BigInteger result = make(words, words.length);
2215     this.ival = result.ival;
2216     this.words = result.words;
2217   }
2218
2219   private void writeObject(ObjectOutputStream s)
2220     throws IOException, ClassNotFoundException
2221   {
2222     signum = signum();
2223     magnitude = toByteArray();
2224     s.defaultWriteObject();
2225   }
2226 }