001/* java.math.BigInteger -- Arbitary precision integers
002   Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007  Free Software Foundation, Inc.
003
004This file is part of GNU Classpath.
005
006GNU Classpath is free software; you can redistribute it and/or modify
007it under the terms of the GNU General Public License as published by
008the Free Software Foundation; either version 2, or (at your option)
009any later version.
010
011GNU Classpath is distributed in the hope that it will be useful, but
012WITHOUT ANY WARRANTY; without even the implied warranty of
013MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
014General Public License for more details.
015
016You should have received a copy of the GNU General Public License
017along with GNU Classpath; see the file COPYING.  If not, write to the
018Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
01902110-1301 USA.
020
021Linking this library statically or dynamically with other modules is
022making a combined work based on this library.  Thus, the terms and
023conditions of the GNU General Public License cover the whole
024combination.
025
026As a special exception, the copyright holders of this library give you
027permission to link this library with independent modules to produce an
028executable, regardless of the license terms of these independent
029modules, and to copy and distribute the resulting executable under
030terms of your choice, provided that you also meet, for each linked
031independent module, the terms and conditions of the license of that
032module.  An independent module is a module which is not derived from
033or based on this library.  If you modify this library, you may extend
034this exception to your version of the library, but you are not
035obligated to do so.  If you do not wish to do so, delete this
036exception statement from your version. */
037
038
039package java.math;
040
041import gnu.classpath.Configuration;
042
043import gnu.java.lang.CPStringBuilder;
044import gnu.java.math.GMP;
045import gnu.java.math.MPN;
046
047import java.io.IOException;
048import java.io.ObjectInputStream;
049import java.io.ObjectOutputStream;
050import java.util.Random;
051import 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 */
065public 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}