001    /* java.math.BigInteger -- Arbitary precision integers
002       Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007  Free Software Foundation, Inc.
003    
004    This file is part of GNU Classpath.
005    
006    GNU Classpath is free software; you can redistribute it and/or modify
007    it under the terms of the GNU General Public License as published by
008    the Free Software Foundation; either version 2, or (at your option)
009    any later version.
010     
011    GNU Classpath is distributed in the hope that it will be useful, but
012    WITHOUT ANY WARRANTY; without even the implied warranty of
013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
014    General Public License for more details.
015    
016    You should have received a copy of the GNU General Public License
017    along with GNU Classpath; see the file COPYING.  If not, write to the
018    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
019    02110-1301 USA.
020    
021    Linking this library statically or dynamically with other modules is
022    making a combined work based on this library.  Thus, the terms and
023    conditions of the GNU General Public License cover the whole
024    combination.
025    
026    As a special exception, the copyright holders of this library give you
027    permission to link this library with independent modules to produce an
028    executable, regardless of the license terms of these independent
029    modules, and to copy and distribute the resulting executable under
030    terms of your choice, provided that you also meet, for each linked
031    independent module, the terms and conditions of the license of that
032    module.  An independent module is a module which is not derived from
033    or based on this library.  If you modify this library, you may extend
034    this exception to your version of the library, but you are not
035    obligated to do so.  If you do not wish to do so, delete this
036    exception statement from your version. */
037    
038    
039    package java.math;
040    
041    import gnu.classpath.Configuration;
042    
043    import gnu.java.lang.CPStringBuilder;
044    import gnu.java.math.GMP;
045    import gnu.java.math.MPN;
046    
047    import java.io.IOException;
048    import java.io.ObjectInputStream;
049    import java.io.ObjectOutputStream;
050    import java.util.Random;
051    import java.util.logging.Logger;
052    
053    /**
054     * Written using on-line Java Platform 1.2 API Specification, as well
055     * as "The Java Class Libraries", 2nd edition (Addison-Wesley, 1998) and
056     * "Applied Cryptography, Second Edition" by Bruce Schneier (Wiley, 1996).
057     * 
058     * Based primarily on IntNum.java BitOps.java by Per Bothner (per@bothner.com)
059     * (found in Kawa 1.6.62).
060     *
061     * @author Warren Levy (warrenl@cygnus.com)
062     * @date December 20, 1999.
063     * @status believed complete and correct.
064     */
065    public class BigInteger extends Number implements Comparable<BigInteger>
066    {
067      private static final Logger log = Logger.getLogger(BigInteger.class.getName());
068    
069      /** All integers are stored in 2's-complement form.
070       * If words == null, the ival is the value of this BigInteger.
071       * Otherwise, the first ival elements of words make the value
072       * of this BigInteger, stored in little-endian order, 2's-complement form. */
073      private transient int ival;
074      private transient int[] words;
075    
076      // Serialization fields.
077      // the first three, although not used in the code, are present for
078      // compatibility with older RI versions of this class. DO NOT REMOVE.
079      private int bitCount = -1;
080      private int bitLength = -1;
081      private int lowestSetBit = -2;
082      private byte[] magnitude;
083      private int signum;
084      private static final long serialVersionUID = -8287574255936472291L;
085    
086    
087      /** We pre-allocate integers in the range minFixNum..maxFixNum. 
088       * Note that we must at least preallocate 0, 1, and 10.  */
089      private static final int minFixNum = -100;
090      private static final int maxFixNum = 1024;
091      private static final int numFixNum = maxFixNum-minFixNum+1;
092      private static final BigInteger[] smallFixNums;
093      
094      /** The alter-ego GMP instance for this. */
095      private transient GMP mpz;
096    
097      private static final boolean USING_NATIVE = Configuration.WANT_NATIVE_BIG_INTEGER
098                                                  && initializeLibrary();
099    
100      static
101      {
102        if (USING_NATIVE)
103          {
104            smallFixNums = null;
105            ZERO = valueOf(0L);
106            ONE = valueOf(1L);
107            TEN = valueOf(10L);
108          }
109        else
110          {
111            smallFixNums = new BigInteger[numFixNum];
112            for (int i = numFixNum;  --i >= 0; )
113              smallFixNums[i] = new BigInteger(i + minFixNum);
114    
115            ZERO = smallFixNums[-minFixNum];
116            ONE = smallFixNums[1 - minFixNum];
117            TEN = smallFixNums[10 - minFixNum];
118          }
119      }
120    
121      /**
122       * The constant zero as a BigInteger.
123       * @since 1.2
124       */
125      public static final BigInteger ZERO;
126    
127      /**
128       * The constant one as a BigInteger.
129       * @since 1.2
130       */
131      public static final BigInteger ONE;
132    
133      /**
134       * The constant ten as a BigInteger.
135       * @since 1.5
136       */
137      public static final BigInteger TEN;
138    
139      /* Rounding modes: */
140      private static final int FLOOR = 1;
141      private static final int CEILING = 2;
142      private static final int TRUNCATE = 3;
143      private static final int ROUND = 4;
144    
145      /** When checking the probability of primes, it is most efficient to
146       * first check the factoring of small primes, so we'll use this array.
147       */
148      private static final int[] primes =
149        {   2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,
150           47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101, 103, 107,
151          109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
152          191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251 };
153    
154      /** HAC (Handbook of Applied Cryptography), Alfred Menezes & al. Table 4.4. */
155      private static final int[] k =
156          {100,150,200,250,300,350,400,500,600,800,1250, Integer.MAX_VALUE};
157      private static final int[] t =
158          { 27, 18, 15, 12,  9,  8,  7,  6,  5,  4,   3, 2};
159    
160      private BigInteger()
161      {
162        super();
163    
164        if (USING_NATIVE)
165          mpz = new GMP();
166      }
167    
168      /* Create a new (non-shared) BigInteger, and initialize to an int. */
169      private BigInteger(int value)
170      {
171        super();
172    
173        ival = value;
174      }
175    
176      public BigInteger(String s, int radix)
177      {
178        this();
179    
180        int len = s.length();
181        int i, digit;
182        boolean negative;
183        byte[] bytes;
184        char ch = s.charAt(0);
185        if (ch == '-')
186          {
187            negative = true;
188            i = 1;
189            bytes = new byte[len - 1];
190          }
191        else
192          {
193            negative = false;
194            i = 0;
195            bytes = new byte[len];
196          }
197        int byte_len = 0;
198        for ( ; i < len;  i++)
199          {
200            ch = s.charAt(i);
201            digit = Character.digit(ch, radix);
202            if (digit < 0)
203              throw new NumberFormatException("Invalid character at position #" + i);
204            bytes[byte_len++] = (byte) digit;
205          }
206    
207        if (USING_NATIVE)
208          {
209            bytes = null;
210            if (mpz.fromString(s, radix) != 0)
211              throw new NumberFormatException("String \"" + s
212                                              + "\" is NOT a valid number in base "
213                                              + radix);
214          }
215        else
216          {
217            BigInteger result;
218            // Testing (len < MPN.chars_per_word(radix)) would be more accurate,
219            // but slightly more expensive, for little practical gain.
220            if (len <= 15 && radix <= 16)
221              result = valueOf(Long.parseLong(s, radix));
222            else
223              result = valueOf(bytes, byte_len, negative, radix);
224    
225            this.ival = result.ival;
226            this.words = result.words;
227          }
228      }
229    
230      public BigInteger(String val)
231      {
232        this(val, 10);
233      }
234    
235      /* Create a new (non-shared) BigInteger, and initialize from a byte array. */
236      public BigInteger(byte[] val)
237      {
238        this();
239    
240        if (val == null || val.length < 1)
241          throw new NumberFormatException();
242    
243        if (USING_NATIVE)
244          mpz.fromByteArray(val);
245        else
246          {
247            words = byteArrayToIntArray(val, val[0] < 0 ? -1 : 0);
248            BigInteger result = make(words, words.length);
249            this.ival = result.ival;
250            this.words = result.words;
251          }
252      }
253    
254      public BigInteger(int signum, byte[] magnitude)
255      {
256        this();
257    
258        if (magnitude == null || signum > 1 || signum < -1)
259          throw new NumberFormatException();
260    
261        if (signum == 0)
262          {
263            int i;
264            for (i = magnitude.length - 1; i >= 0 && magnitude[i] == 0; --i)
265              ;
266            if (i >= 0)
267              throw new NumberFormatException();
268            return;
269          }
270    
271        if (USING_NATIVE)
272          mpz.fromSignedMagnitude(magnitude, signum == -1);
273        else
274          {
275            // Magnitude is always positive, so don't ever pass a sign of -1.
276            words = byteArrayToIntArray(magnitude, 0);
277            BigInteger result = make(words, words.length);
278            this.ival = result.ival;
279            this.words = result.words;
280    
281            if (signum < 0)
282              setNegative();
283          }
284      }
285    
286      public BigInteger(int numBits, Random rnd)
287      {
288        this();
289    
290        if (numBits < 0)
291          throw new IllegalArgumentException();
292    
293        init(numBits, rnd);
294      }
295    
296      private void init(int numBits, Random rnd)
297      {
298        if (USING_NATIVE)
299          {
300            int length = (numBits + 7) / 8;
301            byte[] magnitude = new byte[length];
302            rnd.nextBytes(magnitude);
303            int discardedBitCount = numBits % 8;
304            if (discardedBitCount != 0)
305              {
306                discardedBitCount = 8 - discardedBitCount;
307                magnitude[0] = (byte)((magnitude[0] & 0xFF) >>> discardedBitCount);
308              }
309            mpz.fromSignedMagnitude(magnitude, false);
310            magnitude = null;
311            return;
312          }
313    
314        int highbits = numBits & 31;
315        // minimum number of bytes to store the above number of bits
316        int highBitByteCount = (highbits + 7) / 8;
317        // number of bits to discard from the last byte
318        int discardedBitCount = highbits % 8;
319        if (discardedBitCount != 0)
320          discardedBitCount = 8 - discardedBitCount;
321        byte[] highBitBytes = new byte[highBitByteCount];
322        if (highbits > 0)
323          {
324            rnd.nextBytes(highBitBytes);
325            highbits = (highBitBytes[highBitByteCount - 1] & 0xFF) >>> discardedBitCount;
326            for (int i = highBitByteCount - 2; i >= 0; i--)
327              highbits = (highbits << 8) | (highBitBytes[i] & 0xFF);
328          }
329        int nwords = numBits / 32;
330    
331        while (highbits == 0 && nwords > 0)
332          {
333            highbits = rnd.nextInt();
334            --nwords;
335          }
336        if (nwords == 0 && highbits >= 0)
337          {
338            ival = highbits;
339          }
340        else
341          {
342            ival = highbits < 0 ? nwords + 2 : nwords + 1;
343            words = new int[ival];
344            words[nwords] = highbits;
345            while (--nwords >= 0)
346              words[nwords] = rnd.nextInt();
347          }
348      }
349    
350      public BigInteger(int bitLength, int certainty, Random rnd)
351      {
352        this();
353    
354        BigInteger result = new BigInteger();
355        while (true)
356          {
357            result.init(bitLength, rnd);
358            result = result.setBit(bitLength - 1);
359            if (result.isProbablePrime(certainty))
360              break;
361          }
362    
363        if (USING_NATIVE)
364          mpz.fromBI(result.mpz);
365        else
366          {
367            this.ival = result.ival;
368            this.words = result.words;
369          }
370      }
371    
372      /** 
373       *  Return a BigInteger that is bitLength bits long with a
374       *  probability < 2^-100 of being composite.
375       *
376       *  @param bitLength length in bits of resulting number
377       *  @param rnd random number generator to use
378       *  @throws ArithmeticException if bitLength < 2
379       *  @since 1.4
380       */
381      public static BigInteger probablePrime(int bitLength, Random rnd)
382      {
383        if (bitLength < 2)
384          throw new ArithmeticException();
385    
386        return new BigInteger(bitLength, 100, rnd);
387      }
388    
389      /** Return a (possibly-shared) BigInteger with a given long value. */
390      public static BigInteger valueOf(long val)
391      {
392        if (USING_NATIVE)
393          {
394            BigInteger result = new BigInteger();
395            result.mpz.fromLong(val);
396            return result;
397          }
398    
399        if (val >= minFixNum && val <= maxFixNum)
400          return smallFixNums[(int) val - minFixNum];
401        int i = (int) val;
402        if ((long) i == val)
403          return new BigInteger(i);
404        BigInteger result = alloc(2);
405        result.ival = 2;
406        result.words[0] = i;
407        result.words[1] = (int)(val >> 32);
408        return result;
409      }
410    
411      /**
412       * @return <code>true</code> if the GMP-based native implementation library
413       *         was successfully loaded. Returns <code>false</code> otherwise.
414       */
415      private static boolean initializeLibrary()
416      {
417        boolean result;
418        try
419        {
420          System.loadLibrary("javamath");
421          GMP.natInitializeLibrary();
422          result = true;
423        }
424        catch (Throwable x)
425        {
426          result = false;
427          if (Configuration.DEBUG)
428            {
429              log.info("Unable to use native BigInteger: " + x);
430              log.info("Will use a pure Java implementation instead");
431            }
432        }
433        return result;
434      }
435    
436      /** Make a canonicalized BigInteger from an array of words.
437       * The array may be reused (without copying). */
438      private static BigInteger make(int[] words, int len)
439      {
440        if (words == null)
441          return valueOf(len);
442        len = BigInteger.wordsNeeded(words, len);
443        if (len <= 1)
444          return len == 0 ? ZERO : valueOf(words[0]);
445        BigInteger num = new BigInteger();
446        num.words = words;
447        num.ival = len;
448        return num;
449      }
450    
451      /** Convert a big-endian byte array to a little-endian array of words. */
452      private static int[] byteArrayToIntArray(byte[] bytes, int sign)
453      {
454        // Determine number of words needed.
455        int[] words = new int[bytes.length/4 + 1];
456        int nwords = words.length;
457    
458        // Create a int out of modulo 4 high order bytes.
459        int bptr = 0;
460        int word = sign;
461        for (int i = bytes.length % 4; i > 0; --i, bptr++)
462          word = (word << 8) | (bytes[bptr] & 0xff);
463        words[--nwords] = word;
464    
465        // Elements remaining in byte[] are a multiple of 4.
466        while (nwords > 0)
467          words[--nwords] = bytes[bptr++] << 24 |
468                            (bytes[bptr++] & 0xff) << 16 |
469                            (bytes[bptr++] & 0xff) << 8 |
470                            (bytes[bptr++] & 0xff);
471        return words;
472      }
473    
474      /** Allocate a new non-shared BigInteger.
475       * @param nwords number of words to allocate
476       */
477      private static BigInteger alloc(int nwords)
478      {
479        BigInteger result = new BigInteger();
480        if (nwords > 1)
481        result.words = new int[nwords];
482        return result;
483      }
484    
485      /** Change words.length to nwords.
486       * We allow words.length to be upto nwords+2 without reallocating.
487       */
488      private void realloc(int nwords)
489      {
490        if (nwords == 0)
491          {
492            if (words != null)
493              {
494                if (ival > 0)
495                  ival = words[0];
496                words = null;
497              }
498          }
499        else if (words == null
500                 || words.length < nwords
501                 || words.length > nwords + 2)
502          {
503            int[] new_words = new int [nwords];
504            if (words == null)
505              {
506                new_words[0] = ival;
507                ival = 1;
508              }
509            else
510              {
511                if (nwords < ival)
512                  ival = nwords;
513                System.arraycopy(words, 0, new_words, 0, ival);
514              }
515            words = new_words;
516          }
517      }
518    
519      private boolean isNegative()
520      {
521        return (words == null ? ival : words[ival - 1]) < 0;
522      }
523    
524      public int signum()
525      {
526        if (USING_NATIVE)
527          return mpz.compare(ZERO.mpz);
528    
529        if (ival == 0 && words == null)
530          return 0;
531        int top = words == null ? ival : words[ival-1];
532        return top < 0 ? -1 : 1;
533      }
534    
535      private static int compareTo(BigInteger x, BigInteger y)
536      {
537        if (USING_NATIVE)
538          {
539            int dummy = y.signum; // force NPE check
540            return x.mpz.compare(y.mpz);
541          }
542    
543        if (x.words == null && y.words == null)
544          return x.ival < y.ival ? -1 : x.ival > y.ival ? 1 : 0;
545        boolean x_negative = x.isNegative();
546        boolean y_negative = y.isNegative();
547        if (x_negative != y_negative)
548          return x_negative ? -1 : 1;
549        int x_len = x.words == null ? 1 : x.ival;
550        int y_len = y.words == null ? 1 : y.ival;
551        if (x_len != y_len)
552          return (x_len > y_len) != x_negative ? 1 : -1;
553        return MPN.cmp(x.words, y.words, x_len);
554      }
555    
556      /** @since 1.2 */
557      public int compareTo(BigInteger val)
558      {
559        return compareTo(this, val);
560      }
561    
562      public BigInteger min(BigInteger val)
563      {
564        return compareTo(this, val) < 0 ? this : val;
565      }
566    
567      public BigInteger max(BigInteger val)
568      {
569        return compareTo(this, val) > 0 ? this : val;
570      }
571    
572      private boolean isZero()
573      {
574        return words == null && ival == 0;
575      }
576    
577      private boolean isOne()
578      {
579        return words == null && ival == 1;
580      }
581    
582      /** Calculate how many words are significant in words[0:len-1].
583       * Returns the least value x such that x>0 && words[0:x-1]==words[0:len-1],
584       * when words is viewed as a 2's complement integer.
585       */
586      private static int wordsNeeded(int[] words, int len)
587      {
588        int i = len;
589        if (i > 0)
590          {
591            int word = words[--i];
592            if (word == -1)
593              {
594                while (i > 0 && (word = words[i - 1]) < 0)
595                  {
596                    i--;
597                    if (word != -1) break;
598                  }
599              }
600            else
601              {
602                while (word == 0 && i > 0 && (word = words[i - 1]) >= 0)  i--;
603              }
604          }
605        return i + 1;
606      }
607    
608      private BigInteger canonicalize()
609      {
610        if (words != null
611            && (ival = BigInteger.wordsNeeded(words, ival)) <= 1)
612          {
613            if (ival == 1)
614              ival = words[0];
615            words = null;
616          }
617        if (words == null && ival >= minFixNum && ival <= maxFixNum)
618          return smallFixNums[ival - minFixNum];
619        return this;
620      }
621    
622      /** Add two ints, yielding a BigInteger. */
623      private static BigInteger add(int x, int y)
624      {
625        return valueOf((long) x + (long) y);
626      }
627    
628      /** Add a BigInteger and an int, yielding a new BigInteger. */
629      private static BigInteger add(BigInteger x, int y)
630      {
631        if (x.words == null)
632          return BigInteger.add(x.ival, y);
633        BigInteger result = new BigInteger(0);
634        result.setAdd(x, y);
635        return result.canonicalize();
636      }
637    
638      /** Set this to the sum of x and y.
639       * OK if x==this. */
640      private void setAdd(BigInteger x, int y)
641      {
642        if (x.words == null)
643          {
644            set((long) x.ival + (long) y);
645            return;
646          }
647        int len = x.ival;
648        realloc(len + 1);
649        long carry = y;
650        for (int i = 0;  i < len;  i++)
651          {
652            carry += ((long) x.words[i] & 0xffffffffL);
653            words[i] = (int) carry;
654            carry >>= 32;
655          }
656        if (x.words[len - 1] < 0)
657          carry--;
658        words[len] = (int) carry;
659        ival = wordsNeeded(words, len + 1);
660      }
661    
662      /** Destructively add an int to this. */
663      private void setAdd(int y)
664      {
665        setAdd(this, y);
666      }
667    
668      /** Destructively set the value of this to a long. */
669      private void set(long y)
670      {
671        int i = (int) y;
672        if ((long) i == y)
673          {
674            ival = i;
675            words = null;
676          }
677        else
678          {
679            realloc(2);
680            words[0] = i;
681            words[1] = (int) (y >> 32);
682            ival = 2;
683          }
684      }
685    
686      /** Destructively set the value of this to the given words.
687      * The words array is reused, not copied. */
688      private void set(int[] words, int length)
689      {
690        this.ival = length;
691        this.words = words;
692      }
693    
694      /** Destructively set the value of this to that of y. */
695      private void set(BigInteger y)
696      {
697        if (y.words == null)
698          set(y.ival);
699        else if (this != y)
700          {
701            realloc(y.ival);
702            System.arraycopy(y.words, 0, words, 0, y.ival);
703            ival = y.ival;
704          }
705      }
706    
707      /** Add two BigIntegers, yielding their sum as another BigInteger. */
708      private static BigInteger add(BigInteger x, BigInteger y, int k)
709      {
710        if (x.words == null && y.words == null)
711          return valueOf((long) k * (long) y.ival + (long) x.ival);
712        if (k != 1)
713          {
714            if (k == -1)
715              y = BigInteger.neg(y);
716            else
717              y = BigInteger.times(y, valueOf(k));
718          }
719        if (x.words == null)
720          return BigInteger.add(y, x.ival);
721        if (y.words == null)
722          return BigInteger.add(x, y.ival);
723        // Both are big
724        if (y.ival > x.ival)
725          { // Swap so x is longer then y.
726            BigInteger tmp = x;  x = y;  y = tmp;
727          }
728        BigInteger result = alloc(x.ival + 1);
729        int i = y.ival;
730        long carry = MPN.add_n(result.words, x.words, y.words, i);
731        long y_ext = y.words[i - 1] < 0 ? 0xffffffffL : 0;
732        for (; i < x.ival;  i++)
733          {
734            carry += ((long) x.words[i] & 0xffffffffL) + y_ext;
735            result.words[i] = (int) carry;
736            carry >>>= 32;
737          }
738        if (x.words[i - 1] < 0)
739          y_ext--;
740        result.words[i] = (int) (carry + y_ext);
741        result.ival = i+1;
742        return result.canonicalize();
743      }
744    
745      public BigInteger add(BigInteger val)
746      {
747        if (USING_NATIVE)
748          {
749            int dummy = val.signum; // force NPE check
750            BigInteger result = new BigInteger();
751            mpz.add(val.mpz, result.mpz);
752            return result;
753          }
754    
755        return add(this, val, 1);
756      }
757    
758      public BigInteger subtract(BigInteger val)
759      {
760        if (USING_NATIVE)
761          {
762            int dummy = val.signum; // force NPE check
763            BigInteger result = new BigInteger();
764            mpz.subtract(val.mpz, result.mpz);
765            return result;
766          }
767    
768        return add(this, val, -1);
769      }
770    
771      private static BigInteger times(BigInteger x, int y)
772      {
773        if (y == 0)
774          return ZERO;
775        if (y == 1)
776          return x;
777        int[] xwords = x.words;
778        int xlen = x.ival;
779        if (xwords == null)
780          return valueOf((long) xlen * (long) y);
781        boolean negative;
782        BigInteger result = BigInteger.alloc(xlen + 1);
783        if (xwords[xlen - 1] < 0)
784          {
785            negative = true;
786            negate(result.words, xwords, xlen);
787            xwords = result.words;
788          }
789        else
790          negative = false;
791        if (y < 0)
792          {
793            negative = !negative;
794            y = -y;
795          }
796        result.words[xlen] = MPN.mul_1(result.words, xwords, xlen, y);
797        result.ival = xlen + 1;
798        if (negative)
799          result.setNegative();
800        return result.canonicalize();
801      }
802    
803      private static BigInteger times(BigInteger x, BigInteger y)
804      {
805        if (y.words == null)
806          return times(x, y.ival);
807        if (x.words == null)
808          return times(y, x.ival);
809        boolean negative = false;
810        int[] xwords;
811        int[] ywords;
812        int xlen = x.ival;
813        int ylen = y.ival;
814        if (x.isNegative())
815          {
816            negative = true;
817            xwords = new int[xlen];
818            negate(xwords, x.words, xlen);
819          }
820        else
821          {
822            negative = false;
823            xwords = x.words;
824          }
825        if (y.isNegative())
826          {
827            negative = !negative;
828            ywords = new int[ylen];
829            negate(ywords, y.words, ylen);
830          }
831        else
832          ywords = y.words;
833        // Swap if x is shorter then y.
834        if (xlen < ylen)
835          {
836            int[] twords = xwords;  xwords = ywords;  ywords = twords;
837            int tlen = xlen;  xlen = ylen;  ylen = tlen;
838          }
839        BigInteger result = BigInteger.alloc(xlen+ylen);
840        MPN.mul(result.words, xwords, xlen, ywords, ylen);
841        result.ival = xlen+ylen;
842        if (negative)
843          result.setNegative();
844        return result.canonicalize();
845      }
846    
847      public BigInteger multiply(BigInteger y)
848      {
849        if (USING_NATIVE)
850          {
851            int dummy = y.signum; // force NPE check
852            BigInteger result = new BigInteger();
853            mpz.multiply(y.mpz, result.mpz);
854            return result;
855          }
856    
857        return times(this, y);
858      }
859    
860      private static void divide(long x, long y,
861                                 BigInteger quotient, BigInteger remainder,
862                                 int rounding_mode)
863      {
864        boolean xNegative, yNegative;
865        if (x < 0)
866          {
867            xNegative = true;
868            if (x == Long.MIN_VALUE)
869              {
870                divide(valueOf(x), valueOf(y),
871                       quotient, remainder, rounding_mode);
872                return;
873              }
874            x = -x;
875          }
876        else
877          xNegative = false;
878    
879        if (y < 0)
880          {
881            yNegative = true;
882            if (y == Long.MIN_VALUE)
883              {
884                if (rounding_mode == TRUNCATE)
885                  { // x != Long.Min_VALUE implies abs(x) < abs(y)
886                    if (quotient != null)
887                      quotient.set(0);
888                    if (remainder != null)
889                      remainder.set(x);
890                  }
891                else
892                  divide(valueOf(x), valueOf(y),
893                          quotient, remainder, rounding_mode);
894                return;
895              }
896            y = -y;
897          }
898        else
899          yNegative = false;
900    
901        long q = x / y;
902        long r = x % y;
903        boolean qNegative = xNegative ^ yNegative;
904    
905        boolean add_one = false;
906        if (r != 0)
907          {
908            switch (rounding_mode)
909              {
910              case TRUNCATE:
911                break;
912              case CEILING:
913              case FLOOR:
914                if (qNegative == (rounding_mode == FLOOR))
915                  add_one = true;
916                break;
917              case ROUND:
918                add_one = r > ((y - (q & 1)) >> 1);
919                break;
920              }
921          }
922        if (quotient != null)
923          {
924            if (add_one)
925              q++;
926            if (qNegative)
927              q = -q;
928            quotient.set(q);
929          }
930        if (remainder != null)
931          {
932            // The remainder is by definition: X-Q*Y
933            if (add_one)
934              {
935                // Subtract the remainder from Y.
936                r = y - r;
937                // In this case, abs(Q*Y) > abs(X).
938                // So sign(remainder) = -sign(X).
939                xNegative = ! xNegative;
940              }
941            else
942              {
943                // If !add_one, then: abs(Q*Y) <= abs(X).
944                // So sign(remainder) = sign(X).
945              }
946            if (xNegative)
947              r = -r;
948            remainder.set(r);
949          }
950      }
951    
952      /** Divide two integers, yielding quotient and remainder.
953       * @param x the numerator in the division
954       * @param y the denominator in the division
955       * @param quotient is set to the quotient of the result (iff quotient!=null)
956       * @param remainder is set to the remainder of the result
957       *  (iff remainder!=null)
958       * @param rounding_mode one of FLOOR, CEILING, TRUNCATE, or ROUND.
959       */
960      private static void divide(BigInteger x, BigInteger y,
961                                 BigInteger quotient, BigInteger remainder,
962                                 int rounding_mode)
963      {
964        if ((x.words == null || x.ival <= 2)
965            && (y.words == null || y.ival <= 2))
966          {
967            long x_l = x.longValue();
968            long y_l = y.longValue();
969            if (x_l != Long.MIN_VALUE && y_l != Long.MIN_VALUE)
970              {
971                divide(x_l, y_l, quotient, remainder, rounding_mode);
972                return;
973              }
974          }
975    
976        boolean xNegative = x.isNegative();
977        boolean yNegative = y.isNegative();
978        boolean qNegative = xNegative ^ yNegative;
979    
980        int ylen = y.words == null ? 1 : y.ival;
981        int[] ywords = new int[ylen];
982        y.getAbsolute(ywords);
983        while (ylen > 1 && ywords[ylen - 1] == 0)  ylen--;
984    
985        int xlen = x.words == null ? 1 : x.ival;
986        int[] xwords = new int[xlen+2];
987        x.getAbsolute(xwords);
988        while (xlen > 1 && xwords[xlen-1] == 0)  xlen--;
989    
990        int qlen, rlen;
991    
992        int cmpval = MPN.cmp(xwords, xlen, ywords, ylen);
993        if (cmpval < 0)  // abs(x) < abs(y)
994          { // quotient = 0;  remainder = num.
995            int[] rwords = xwords;  xwords = ywords;  ywords = rwords;
996            rlen = xlen;  qlen = 1;  xwords[0] = 0;
997          }
998        else if (cmpval == 0)  // abs(x) == abs(y)
999          {
1000            xwords[0] = 1;  qlen = 1;  // quotient = 1
1001            ywords[0] = 0;  rlen = 1;  // remainder = 0;
1002          }
1003        else if (ylen == 1)
1004          {
1005            qlen = xlen;
1006            // Need to leave room for a word of leading zeros if dividing by 1
1007            // and the dividend has the high bit set.  It might be safe to
1008            // increment qlen in all cases, but it certainly is only necessary
1009            // in the following case.
1010            if (ywords[0] == 1 && xwords[xlen-1] < 0)
1011              qlen++;
1012            rlen = 1;
1013            ywords[0] = MPN.divmod_1(xwords, xwords, xlen, ywords[0]);
1014          }
1015        else  // abs(x) > abs(y)
1016          {
1017            // Normalize the denominator, i.e. make its most significant bit set by
1018            // shifting it normalization_steps bits to the left.  Also shift the
1019            // numerator the same number of steps (to keep the quotient the same!).
1020    
1021            int nshift = MPN.count_leading_zeros(ywords[ylen - 1]);
1022            if (nshift != 0)
1023              {
1024                // Shift up the denominator setting the most significant bit of
1025                // the most significant word.
1026                MPN.lshift(ywords, 0, ywords, ylen, nshift);
1027    
1028                // Shift up the numerator, possibly introducing a new most
1029                // significant word.
1030                int x_high = MPN.lshift(xwords, 0, xwords, xlen, nshift);
1031                xwords[xlen++] = x_high;
1032              }
1033    
1034            if (xlen == ylen)
1035              xwords[xlen++] = 0;
1036            MPN.divide(xwords, xlen, ywords, ylen);
1037            rlen = ylen;
1038            MPN.rshift0 (ywords, xwords, 0, rlen, nshift);
1039    
1040            qlen = xlen + 1 - ylen;
1041            if (quotient != null)
1042              {
1043                for (int i = 0;  i < qlen;  i++)
1044                  xwords[i] = xwords[i+ylen];
1045              }
1046          }
1047    
1048        if (ywords[rlen-1] < 0)
1049          {
1050            ywords[rlen] = 0;
1051            rlen++;
1052          }
1053    
1054        // Now the quotient is in xwords, and the remainder is in ywords.
1055    
1056        boolean add_one = false;
1057        if (rlen > 1 || ywords[0] != 0)
1058          { // Non-zero remainder i.e. in-exact quotient.
1059            switch (rounding_mode)
1060              {
1061              case TRUNCATE:
1062                break;
1063              case CEILING:
1064              case FLOOR:
1065                if (qNegative == (rounding_mode == FLOOR))
1066                  add_one = true;
1067                break;
1068              case ROUND:
1069                // int cmp = compareTo(remainder<<1, abs(y));
1070                BigInteger tmp = remainder == null ? new BigInteger() : remainder;
1071                tmp.set(ywords, rlen);
1072                tmp = shift(tmp, 1);
1073                if (yNegative)
1074                  tmp.setNegative();
1075                int cmp = compareTo(tmp, y);
1076                // Now cmp == compareTo(sign(y)*(remainder<<1), y)
1077                if (yNegative)
1078                  cmp = -cmp;
1079                add_one = (cmp == 1) || (cmp == 0 && (xwords[0]&1) != 0);
1080              }
1081          }
1082        if (quotient != null)
1083          {
1084            quotient.set(xwords, qlen);
1085            if (qNegative)
1086              {
1087                if (add_one)  // -(quotient + 1) == ~(quotient)
1088                  quotient.setInvert();
1089                else
1090                  quotient.setNegative();
1091              }
1092            else if (add_one)
1093              quotient.setAdd(1);
1094          }
1095        if (remainder != null)
1096          {
1097            // The remainder is by definition: X-Q*Y
1098            remainder.set(ywords, rlen);
1099            if (add_one)
1100              {
1101                // Subtract the remainder from Y:
1102                // abs(R) = abs(Y) - abs(orig_rem) = -(abs(orig_rem) - abs(Y)).
1103                BigInteger tmp;
1104                if (y.words == null)
1105                  {
1106                    tmp = remainder;
1107                    tmp.set(yNegative ? ywords[0] + y.ival : ywords[0] - y.ival);
1108                  }
1109                else
1110                  tmp = BigInteger.add(remainder, y, yNegative ? 1 : -1);
1111                // Now tmp <= 0.
1112                // In this case, abs(Q) = 1 + floor(abs(X)/abs(Y)).
1113                // Hence, abs(Q*Y) > abs(X).
1114                // So sign(remainder) = -sign(X).
1115                if (xNegative)
1116                  remainder.setNegative(tmp);
1117                else
1118                  remainder.set(tmp);
1119              }
1120            else
1121              {
1122                // If !add_one, then: abs(Q*Y) <= abs(X).
1123                // So sign(remainder) = sign(X).
1124                if (xNegative)
1125                  remainder.setNegative();
1126              }
1127          }
1128      }
1129    
1130      public BigInteger divide(BigInteger val)
1131      {
1132        if (USING_NATIVE)
1133          {
1134            if (val.compareTo(ZERO) == 0)
1135              throw new ArithmeticException("divisor is zero");
1136    
1137            BigInteger result = new BigInteger();
1138            mpz.quotient(val.mpz, result.mpz);
1139            return result;
1140          }
1141    
1142        if (val.isZero())
1143          throw new ArithmeticException("divisor is zero");
1144    
1145        BigInteger quot = new BigInteger();
1146        divide(this, val, quot, null, TRUNCATE);
1147        return quot.canonicalize();
1148      }
1149    
1150      public BigInteger remainder(BigInteger val)
1151      {
1152        if (USING_NATIVE)
1153          {
1154            if (val.compareTo(ZERO) == 0)
1155              throw new ArithmeticException("divisor is zero");
1156    
1157            BigInteger result = new BigInteger();
1158            mpz.remainder(val.mpz, result.mpz);
1159            return result;
1160          }
1161    
1162        if (val.isZero())
1163          throw new ArithmeticException("divisor is zero");
1164    
1165        BigInteger rem = new BigInteger();
1166        divide(this, val, null, rem, TRUNCATE);
1167        return rem.canonicalize();
1168      }
1169    
1170      public BigInteger[] divideAndRemainder(BigInteger val)
1171      {
1172        if (USING_NATIVE)
1173          {
1174            if (val.compareTo(ZERO) == 0)
1175              throw new ArithmeticException("divisor is zero");
1176    
1177            BigInteger q = new BigInteger();
1178            BigInteger r = new BigInteger();
1179            mpz.quotientAndRemainder(val.mpz, q.mpz, r.mpz);
1180            return new BigInteger[] { q, r };
1181          }
1182    
1183        if (val.isZero())
1184          throw new ArithmeticException("divisor is zero");
1185    
1186        BigInteger[] result = new BigInteger[2];
1187        result[0] = new BigInteger();
1188        result[1] = new BigInteger();
1189        divide(this, val, result[0], result[1], TRUNCATE);
1190        result[0].canonicalize();
1191        result[1].canonicalize();
1192        return result;
1193      }
1194    
1195      public BigInteger mod(BigInteger m)
1196      {
1197        if (USING_NATIVE)
1198          {
1199            int dummy = m.signum; // force NPE check
1200            if (m.compareTo(ZERO) < 1)
1201              throw new ArithmeticException("non-positive modulus");
1202    
1203            BigInteger result = new BigInteger();
1204            mpz.modulo(m.mpz, result.mpz);
1205            return result;
1206          }
1207    
1208        if (m.isNegative() || m.isZero())
1209          throw new ArithmeticException("non-positive modulus");
1210    
1211        BigInteger rem = new BigInteger();
1212        divide(this, m, null, rem, FLOOR);
1213        return rem.canonicalize();
1214      }
1215    
1216      /** Calculate the integral power of a BigInteger.
1217       * @param exponent the exponent (must be non-negative)
1218       */
1219      public BigInteger pow(int exponent)
1220      {
1221        if (exponent <= 0)
1222          {
1223            if (exponent == 0)
1224              return ONE;
1225              throw new ArithmeticException("negative exponent");
1226          }
1227    
1228        if (USING_NATIVE)
1229          {
1230            BigInteger result = new BigInteger();
1231            mpz.pow(exponent, result.mpz);
1232            return result;
1233          }
1234    
1235        if (isZero())
1236          return this;
1237        int plen = words == null ? 1 : ival;  // Length of pow2.
1238        int blen = ((bitLength() * exponent) >> 5) + 2 * plen;
1239        boolean negative = isNegative() && (exponent & 1) != 0;
1240        int[] pow2 = new int [blen];
1241        int[] rwords = new int [blen];
1242        int[] work = new int [blen];
1243        getAbsolute(pow2);  // pow2 = abs(this);
1244        int rlen = 1;
1245        rwords[0] = 1; // rwords = 1;
1246        for (;;)  // for (i = 0;  ; i++)
1247          {
1248            // pow2 == this**(2**i)
1249            // prod = this**(sum(j=0..i-1, (exponent>>j)&1))
1250            if ((exponent & 1) != 0)
1251              { // r *= pow2
1252                MPN.mul(work, pow2, plen, rwords, rlen);
1253                int[] temp = work;  work = rwords;  rwords = temp;
1254                rlen += plen;
1255                while (rwords[rlen - 1] == 0)  rlen--;
1256              }
1257            exponent >>= 1;
1258            if (exponent == 0)
1259              break;
1260            // pow2 *= pow2;
1261            MPN.mul(work, pow2, plen, pow2, plen);
1262            int[] temp = work;  work = pow2;  pow2 = temp;  // swap to avoid a copy
1263            plen *= 2;
1264            while (pow2[plen - 1] == 0)  plen--;
1265          }
1266        if (rwords[rlen - 1] < 0)
1267          rlen++;
1268        if (negative)
1269          negate(rwords, rwords, rlen);
1270        return BigInteger.make(rwords, rlen);
1271      }
1272    
1273      private static int[] euclidInv(int a, int b, int prevDiv)
1274      {
1275        if (b == 0)
1276          throw new ArithmeticException("not invertible");
1277    
1278        if (b == 1)
1279            // Success:  values are indeed invertible!
1280            // Bottom of the recursion reached; start unwinding.
1281            return new int[] { -prevDiv, 1 };
1282    
1283        int[] xy = euclidInv(b, a % b, a / b);      // Recursion happens here.
1284        a = xy[0]; // use our local copy of 'a' as a work var
1285        xy[0] = a * -prevDiv + xy[1];
1286        xy[1] = a;
1287        return xy;
1288      }
1289    
1290      private static void euclidInv(BigInteger a, BigInteger b,
1291                                    BigInteger prevDiv, BigInteger[] xy)
1292      {
1293        if (b.isZero())
1294          throw new ArithmeticException("not invertible");
1295    
1296        if (b.isOne())
1297          {
1298            // Success:  values are indeed invertible!
1299            // Bottom of the recursion reached; start unwinding.
1300            xy[0] = neg(prevDiv);
1301            xy[1] = ONE;
1302            return;
1303          }
1304    
1305        // Recursion happens in the following conditional!
1306    
1307        // If a just contains an int, then use integer math for the rest.
1308        if (a.words == null)
1309          {
1310            int[] xyInt = euclidInv(b.ival, a.ival % b.ival, a.ival / b.ival);
1311            xy[0] = new BigInteger(xyInt[0]);
1312            xy[1] = new BigInteger(xyInt[1]);
1313          }
1314        else
1315          {
1316            BigInteger rem = new BigInteger();
1317            BigInteger quot = new BigInteger();
1318            divide(a, b, quot, rem, FLOOR);
1319            // quot and rem may not be in canonical form. ensure
1320            rem.canonicalize();
1321            quot.canonicalize();
1322            euclidInv(b, rem, quot, xy);
1323          }
1324    
1325        BigInteger t = xy[0];
1326        xy[0] = add(xy[1], times(t, prevDiv), -1);
1327        xy[1] = t;
1328      }
1329    
1330      public BigInteger modInverse(BigInteger y)
1331      {
1332        if (USING_NATIVE)
1333          {
1334            int dummy = y.signum; // force NPE check
1335            if (mpz.compare(ZERO.mpz) < 1)
1336              throw new ArithmeticException("non-positive modulo");
1337    
1338            BigInteger result = new BigInteger();
1339            mpz.modInverse(y.mpz, result.mpz);
1340            return result;
1341          }
1342    
1343        if (y.isNegative() || y.isZero())
1344          throw new ArithmeticException("non-positive modulo");
1345    
1346        // Degenerate cases.
1347        if (y.isOne())
1348          return ZERO;
1349        if (isOne())
1350          return ONE;
1351    
1352        // Use Euclid's algorithm as in gcd() but do this recursively
1353        // rather than in a loop so we can use the intermediate results as we
1354        // unwind from the recursion.
1355        // Used http://www.math.nmsu.edu/~crypto/EuclideanAlgo.html as reference.
1356        BigInteger result = new BigInteger();
1357        boolean swapped = false;
1358    
1359        if (y.words == null)
1360          {
1361            // The result is guaranteed to be less than the modulus, y (which is
1362            // an int), so simplify this by working with the int result of this
1363            // modulo y.  Also, if this is negative, make it positive via modulo
1364            // math.  Note that BigInteger.mod() must be used even if this is
1365            // already an int as the % operator would provide a negative result if
1366            // this is negative, BigInteger.mod() never returns negative values.
1367            int xval = (words != null || isNegative()) ? mod(y).ival : ival;
1368            int yval = y.ival;
1369    
1370            // Swap values so x > y.
1371            if (yval > xval)
1372              {
1373                int tmp = xval; xval = yval; yval = tmp;
1374                swapped = true;
1375              }
1376            // Normally, the result is in the 2nd element of the array, but
1377            // if originally x < y, then x and y were swapped and the result
1378            // is in the 1st element of the array.
1379            result.ival =
1380              euclidInv(yval, xval % yval, xval / yval)[swapped ? 0 : 1];
1381    
1382            // Result can't be negative, so make it positive by adding the
1383            // original modulus, y.ival (not the possibly "swapped" yval).
1384            if (result.ival < 0)
1385              result.ival += y.ival;
1386          }
1387        else
1388          {
1389            // As above, force this to be a positive value via modulo math.
1390            BigInteger x = isNegative() ? this.mod(y) : this;
1391    
1392            // Swap values so x > y.
1393            if (x.compareTo(y) < 0)
1394              {
1395                result = x; x = y; y = result; // use 'result' as a work var
1396                swapped = true;
1397              }
1398            // As above (for ints), result will be in the 2nd element unless
1399            // the original x and y were swapped.
1400            BigInteger rem = new BigInteger();
1401            BigInteger quot = new BigInteger();
1402            divide(x, y, quot, rem, FLOOR);
1403            // quot and rem may not be in canonical form. ensure
1404            rem.canonicalize();
1405            quot.canonicalize();
1406            BigInteger[] xy = new BigInteger[2];
1407            euclidInv(y, rem, quot, xy);
1408            result = swapped ? xy[0] : xy[1];
1409    
1410            // Result can't be negative, so make it positive by adding the
1411            // original modulus, y (which is now x if they were swapped).
1412            if (result.isNegative())
1413              result = add(result, swapped ? x : y, 1);
1414          }
1415        
1416        return result;
1417      }
1418    
1419      public BigInteger modPow(BigInteger exponent, BigInteger m)
1420      {
1421        if (USING_NATIVE)
1422          {
1423            int dummy = exponent.signum; // force NPE check
1424            if (m.mpz.compare(ZERO.mpz) < 1)
1425              throw new ArithmeticException("non-positive modulo");
1426    
1427            BigInteger result = new BigInteger();
1428            mpz.modPow(exponent.mpz, m.mpz, result.mpz);
1429            return result;
1430          }
1431    
1432        if (m.isNegative() || m.isZero())
1433          throw new ArithmeticException("non-positive modulo");
1434    
1435        if (exponent.isNegative())
1436          return modInverse(m).modPow(exponent.negate(), m);
1437        if (exponent.isOne())
1438          return mod(m);
1439    
1440        // To do this naively by first raising this to the power of exponent
1441        // and then performing modulo m would be extremely expensive, especially
1442        // for very large numbers.  The solution is found in Number Theory
1443        // where a combination of partial powers and moduli can be done easily.
1444        //
1445        // We'll use the algorithm for Additive Chaining which can be found on
1446        // p. 244 of "Applied Cryptography, Second Edition" by Bruce Schneier.
1447        BigInteger s = ONE;
1448        BigInteger t = this;
1449        BigInteger u = exponent;
1450    
1451        while (!u.isZero())
1452          {
1453            if (u.and(ONE).isOne())
1454              s = times(s, t).mod(m);
1455            u = u.shiftRight(1);
1456            t = times(t, t).mod(m);
1457          }
1458    
1459        return s;
1460      }
1461    
1462      /** Calculate Greatest Common Divisor for non-negative ints. */
1463      private static int gcd(int a, int b)
1464      {
1465        // Euclid's algorithm, copied from libg++.
1466        int tmp;
1467        if (b > a)
1468          {
1469            tmp = a; a = b; b = tmp;
1470          }
1471        for(;;)
1472          {
1473            if (b == 0)
1474              return a;
1475            if (b == 1)
1476              return b;
1477            tmp = b;
1478                b = a % b;
1479                a = tmp;
1480              }
1481          }
1482    
1483      public BigInteger gcd(BigInteger y)
1484      {
1485        if (USING_NATIVE)
1486          {
1487            int dummy = y.signum; // force NPE check
1488            BigInteger result = new BigInteger();
1489            mpz.gcd(y.mpz, result.mpz);
1490            return result;
1491          }
1492    
1493        int xval = ival;
1494        int yval = y.ival;
1495        if (words == null)
1496          {
1497            if (xval == 0)
1498              return abs(y);
1499            if (y.words == null
1500                && xval != Integer.MIN_VALUE && yval != Integer.MIN_VALUE)
1501              {
1502                if (xval < 0)
1503                  xval = -xval;
1504                if (yval < 0)
1505                  yval = -yval;
1506                return valueOf(gcd(xval, yval));
1507              }
1508            xval = 1;
1509          }
1510        if (y.words == null)
1511          {
1512            if (yval == 0)
1513              return abs(this);
1514            yval = 1;
1515          }
1516        int len = (xval > yval ? xval : yval) + 1;
1517        int[] xwords = new int[len];
1518        int[] ywords = new int[len];
1519        getAbsolute(xwords);
1520        y.getAbsolute(ywords);
1521        len = MPN.gcd(xwords, ywords, len);
1522        BigInteger result = new BigInteger(0);
1523        result.ival = len;
1524        result.words = xwords;
1525        return result.canonicalize();
1526      }
1527    
1528      /**
1529       * <p>Returns <code>true</code> if this BigInteger is probably prime,
1530       * <code>false</code> if it's definitely composite. If <code>certainty</code>
1531       * is <code><= 0</code>, <code>true</code> is returned.</p>
1532       *
1533       * @param certainty a measure of the uncertainty that the caller is willing
1534       * to tolerate: if the call returns <code>true</code> the probability that
1535       * this BigInteger is prime exceeds <code>(1 - 1/2<sup>certainty</sup>)</code>.
1536       * The execution time of this method is proportional to the value of this
1537       * parameter.
1538       * @return <code>true</code> if this BigInteger is probably prime,
1539       * <code>false</code> if it's definitely composite.
1540       */
1541      public boolean isProbablePrime(int certainty)
1542      {
1543        if (certainty < 1)
1544          return true;
1545    
1546        if (USING_NATIVE)
1547          return mpz.testPrimality(certainty) != 0;
1548    
1549        /** We'll use the Rabin-Miller algorithm for doing a probabilistic
1550         * primality test.  It is fast, easy and has faster decreasing odds of a
1551         * composite passing than with other tests.  This means that this
1552         * method will actually have a probability much greater than the
1553         * 1 - .5^certainty specified in the JCL (p. 117), but I don't think
1554         * anyone will complain about better performance with greater certainty.
1555         *
1556         * The Rabin-Miller algorithm can be found on pp. 259-261 of "Applied
1557         * Cryptography, Second Edition" by Bruce Schneier.
1558         */
1559    
1560        // First rule out small prime factors
1561        BigInteger rem = new BigInteger();
1562        int i;
1563        for (i = 0; i < primes.length; i++)
1564          {
1565            if (words == null && ival == primes[i])
1566              return true;
1567    
1568            divide(this, smallFixNums[primes[i] - minFixNum], null, rem, TRUNCATE);
1569            if (rem.canonicalize().isZero())
1570              return false;
1571          }
1572    
1573        // Now perform the Rabin-Miller test.
1574    
1575        // Set b to the number of times 2 evenly divides (this - 1).
1576        // I.e. 2^b is the largest power of 2 that divides (this - 1).
1577        BigInteger pMinus1 = add(this, -1);
1578        int b = pMinus1.getLowestSetBit();
1579    
1580        // Set m such that this = 1 + 2^b * m.
1581        BigInteger m = pMinus1.divide(valueOf(2L).pow(b));
1582    
1583        // The HAC (Handbook of Applied Cryptography), Alfred Menezes & al. Note
1584        // 4.49 (controlling the error probability) gives the number of trials
1585        // for an error probability of 1/2**80, given the number of bits in the
1586        // number to test.  we shall use these numbers as is if/when 'certainty'
1587        // is less or equal to 80, and twice as much if it's greater.
1588        int bits = this.bitLength();
1589        for (i = 0; i < k.length; i++)
1590          if (bits <= k[i])
1591            break;
1592        int trials = t[i];
1593        if (certainty > 80)
1594          trials *= 2;
1595        BigInteger z;
1596        for (int t = 0; t < trials; t++)
1597          {
1598            // The HAC (Handbook of Applied Cryptography), Alfred Menezes & al.
1599            // Remark 4.28 states: "...A strategy that is sometimes employed
1600            // is to fix the bases a to be the first few primes instead of
1601            // choosing them at random.
1602            z = smallFixNums[primes[t] - minFixNum].modPow(m, this);
1603            if (z.isOne() || z.equals(pMinus1))
1604              continue;                     // Passes the test; may be prime.
1605    
1606            for (i = 0; i < b; )
1607              {
1608                if (z.isOne())
1609                  return false;
1610                i++;
1611                if (z.equals(pMinus1))
1612                  break;                    // Passes the test; may be prime.
1613    
1614                z = z.modPow(valueOf(2), this);
1615              }
1616    
1617            if (i == b && !z.equals(pMinus1))
1618              return false;
1619          }
1620        return true;
1621      }
1622    
1623      private void setInvert()
1624      {
1625        if (words == null)
1626          ival = ~ival;
1627        else
1628          {
1629            for (int i = ival;  --i >= 0; )
1630              words[i] = ~words[i];
1631          }
1632      }
1633    
1634      private void setShiftLeft(BigInteger x, int count)
1635      {
1636        int[] xwords;
1637        int xlen;
1638        if (x.words == null)
1639          {
1640            if (count < 32)
1641              {
1642                set((long) x.ival << count);
1643                return;
1644              }
1645            xwords = new int[1];
1646            xwords[0] = x.ival;
1647            xlen = 1;
1648          }
1649        else
1650          {
1651            xwords = x.words;
1652            xlen = x.ival;
1653          }
1654        int word_count = count >> 5;
1655        count &= 31;
1656        int new_len = xlen + word_count;
1657        if (count == 0)
1658          {
1659            realloc(new_len);
1660            for (int i = xlen;  --i >= 0; )
1661              words[i+word_count] = xwords[i];
1662          }
1663        else
1664          {
1665            new_len++;
1666            realloc(new_len);
1667            int shift_out = MPN.lshift(words, word_count, xwords, xlen, count);
1668            count = 32 - count;
1669            words[new_len-1] = (shift_out << count) >> count;  // sign-extend.
1670          }
1671        ival = new_len;
1672        for (int i = word_count;  --i >= 0; )
1673          words[i] = 0;
1674      }
1675    
1676      private void setShiftRight(BigInteger x, int count)
1677      {
1678        if (x.words == null)
1679          set(count < 32 ? x.ival >> count : x.ival < 0 ? -1 : 0);
1680        else if (count == 0)
1681          set(x);
1682        else
1683          {
1684            boolean neg = x.isNegative();
1685            int word_count = count >> 5;
1686            count &= 31;
1687            int d_len = x.ival - word_count;
1688            if (d_len <= 0)
1689              set(neg ? -1 : 0);
1690            else
1691              {
1692                if (words == null || words.length < d_len)
1693                  realloc(d_len);
1694                MPN.rshift0 (words, x.words, word_count, d_len, count);
1695                ival = d_len;
1696                if (neg)
1697                  words[d_len-1] |= -2 << (31 - count);
1698              }
1699          }
1700      }
1701    
1702      private void setShift(BigInteger x, int count)
1703      {
1704        if (count > 0)
1705          setShiftLeft(x, count);
1706        else
1707          setShiftRight(x, -count);
1708      }
1709    
1710      private static BigInteger shift(BigInteger x, int count)
1711      {
1712        if (x.words == null)
1713          {
1714            if (count <= 0)
1715              return valueOf(count > -32 ? x.ival >> (-count) : x.ival < 0 ? -1 : 0);
1716            if (count < 32)
1717              return valueOf((long) x.ival << count);
1718          }
1719        if (count == 0)
1720          return x;
1721        BigInteger result = new BigInteger(0);
1722        result.setShift(x, count);
1723        return result.canonicalize();
1724      }
1725    
1726      public BigInteger shiftLeft(int n)
1727      {
1728        if (n == 0)
1729          return this;
1730    
1731        if (USING_NATIVE)
1732          {
1733            BigInteger result = new BigInteger();
1734            if (n < 0)
1735              mpz.shiftRight(-n, result.mpz);
1736            else
1737              mpz.shiftLeft(n, result.mpz);
1738            return result;
1739          }
1740    
1741        return shift(this, n);
1742      }
1743    
1744      public BigInteger shiftRight(int n)
1745      {
1746        if (n == 0)
1747          return this;
1748    
1749        if (USING_NATIVE)
1750          {
1751            BigInteger result = new BigInteger();
1752            if (n < 0)
1753              mpz.shiftLeft(-n, result.mpz);
1754            else
1755              mpz.shiftRight(n, result.mpz);
1756            return result;
1757          }
1758    
1759        return shift(this, -n);
1760      }
1761    
1762      private void format(int radix, CPStringBuilder buffer)
1763      {
1764        if (words == null)
1765          buffer.append(Integer.toString(ival, radix));
1766        else if (ival <= 2)
1767          buffer.append(Long.toString(longValue(), radix));
1768        else
1769          {
1770            boolean neg = isNegative();
1771            int[] work;
1772            if (neg || radix != 16)
1773              {
1774                work = new int[ival];
1775                getAbsolute(work);
1776              }
1777            else
1778              work = words;
1779            int len = ival;
1780    
1781            if (radix == 16)
1782              {
1783                if (neg)
1784                  buffer.append('-');
1785                int buf_start = buffer.length();
1786                for (int i = len;  --i >= 0; )
1787                  {
1788                    int word = work[i];
1789                    for (int j = 8;  --j >= 0; )
1790                      {
1791                        int hex_digit = (word >> (4 * j)) & 0xF;
1792                        // Suppress leading zeros:
1793                        if (hex_digit > 0 || buffer.length() > buf_start)
1794                          buffer.append(Character.forDigit(hex_digit, 16));
1795                      }
1796                  }
1797              }
1798            else
1799              {
1800                int i = buffer.length();
1801                for (;;)
1802                  {
1803                    int digit = MPN.divmod_1(work, work, len, radix);
1804                    buffer.append(Character.forDigit(digit, radix));
1805                    while (len > 0 && work[len-1] == 0) len--;
1806                    if (len == 0)
1807                      break;
1808                  }
1809                if (neg)
1810                  buffer.append('-');
1811                /* Reverse buffer. */
1812                int j = buffer.length() - 1;
1813                while (i < j)
1814                  {
1815                    char tmp = buffer.charAt(i);
1816                    buffer.setCharAt(i, buffer.charAt(j));
1817                    buffer.setCharAt(j, tmp);
1818                    i++;  j--;
1819                  }
1820              }
1821          }
1822      }
1823    
1824      public String toString()
1825      {
1826        return toString(10);
1827      }
1828    
1829      public String toString(int radix)
1830      {
1831        if (USING_NATIVE)
1832          return mpz.toString(radix);
1833    
1834        if (words == null)
1835          return Integer.toString(ival, radix);
1836        if (ival <= 2)
1837          return Long.toString(longValue(), radix);
1838        int buf_size = ival * (MPN.chars_per_word(radix) + 1);
1839        CPStringBuilder buffer = new CPStringBuilder(buf_size);
1840        format(radix, buffer);
1841        return buffer.toString();
1842      }
1843    
1844      public int intValue()
1845      {
1846        if (USING_NATIVE)
1847          {
1848            int result = mpz.absIntValue();
1849            return mpz.compare(ZERO.mpz) < 0 ? - result : result;
1850          }
1851    
1852        if (words == null)
1853          return ival;
1854        return words[0];
1855      }
1856    
1857      public long longValue()
1858      {
1859        if (USING_NATIVE)
1860          {
1861            long result;
1862            result = (abs().shiftRight(32)).mpz.absIntValue();
1863            result <<= 32;
1864            result |= mpz.absIntValue() & 0xFFFFFFFFL;
1865            return this.compareTo(ZERO) < 0 ? - result : result;
1866          }
1867    
1868        if (words == null)
1869          return ival;
1870        if (ival == 1)
1871          return words[0];
1872        return ((long)words[1] << 32) + ((long)words[0] & 0xffffffffL);
1873      }
1874    
1875      public int hashCode()
1876      {
1877        // FIXME: May not match hashcode of JDK.
1878        if (USING_NATIVE)
1879          {
1880            // TODO: profile to decide whether to make it native
1881            byte[] bytes = this.toByteArray();
1882            int result = 0;
1883            for (int i = 0; i < bytes.length; i++)
1884              result ^= (bytes[i] & 0xFF) << (8 * (i % 4));
1885            return result;
1886          }
1887    
1888        return words == null ? ival : (words[0] + words[ival - 1]);
1889      }
1890    
1891      /* Assumes x and y are both canonicalized. */
1892      private static boolean equals(BigInteger x, BigInteger y)
1893      {
1894        if (USING_NATIVE)
1895          return x.mpz.compare(y.mpz) == 0;
1896    
1897        if (x.words == null && y.words == null)
1898          return x.ival == y.ival;
1899        if (x.words == null || y.words == null || x.ival != y.ival)
1900          return false;
1901        for (int i = x.ival; --i >= 0; )
1902          {
1903            if (x.words[i] != y.words[i])
1904              return false;
1905          }
1906        return true;
1907      }
1908    
1909      /* Assumes this and obj are both canonicalized. */
1910      public boolean equals(Object obj)
1911      {
1912        if (! (obj instanceof BigInteger))
1913          return false;
1914        return equals(this, (BigInteger) obj);
1915      }
1916    
1917      private static BigInteger valueOf(byte[] digits, int byte_len,
1918                                        boolean negative, int radix)
1919      {
1920        int chars_per_word = MPN.chars_per_word(radix);
1921        int[] words = new int[byte_len / chars_per_word + 1];
1922        int size = MPN.set_str(words, digits, byte_len, radix);
1923        if (size == 0)
1924          return ZERO;
1925        if (words[size-1] < 0)
1926          words[size++] = 0;
1927        if (negative)
1928          negate(words, words, size);
1929        return make(words, size);
1930      }
1931    
1932      public double doubleValue()
1933      {
1934        if (USING_NATIVE)
1935          return mpz.doubleValue();
1936    
1937        if (words == null)
1938          return (double) ival;
1939        if (ival <= 2)
1940          return (double) longValue();
1941        if (isNegative())
1942          return neg(this).roundToDouble(0, true, false);
1943          return roundToDouble(0, false, false);
1944      }
1945    
1946      public float floatValue()
1947      {
1948        return (float) doubleValue();
1949      }
1950    
1951      /** Return true if any of the lowest n bits are one.
1952       * (false if n is negative).  */
1953      private boolean checkBits(int n)
1954      {
1955        if (n <= 0)
1956          return false;
1957        if (words == null)
1958          return n > 31 || ((ival & ((1 << n) - 1)) != 0);
1959        int i;
1960        for (i = 0; i < (n >> 5) ; i++)
1961          if (words[i] != 0)
1962            return true;
1963        return (n & 31) != 0 && (words[i] & ((1 << (n & 31)) - 1)) != 0;
1964      }
1965    
1966      /** Convert a semi-processed BigInteger to double.
1967       * Number must be non-negative.  Multiplies by a power of two, applies sign,
1968       * and converts to double, with the usual java rounding.
1969       * @param exp power of two, positive or negative, by which to multiply
1970       * @param neg true if negative
1971       * @param remainder true if the BigInteger is the result of a truncating
1972       * division that had non-zero remainder.  To ensure proper rounding in
1973       * this case, the BigInteger must have at least 54 bits.  */
1974      private double roundToDouble(int exp, boolean neg, boolean remainder)
1975      {
1976        // Compute length.
1977        int il = bitLength();
1978    
1979        // Exponent when normalized to have decimal point directly after
1980        // leading one.  This is stored excess 1023 in the exponent bit field.
1981        exp += il - 1;
1982    
1983        // Gross underflow.  If exp == -1075, we let the rounding
1984        // computation determine whether it is minval or 0 (which are just
1985        // 0x0000 0000 0000 0001 and 0x0000 0000 0000 0000 as bit
1986        // patterns).
1987        if (exp < -1075)
1988          return neg ? -0.0 : 0.0;
1989    
1990        // gross overflow
1991        if (exp > 1023)
1992          return neg ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1993    
1994        // number of bits in mantissa, including the leading one.
1995        // 53 unless it's denormalized
1996        int ml = (exp >= -1022 ? 53 : 53 + exp + 1022);
1997    
1998        // Get top ml + 1 bits.  The extra one is for rounding.
1999        long m;
2000        int excess_bits = il - (ml + 1);
2001        if (excess_bits > 0)
2002          m = ((words == null) ? ival >> excess_bits
2003               : MPN.rshift_long(words, ival, excess_bits));
2004        else
2005          m = longValue() << (- excess_bits);
2006    
2007        // Special rounding for maxval.  If the number exceeds maxval by
2008        // any amount, even if it's less than half a step, it overflows.
2009        if (exp == 1023 && ((m >> 1) == (1L << 53) - 1))
2010          {
2011            if (remainder || checkBits(il - ml))
2012              return neg ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
2013            else
2014              return neg ? - Double.MAX_VALUE : Double.MAX_VALUE;
2015          }
2016    
2017        // Normal round-to-even rule: round up if the bit dropped is a one, and
2018        // the bit above it or any of the bits below it is a one.
2019        if ((m & 1) == 1
2020            && ((m & 2) == 2 || remainder || checkBits(excess_bits)))
2021          {
2022            m += 2;
2023            // Check if we overflowed the mantissa
2024            if ((m & (1L << 54)) != 0)
2025              {
2026                exp++;
2027                // renormalize
2028                m >>= 1;
2029              }
2030            // Check if a denormalized mantissa was just rounded up to a
2031            // normalized one.
2032            else if (ml == 52 && (m & (1L << 53)) != 0)
2033              exp++;
2034          }
2035            
2036        // Discard the rounding bit
2037        m >>= 1;
2038    
2039        long bits_sign = neg ? (1L << 63) : 0;
2040        exp += 1023;
2041        long bits_exp = (exp <= 0) ? 0 : ((long)exp) << 52;
2042        long bits_mant = m & ~(1L << 52);
2043        return Double.longBitsToDouble(bits_sign | bits_exp | bits_mant);
2044      }
2045    
2046      /** Copy the abolute value of this into an array of words.
2047       * Assumes words.length >= (this.words == null ? 1 : this.ival).
2048       * Result is zero-extended, but need not be a valid 2's complement number.
2049       */
2050      private void getAbsolute(int[] words)
2051      {
2052        int len;
2053        if (this.words == null)
2054          {
2055            len = 1;
2056            words[0] = this.ival;
2057          }
2058        else
2059          {
2060            len = this.ival;
2061            for (int i = len;  --i >= 0; )
2062              words[i] = this.words[i];
2063          }
2064        if (words[len - 1] < 0)
2065          negate(words, words, len);
2066        for (int i = words.length;  --i > len; )
2067          words[i] = 0;
2068      }
2069    
2070      /** Set dest[0:len-1] to the negation of src[0:len-1].
2071       * Return true if overflow (i.e. if src is -2**(32*len-1)).
2072       * Ok for src==dest. */
2073      private static boolean negate(int[] dest, int[] src, int len)
2074      {
2075        long carry = 1;
2076        boolean negative = src[len-1] < 0;
2077        for (int i = 0;  i < len;  i++)
2078          {
2079            carry += ((long) (~src[i]) & 0xffffffffL);
2080            dest[i] = (int) carry;
2081            carry >>= 32;
2082          }
2083        return (negative && dest[len-1] < 0);
2084      }
2085    
2086      /** Destructively set this to the negative of x.
2087       * It is OK if x==this.*/
2088      private void setNegative(BigInteger x)
2089      {
2090        int len = x.ival;
2091        if (x.words == null)
2092          {
2093            if (len == Integer.MIN_VALUE)
2094              set(- (long) len);
2095            else
2096              set(-len);
2097            return;
2098          }
2099        realloc(len + 1);
2100        if (negate(words, x.words, len))
2101          words[len++] = 0;
2102        ival = len;
2103      }
2104    
2105      /** Destructively negate this. */
2106      private void setNegative()
2107      {
2108        setNegative(this);
2109      }
2110    
2111      private static BigInteger abs(BigInteger x)
2112      {
2113        return x.isNegative() ? neg(x) : x;
2114      }
2115    
2116      public BigInteger abs()
2117      {
2118        if (USING_NATIVE)
2119          {
2120            BigInteger result = new BigInteger();
2121            mpz.abs(result.mpz);
2122            return result;
2123          }
2124    
2125        return abs(this);
2126      }
2127    
2128      private static BigInteger neg(BigInteger x)
2129      {
2130        if (x.words == null && x.ival != Integer.MIN_VALUE)
2131          return valueOf(- x.ival);
2132        BigInteger result = new BigInteger(0);
2133        result.setNegative(x);
2134        return result.canonicalize();
2135      }
2136    
2137      public BigInteger negate()
2138      {
2139        if (USING_NATIVE)
2140          {
2141            BigInteger result = new BigInteger();
2142            mpz.negate(result.mpz);
2143            return result;
2144          }
2145    
2146        return neg(this);
2147      }
2148    
2149      /** Calculates ceiling(log2(this < 0 ? -this : this+1))
2150       * See Common Lisp: the Language, 2nd ed, p. 361.
2151       */
2152      public int bitLength()
2153      {
2154        if (USING_NATIVE)
2155          return mpz.bitLength();
2156    
2157        if (words == null)
2158          return MPN.intLength(ival);
2159          return MPN.intLength(words, ival);
2160      }
2161    
2162      public byte[] toByteArray()
2163      {
2164        if (signum() == 0)
2165          return new byte[1];
2166    
2167        if (USING_NATIVE)
2168          {
2169            // the minimal number of bytes required to represent the MPI is function
2170            // of (a) its bit-length, and (b) its sign.  only when this MPI is both
2171            // positive, and its bit-length is a multiple of 8 do we add one zero
2172            // bit for its sign.  we do this so if we construct a new MPI from the
2173            // resulting byte array, we wouldn't mistake a positive number, whose
2174            // bit-length is a multiple of 8, for a similar-length negative one.
2175            int bits = bitLength();
2176            if (bits % 8 == 0 || this.signum() == 1)
2177              bits++;
2178            byte[] bytes = new byte[(bits + 7) / 8];
2179            mpz.toByteArray(bytes);
2180            return bytes;
2181          }
2182    
2183        // Determine number of bytes needed.  The method bitlength returns
2184        // the size without the sign bit, so add one bit for that and then
2185        // add 7 more to emulate the ceil function using integer math.
2186        byte[] bytes = new byte[(bitLength() + 1 + 7) / 8];
2187        int nbytes = bytes.length;
2188    
2189        int wptr = 0;
2190        int word;
2191    
2192        // Deal with words array until one word or less is left to process.
2193        // If BigInteger is an int, then it is in ival and nbytes will be <= 4.
2194        while (nbytes > 4)
2195          {
2196            word = words[wptr++];
2197            for (int i = 4; i > 0; --i, word >>= 8)
2198              bytes[--nbytes] = (byte) word;
2199          }
2200    
2201        // Deal with the last few bytes.  If BigInteger is an int, use ival.
2202        word = (words == null) ? ival : words[wptr];
2203        for ( ; nbytes > 0; word >>= 8)
2204          bytes[--nbytes] = (byte) word;
2205    
2206        return bytes;
2207      }
2208    
2209      /** Return the boolean opcode (for bitOp) for swapped operands.
2210       * I.e. bitOp(swappedOp(op), x, y) == bitOp(op, y, x).
2211       */
2212      private static int swappedOp(int op)
2213      {
2214        return
2215        "\000\001\004\005\002\003\006\007\010\011\014\015\012\013\016\017"
2216        .charAt(op);
2217      }
2218    
2219      /** Do one the the 16 possible bit-wise operations of two BigIntegers. */
2220      private static BigInteger bitOp(int op, BigInteger x, BigInteger y)
2221      {
2222        switch (op)
2223          {
2224            case 0:  return ZERO;
2225            case 1:  return x.and(y);
2226            case 3:  return x;
2227            case 5:  return y;
2228            case 15: return valueOf(-1);
2229          }
2230        BigInteger result = new BigInteger();
2231        setBitOp(result, op, x, y);
2232        return result.canonicalize();
2233      }
2234    
2235      /** Do one the the 16 possible bit-wise operations of two BigIntegers. */
2236      private static void setBitOp(BigInteger result, int op,
2237                                   BigInteger x, BigInteger y)
2238      {
2239        if ((y.words != null) && (x.words == null || x.ival < y.ival))
2240          {
2241            BigInteger temp = x;  x = y;  y = temp;
2242            op = swappedOp(op);
2243          }
2244        int xi;
2245        int yi;
2246        int xlen, ylen;
2247        if (y.words == null)
2248          {
2249            yi = y.ival;
2250            ylen = 1;
2251          }
2252        else
2253          {
2254            yi = y.words[0];
2255            ylen = y.ival;
2256          }
2257        if (x.words == null)
2258          {
2259            xi = x.ival;
2260            xlen = 1;
2261          }
2262        else
2263          {
2264            xi = x.words[0];
2265            xlen = x.ival;
2266          }
2267        if (xlen > 1)
2268          result.realloc(xlen);
2269        int[] w = result.words;
2270        int i = 0;
2271        // Code for how to handle the remainder of x.
2272        // 0:  Truncate to length of y.
2273        // 1:  Copy rest of x.
2274        // 2:  Invert rest of x.
2275        int finish = 0;
2276        int ni;
2277        switch (op)
2278          {
2279          case 0:  // clr
2280            ni = 0;
2281            break;
2282          case 1: // and
2283            for (;;)
2284              {
2285                ni = xi & yi;
2286                if (i+1 >= ylen) break;
2287                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2288              }
2289            if (yi < 0) finish = 1;
2290            break;
2291          case 2: // andc2
2292            for (;;)
2293              {
2294                ni = xi & ~yi;
2295                if (i+1 >= ylen) break;
2296                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2297              }
2298            if (yi >= 0) finish = 1;
2299            break;
2300          case 3:  // copy x
2301            ni = xi;
2302            finish = 1;  // Copy rest
2303            break;
2304          case 4: // andc1
2305            for (;;)
2306              {
2307                ni = ~xi & yi;
2308                if (i+1 >= ylen) break;
2309                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2310              }
2311            if (yi < 0) finish = 2;
2312            break;
2313          case 5: // copy y
2314            for (;;)
2315              {
2316                ni = yi;
2317                if (i+1 >= ylen) break;
2318                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2319              }
2320            break;
2321          case 6:  // xor
2322            for (;;)
2323              {
2324                ni = xi ^ yi;
2325                if (i+1 >= ylen) break;
2326                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2327              }
2328            finish = yi < 0 ? 2 : 1;
2329            break;
2330          case 7:  // ior
2331            for (;;)
2332              {
2333                ni = xi | yi;
2334                if (i+1 >= ylen) break;
2335                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2336              }
2337            if (yi >= 0) finish = 1;
2338            break;
2339          case 8:  // nor
2340            for (;;)
2341              {
2342                ni = ~(xi | yi);
2343                if (i+1 >= ylen) break;
2344                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2345              }
2346            if (yi >= 0)  finish = 2;
2347            break;
2348          case 9:  // eqv [exclusive nor]
2349            for (;;)
2350              {
2351                ni = ~(xi ^ yi);
2352                if (i+1 >= ylen) break;
2353                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2354              }
2355            finish = yi >= 0 ? 2 : 1;
2356            break;
2357          case 10:  // c2
2358            for (;;)
2359              {
2360                ni = ~yi;
2361                if (i+1 >= ylen) break;
2362                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2363              }
2364            break;
2365          case 11:  // orc2
2366            for (;;)
2367              {
2368                ni = xi | ~yi;
2369                if (i+1 >= ylen) break;
2370                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2371              }
2372            if (yi < 0)  finish = 1;
2373            break;
2374          case 12:  // c1
2375            ni = ~xi;
2376            finish = 2;
2377            break;
2378          case 13:  // orc1
2379            for (;;)
2380              {
2381                ni = ~xi | yi;
2382                if (i+1 >= ylen) break;
2383                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2384              }
2385            if (yi >= 0) finish = 2;
2386            break;
2387          case 14:  // nand
2388            for (;;)
2389              {
2390                ni = ~(xi & yi);
2391                if (i+1 >= ylen) break;
2392                w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
2393              }
2394            if (yi < 0) finish = 2;
2395            break;
2396          default:
2397          case 15:  // set
2398            ni = -1;
2399            break;
2400          }
2401        // Here i==ylen-1; w[0]..w[i-1] have the correct result;
2402        // and ni contains the correct result for w[i+1].
2403        if (i+1 == xlen)
2404          finish = 0;
2405        switch (finish)
2406          {
2407          case 0:
2408            if (i == 0 && w == null)
2409              {
2410                result.ival = ni;
2411                return;
2412              }
2413            w[i++] = ni;
2414            break;
2415          case 1:  w[i] = ni;  while (++i < xlen)  w[i] = x.words[i];  break;
2416          case 2:  w[i] = ni;  while (++i < xlen)  w[i] = ~x.words[i];  break;
2417          }
2418        result.ival = i;
2419      }
2420    
2421      /** Return the logical (bit-wise) "and" of a BigInteger and an int. */
2422      private static BigInteger and(BigInteger x, int y)
2423      {
2424        if (x.words == null)
2425          return valueOf(x.ival & y);
2426        if (y >= 0)
2427          return valueOf(x.words[0] & y);
2428        int len = x.ival;
2429        int[] words = new int[len];
2430        words[0] = x.words[0] & y;
2431        while (--len > 0)
2432          words[len] = x.words[len];
2433        return make(words, x.ival);
2434      }
2435    
2436      /** Return the logical (bit-wise) "and" of two BigIntegers. */
2437      public BigInteger and(BigInteger y)
2438      {
2439        if (USING_NATIVE)
2440          {
2441            int dummy = y.signum; // force NPE check
2442            BigInteger result = new BigInteger();
2443            mpz.and(y.mpz, result.mpz);
2444            return result;
2445          }
2446    
2447        if (y.words == null)
2448          return and(this, y.ival);
2449        else if (words == null)
2450          return and(y, ival);
2451    
2452        BigInteger x = this;
2453        if (ival < y.ival)
2454          {
2455            BigInteger temp = this;  x = y;  y = temp;
2456          }
2457        int i;
2458        int len = y.isNegative() ? x.ival : y.ival;
2459        int[] words = new int[len];
2460        for (i = 0;  i < y.ival;  i++)
2461          words[i] = x.words[i] & y.words[i];
2462        for ( ; i < len;  i++)
2463          words[i] = x.words[i];
2464        return make(words, len);
2465      }
2466    
2467      /** Return the logical (bit-wise) "(inclusive) or" of two BigIntegers. */
2468      public BigInteger or(BigInteger y)
2469      {
2470        if (USING_NATIVE)
2471          {
2472            int dummy = y.signum; // force NPE check
2473            BigInteger result = new BigInteger();
2474            mpz.or(y.mpz, result.mpz);
2475            return result;
2476          }
2477    
2478        return bitOp(7, this, y);
2479      }
2480    
2481      /** Return the logical (bit-wise) "exclusive or" of two BigIntegers. */
2482      public BigInteger xor(BigInteger y)
2483      {
2484        if (USING_NATIVE)
2485          {
2486            int dummy = y.signum; // force NPE check
2487            BigInteger result = new BigInteger();
2488            mpz.xor(y.mpz, result.mpz);
2489            return result;
2490          }
2491    
2492        return bitOp(6, this, y);
2493      }
2494    
2495      /** Return the logical (bit-wise) negation of a BigInteger. */
2496      public BigInteger not()
2497      {
2498        if (USING_NATIVE)
2499          {
2500            BigInteger result = new BigInteger();
2501            mpz.not(result.mpz);
2502            return result;
2503          }
2504    
2505        return bitOp(12, this, ZERO);
2506      }
2507    
2508      public BigInteger andNot(BigInteger val)
2509      {
2510        if (USING_NATIVE)
2511          {
2512            int dummy = val.signum; // force NPE check
2513            BigInteger result = new BigInteger();
2514            mpz.andNot(val.mpz, result.mpz);
2515            return result;
2516          }
2517    
2518        return and(val.not());
2519      }
2520    
2521      public BigInteger clearBit(int n)
2522      {
2523        if (n < 0)
2524          throw new ArithmeticException();
2525    
2526        if (USING_NATIVE)
2527          {
2528            BigInteger result = new BigInteger();
2529            mpz.setBit(n, false, result.mpz);
2530            return result;
2531          }
2532    
2533        return and(ONE.shiftLeft(n).not());
2534      }
2535    
2536      public BigInteger setBit(int n)
2537      {
2538        if (n < 0)
2539          throw new ArithmeticException();
2540    
2541        if (USING_NATIVE)
2542          {
2543            BigInteger result = new BigInteger();
2544            mpz.setBit(n, true, result.mpz);
2545            return result;
2546          }
2547    
2548        return or(ONE.shiftLeft(n));
2549      }
2550    
2551      public boolean testBit(int n)
2552      {
2553        if (n < 0)
2554          throw new ArithmeticException();
2555    
2556        if (USING_NATIVE)
2557          return mpz.testBit(n) != 0;
2558    
2559        return !and(ONE.shiftLeft(n)).isZero();
2560      }
2561    
2562      public BigInteger flipBit(int n)
2563      {
2564        if (n < 0)
2565          throw new ArithmeticException();
2566    
2567        if (USING_NATIVE)
2568          {
2569            BigInteger result = new BigInteger();
2570            mpz.flipBit(n, result.mpz);
2571            return result;
2572          }
2573    
2574        return xor(ONE.shiftLeft(n));
2575      }
2576    
2577      public int getLowestSetBit()
2578      {
2579        if (USING_NATIVE)
2580          return mpz.compare(ZERO.mpz) == 0 ? -1 : mpz.lowestSetBit();
2581    
2582        if (isZero())
2583          return -1;
2584    
2585        if (words == null)
2586          return MPN.findLowestBit(ival);
2587        else
2588          return MPN.findLowestBit(words);
2589      }
2590    
2591      // bit4count[I] is number of '1' bits in I.
2592      private static final byte[] bit4_count = { 0, 1, 1, 2,  1, 2, 2, 3,
2593                                                 1, 2, 2, 3,  2, 3, 3, 4};
2594    
2595      private static int bitCount(int i)
2596      {
2597        int count = 0;
2598        while (i != 0)
2599          {
2600            count += bit4_count[i & 15];
2601            i >>>= 4;
2602          }
2603        return count;
2604      }
2605    
2606      private static int bitCount(int[] x, int len)
2607      {
2608        int count = 0;
2609        while (--len >= 0)
2610          count += bitCount(x[len]);
2611        return count;
2612      }
2613    
2614      /** Count one bits in a BigInteger.
2615       * If argument is negative, count zero bits instead. */
2616      public int bitCount()
2617      {
2618        if (USING_NATIVE)
2619          return mpz.bitCount();
2620    
2621        int i, x_len;
2622        int[] x_words = words;
2623        if (x_words == null)
2624          {
2625            x_len = 1;
2626            i = bitCount(ival);
2627          }
2628        else
2629          {
2630            x_len = ival;
2631            i = bitCount(x_words, x_len);
2632          }
2633        return isNegative() ? x_len * 32 - i : i;
2634      }
2635    
2636      private void readObject(ObjectInputStream s)
2637        throws IOException, ClassNotFoundException
2638      {
2639        if (USING_NATIVE)
2640          {
2641            mpz = new GMP();
2642            s.defaultReadObject();
2643            if (signum != 0)
2644              mpz.fromByteArray(magnitude);
2645            // else it's zero and we need to do nothing
2646          }
2647        else
2648          {
2649            s.defaultReadObject();
2650            if (magnitude.length == 0 || signum == 0)
2651              {
2652                this.ival = 0;
2653                this.words = null;
2654              }
2655            else
2656              {
2657                words = byteArrayToIntArray(magnitude, signum < 0 ? -1 : 0);
2658                BigInteger result = make(words, words.length);
2659                this.ival = result.ival;
2660                this.words = result.words;        
2661              }    
2662          }
2663      }
2664    
2665      private void writeObject(ObjectOutputStream s)
2666        throws IOException, ClassNotFoundException
2667      {
2668        signum = signum();
2669        magnitude = signum == 0 ? new byte[0] : toByteArray();
2670        s.defaultWriteObject();
2671        magnitude = null; // not needed anymore
2672      }
2673    
2674      // inner class(es) ..........................................................
2675    
2676    }