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