libjava/java/math/BigInteger.java

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