libjava/java/math/BigInteger.java

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