Makefile.am: Added MPN.java and BigInteger.java.

	* Makefile.am: Added MPN.java and BigInteger.java.
	* Makefile.in: Rebuilt.
	* gnu/gcj/math/MPN.java: New file.
	* java/math/BigInteger.java: New file.

From-SVN: r31794

libjava/ChangeLog

+2000-02-04  Warren Levy  <warrenl@cygnus.com>
+
+	* Makefile.am: Added MPN.java and BigInteger.java.
+	* Makefile.in: Rebuilt.
+	* gnu/gcj/math/MPN.java: New file.
+	* java/math/BigInteger.java: New file.
+
 2000-02-04  Tom Tromey  <tromey@cygnus.com>
 
 	* defineclass.cc (handleMethodsBegin): Allocate _Jv_MethodBase

libjava/Makefile.am

@@ -521,6 +521,7 @@ gnu/gcj/text/LocaleData_en.java	\
 gnu/gcj/text/LocaleData_en_US.java \
 gnu/gcj/text/SentenceBreakIterator.java	\
 gnu/gcj/text/WordBreakIterator.java \
+gnu/gcj/math/MPN.java \
 gnu/gcj/protocol/file/Connection.java \
 gnu/gcj/protocol/file/Handler.java \
 gnu/gcj/protocol/http/Connection.java \
@@ -661,6 +662,7 @@ java/lang/reflect/InvocationTargetException.java \
 java/lang/reflect/Member.java \
 java/lang/reflect/Method.java \
 java/lang/reflect/Modifier.java	\
+java/math/BigInteger.java \
 java/net/BindException.java \
 java/net/ConnectException.java \
 java/net/ContentHandler.java \

libjava/Makefile.in

@@ -335,6 +335,7 @@ gnu/gcj/text/LocaleData_en.java	\
 gnu/gcj/text/LocaleData_en_US.java \
 gnu/gcj/text/SentenceBreakIterator.java	\
 gnu/gcj/text/WordBreakIterator.java \
+gnu/gcj/math/MPN.java \
 gnu/gcj/protocol/file/Connection.java \
 gnu/gcj/protocol/file/Handler.java \
 gnu/gcj/protocol/http/Connection.java \
@@ -475,6 +476,7 @@ java/lang/reflect/InvocationTargetException.java \
 java/lang/reflect/Member.java \
 java/lang/reflect/Method.java \
 java/lang/reflect/Modifier.java	\
+java/math/BigInteger.java \
 java/net/BindException.java \
 java/net/ConnectException.java \
 java/net/ContentHandler.java \
@@ -750,7 +752,7 @@ DEP_FILES =  .deps/$(srcdir)/$(CONVERT_DIR)/gen-from-JIS.P \
 .deps/gnu/gcj/convert/Output_JavaSrc.P \
 .deps/gnu/gcj/convert/Output_SJIS.P .deps/gnu/gcj/convert/Output_UTF8.P \
 .deps/gnu/gcj/convert/Output_iconv.P \
-.deps/gnu/gcj/convert/UnicodeToBytes.P \
+.deps/gnu/gcj/convert/UnicodeToBytes.P .deps/gnu/gcj/math/MPN.P \
 .deps/gnu/gcj/protocol/file/Connection.P \
 .deps/gnu/gcj/protocol/file/Handler.P \
 .deps/gnu/gcj/protocol/http/Connection.P \
@@ -869,12 +871,12 @@ DEP_FILES =  .deps/$(srcdir)/$(CONVERT_DIR)/gen-from-JIS.P \
 .deps/java/lang/w_exp.P .deps/java/lang/w_fmod.P \
 .deps/java/lang/w_log.P .deps/java/lang/w_pow.P \
 .deps/java/lang/w_remainder.P .deps/java/lang/w_sqrt.P \
-.deps/java/net/BindException.P .deps/java/net/ConnectException.P \
-.deps/java/net/ContentHandler.P .deps/java/net/ContentHandlerFactory.P \
-.deps/java/net/DatagramPacket.P .deps/java/net/DatagramSocket.P \
-.deps/java/net/DatagramSocketImpl.P .deps/java/net/FileNameMap.P \
-.deps/java/net/HttpURLConnection.P .deps/java/net/InetAddress.P \
-.deps/java/net/JarURLConnection.P \
+.deps/java/math/BigInteger.P .deps/java/net/BindException.P \
+.deps/java/net/ConnectException.P .deps/java/net/ContentHandler.P \
+.deps/java/net/ContentHandlerFactory.P .deps/java/net/DatagramPacket.P \
+.deps/java/net/DatagramSocket.P .deps/java/net/DatagramSocketImpl.P \
+.deps/java/net/FileNameMap.P .deps/java/net/HttpURLConnection.P \
+.deps/java/net/InetAddress.P .deps/java/net/JarURLConnection.P \
 .deps/java/net/MalformedURLException.P .deps/java/net/MulticastSocket.P \
 .deps/java/net/NoRouteToHostException.P \
 .deps/java/net/PlainDatagramSocketImpl.P \

libjava/gnu/gcj/math/MPN.java

+/* Copyright (C) 1999, 2000  Red Hat, Inc.
+
+   This file is part of libgcj.
+
+This software is copyrighted work licensed under the terms of the
+Libgcj License.  Please consult the file "LIBGCJ_LICENSE" for
+details.  */
+
+// Included from Kawa 1.6.62 with permission of the author,
+// Per Bothner <per@bothner.com>.
+
+package gnu.gcj.math;
+
+/** This contains various low-level routines for unsigned bigints.
+ * The interfaces match the mpn interfaces in gmp,
+ * so it should be easy to replace them with fast native functions
+ * that are trivial wrappers around the mpn_ functions in gmp
+ * (at least on platforms that use 32-bit "limbs").
+ */
+
+public class MPN
+{
+  /** Add x[0:size-1] and y, and write the size least
+   * significant words of the result to dest.
+   * Return carry, either 0 or 1.
+   * All values are unsigned.
+   * This is basically the same as gmp's mpn_add_1. */
+  public static int add_1 (int[] dest, int[] x, int size, int y)
+  {
+    long carry = (long) y & 0xffffffffL;
+    for (int i = 0;  i < size;  i++)
+      {
+	carry += ((long) x[i] & 0xffffffffL);
+	dest[i] = (int) carry;
+	carry >>= 32;
+      }
+    return (int) carry;
+  }
+
+  /** Add x[0:len-1] and y[0:len-1] and write the len least
+   * significant words of the result to dest[0:len-1].
+   * All words are treated as unsigned.
+   * @return the carry, either 0 or 1
+   * This function is basically the same as gmp's mpn_add_n.
+   */
+  public static int add_n (int dest[], int[] x, int[] y, int len)
+  {
+    long carry = 0;
+    for (int i = 0; i < len;  i++)
+      {
+	carry += ((long) x[i] & 0xffffffffL)
+	  + ((long) y[i] & 0xffffffffL);
+	dest[i] = (int) carry;
+	carry >>>= 32;
+      }
+    return (int) carry;
+  }
+
+  /** Subtract Y[0:size-1] from X[0:size-1], and write
+   * the size least significant words of the result to dest[0:size-1].
+   * Return borrow, either 0 or 1.
+   * This is basically the same as gmp's mpn_sub_n function.
+   */
+
+  public static int sub_n (int[] dest, int[] X, int[] Y, int size)
+  {
+    int cy = 0;
+    for (int i = 0;  i < size;  i++)
+      {
+	int y = Y[i];
+	int x = X[i];
+	y += cy;	/* add previous carry to subtrahend */
+	// Invert the high-order bit, because: (unsigned) X > (unsigned) Y
+	// iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
+	cy = (y^0x80000000) < (cy^0x80000000) ? 1 : 0;
+	y = x - y;
+	cy += (y^0x80000000) > (x ^ 0x80000000) ? 1 : 0;
+	dest[i] = y;
+      }
+    return cy;
+  }
+
+  /** Multiply x[0:len-1] by y, and write the len least
+   * significant words of the product to dest[0:len-1].
+   * Return the most significant word of the product.
+   * All values are treated as if they were unsigned
+   * (i.e. masked with 0xffffffffL).
+   * OK if dest==x (not sure if this is guaranteed for mpn_mul_1).
+   * This function is basically the same as gmp's mpn_mul_1.
+   */
+
+  public static int mul_1 (int[] dest, int[] x, int len, int y)
+  {
+    long yword = (long) y & 0xffffffffL;
+    long carry = 0;
+    for (int j = 0;  j < len; j++)
+      {
+        carry += ((long) x[j] & 0xffffffffL) * yword;
+        dest[j] = (int) carry;
+        carry >>>= 32;
+      }
+    return (int) carry;
+  }
+
+  /**
+   * Multiply x[0:xlen-1] and y[0:ylen-1], and
+   * write the result to dest[0:xlen+ylen-1].
+   * The destination has to have space for xlen+ylen words,
+   * even if the result might be one limb smaller.
+   * This function requires that xlen >= ylen.
+   * The destination must be distinct from either input operands.
+   * All operands are unsigned.
+   * This function is basically the same gmp's mpn_mul. */
+
+  public static void mul (int[] dest,
+			  int[] x, int xlen,
+			  int[] y, int ylen)
+  {
+    dest[xlen] = MPN.mul_1 (dest, x, xlen, y[0]);
+
+    for (int i = 1;  i < ylen; i++)
+      {
+	long yword = (long) y[i] & 0xffffffffL;
+	long carry = 0;
+	for (int j = 0;  j < xlen; j++)
+	  {
+	    carry += ((long) x[j] & 0xffffffffL) * yword
+	      + ((long) dest[i+j] & 0xffffffffL);
+	    dest[i+j] = (int) carry;
+	    carry >>>= 32;
+	  }
+	dest[i+xlen] = (int) carry;
+      }
+  }
+
+  /* Divide (unsigned long) N by (unsigned int) D.
+   * Returns (remainder << 32)+(unsigned int)(quotient).
+   * Assumes (unsigned int)(N>>32) < (unsigned int)D.
+   * Code transcribed from gmp-2.0's mpn_udiv_w_sdiv function.
+   */
+  public static long udiv_qrnnd (long N, int D)
+  {
+    long q, r;
+    long a1 = N >>> 32;
+    long a0 = N & 0xffffffffL;
+    if (D >= 0)
+      {
+	if (a1 < ((D - a1 - (a0 >>> 31)) & 0xffffffffL))
+	  {
+	    /* dividend, divisor, and quotient are nonnegative */
+	    q = N / D;
+	    r = N % D;
+	  }
+	else
+	  {
+	    /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
+	    long c = N - ((long) D << 31);
+	    /* Divide (c1*2^32 + c0) by d */
+	    q = c / D;
+	    r = c % D;
+	    /* Add 2^31 to quotient */
+	    q += 1 << 31;
+	  }
+      }
+    else
+      {
+	long b1 = D >>> 1;	/* d/2, between 2^30 and 2^31 - 1 */
+	//long c1 = (a1 >> 1); /* A/2 */
+	//int c0 = (a1 << 31) + (a0 >> 1);
+	long c = N >>> 1;
+	if (a1 < b1 || (a1 >> 1) < b1)
+	  {
+	    if (a1 < b1)
+	      {
+		q = c / b1;
+		r = c % b1;
+	      }
+	    else /* c1 < b1, so 2^31 <= (A/2)/b1 < 2^32 */
+	      {
+		c = ~(c - (b1 << 32));
+		q = c / b1;  /* (A/2) / (d/2) */
+		r = c % b1;
+		q = (~q) & 0xffffffffL;    /* (A/2)/b1 */
+		r = (b1 - 1) - r; /* r < b1 => new r >= 0 */
+	      }
+	    r = 2 * r + (a0 & 1);
+	    if ((D & 1) != 0)
+	      {
+		if (r >= q) {
+		        r = r - q;
+		} else if (q - r <= ((long) D & 0xffffffffL)) {
+                       r = r - q + D;
+        		q -= 1;
+		} else {
+                       r = r - q + D + D;
+        		q -= 2;
+		}
+	      }
+	  }
+	else				/* Implies c1 = b1 */
+	  {				/* Hence a1 = d - 1 = 2*b1 - 1 */
+	    if (a0 >= ((long)(-D) & 0xffffffffL))
+	      {
+		q = -1;
+	        r = a0 + D;
+ 	      }
+	    else
+	      {
+		q = -2;
+	        r = a0 + D + D;
+	      }
+	  }
+      }
+
+    return (r << 32) | (q & 0xFFFFFFFFl);
+  }
+
+    /** Divide divident[0:len-1] by (unsigned int)divisor.
+     * Write result into quotient[0:len-1.
+     * Return the one-word (unsigned) remainder.
+     * OK for quotient==dividend.
+     */
+
+  public static int divmod_1 (int[] quotient, int[] dividend,
+			      int len, int divisor)
+  {
+    int i = len - 1;
+    long r = dividend[i];
+    if ((r & 0xffffffffL) >= ((long)divisor & 0xffffffffL))
+      r = 0;
+    else
+      {
+	quotient[i--] = 0;
+	r <<= 32;
+      }
+
+    for (;  i >= 0;  i--)
+      {
+	int n0 = dividend[i];
+	r = (r & ~0xffffffffL) | (n0 & 0xffffffffL);
+	r = udiv_qrnnd (r, divisor);
+	quotient[i] = (int) r;
+      }
+    return (int)(r >> 32);
+  }
+
+  /* Subtract x[0:len-1]*y from dest[offset:offset+len-1].
+   * All values are treated as if unsigned.
+   * @return the most significant word of
+   * the product, minus borrow-out from the subtraction.
+   */
+  public static int submul_1 (int[] dest, int offset, int[] x, int len, int y)
+  {
+    long yl = (long) y & 0xffffffffL;
+    int carry = 0;
+    int j = 0;
+    do
+      {
+	long prod = ((long) x[j] & 0xffffffffL) * yl;
+	int prod_low = (int) prod;
+	int prod_high = (int) (prod >> 32);
+	prod_low += carry;
+	// Invert the high-order bit, because: (unsigned) X > (unsigned) Y
+	// iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
+	carry = ((prod_low ^ 0x80000000) < (carry ^ 0x80000000) ? 1 : 0)
+	  + prod_high;
+	int x_j = dest[offset+j];
+	prod_low = x_j - prod_low;
+	if ((prod_low ^ 0x80000000) > (x_j ^ 0x80000000))
+	  carry++;
+	dest[offset+j] = prod_low;
+      }
+    while (++j < len);
+    return carry;
+  }
+
+  /** Divide zds[0:nx] by y[0:ny-1].
+   * The remainder ends up in zds[0:ny-1].
+   * The quotient ends up in zds[ny:nx].
+   * Assumes:  nx>ny.
+   * (int)y[ny-1] < 0  (i.e. most significant bit set)
+   */
+
+  public static void divide (int[] zds, int nx, int[] y, int ny)
+  {
+    // This is basically Knuth's formulation of the classical algorithm,
+    // but translated from in scm_divbigbig in Jaffar's SCM implementation.
+
+    // Correspondance with Knuth's notation:
+    // Knuth's u[0:m+n] == zds[nx:0].
+    // Knuth's v[1:n] == y[ny-1:0]
+    // Knuth's n == ny.
+    // Knuth's m == nx-ny.
+    // Our nx == Knuth's m+n.
+
+    // Could be re-implemented using gmp's mpn_divrem:
+    // zds[nx] = mpn_divrem (&zds[ny], 0, zds, nx, y, ny).
+
+    int j = nx;
+    do
+      {                          // loop over digits of quotient
+	// Knuth's j == our nx-j.
+	// Knuth's u[j:j+n] == our zds[j:j-ny].
+	int qhat;  // treated as unsigned
+	if (zds[j]==y[ny-1])
+	  qhat = -1;  // 0xffffffff
+	else
+	  {
+	    long w = (((long)(zds[j])) << 32) + ((long)zds[j-1] & 0xffffffffL);
+	    qhat = (int) udiv_qrnnd (w, y[ny-1]);
+	  }
+	if (qhat != 0)
+	  {
+	    int borrow = submul_1 (zds, j - ny, y, ny, qhat);
+	    int save = zds[j];
+	    long num = ((long)save&0xffffffffL) - ((long)borrow&0xffffffffL);
+            while (num != 0)
+	      {
+		qhat--;
+		long carry = 0;
+		for (int i = 0;  i < ny; i++)
+		  {
+		    carry += ((long) zds[j-ny+i] & 0xffffffffL)
+		      + ((long) y[i] & 0xffffffffL);
+		    zds[j-ny+i] = (int) carry;
+		    carry >>>= 32;
+		  }
+		zds[j] += carry;
+		num = carry - 1;
+	      }
+	  }
+	zds[j] = qhat;
+      } while (--j >= ny);
+  }
+
+  /** Number of digits in the conversion base that always fits in a word.
+   * For example, for base 10 this is 9, since 10**9 is the
+   * largest number that fits into a words (assuming 32-bit words).
+   * This is the same as gmp's __mp_bases[radix].chars_per_limb.
+   * @param radix the base
+   * @return number of digits */
+  public static int chars_per_word (int radix)
+  {
+    if (radix < 10)
+      {
+	if (radix < 8)
+	  {
+	    if (radix <= 2)
+	      return 32;
+	    else if (radix == 3)
+	      return 20;
+	    else if (radix == 4)
+	      return 16;
+	    else
+	      return 18 - radix;
+	  }
+	else
+	  return 10;
+      }
+    else if (radix < 12)
+      return 9;
+    else if (radix <= 16)
+      return 8;
+    else if (radix <= 23)
+      return 7;
+    else if (radix <= 40)
+      return 6;
+    // The following are conservative, but we don't care.
+    else if (radix <= 256)
+      return 4;
+    else
+      return 1;
+  }
+
+  /** Count the number of leading zero bits in an int. */
+  public static int count_leading_zeros (int i)
+  {
+    if (i == 0)
+      return 32;
+    int count = 0;
+    for (int k = 16;  k > 0;  k = k >> 1) {
+      int j = i >>> k;
+      if (j == 0)
+	count += k;
+      else
+	i = j;
+    }
+    return count;
+  }
+
+  public static int set_str (int dest[], byte[] str, int str_len, int base)
+  {
+    int size = 0;
+    if ((base & (base - 1)) == 0)
+      {
+	// The base is a power of 2.  Read the input string from
+	// least to most significant character/digit.  */
+
+	int next_bitpos = 0;
+	int bits_per_indigit = 0;
+	for (int i = base; (i >>= 1) != 0; ) bits_per_indigit++;
+	int res_digit = 0;
+
+	for (int i = str_len;  --i >= 0; )
+	  {
+	    int inp_digit = str[i];
+	    res_digit |= inp_digit << next_bitpos;
+	    next_bitpos += bits_per_indigit;
+	    if (next_bitpos >= 32)
+	      {
+		dest[size++] = res_digit;
+		next_bitpos -= 32;
+		res_digit = inp_digit >> (bits_per_indigit - next_bitpos);
+	      }
+	  }
+
+	if (res_digit != 0)
+	  dest[size++] = res_digit;
+      }
+    else
+      {
+	// General case.  The base is not a power of 2.
+	int indigits_per_limb = MPN.chars_per_word (base);
+	int str_pos = 0;
+
+	while (str_pos < str_len)
+	  {
+	    int chunk = str_len - str_pos;
+	    if (chunk > indigits_per_limb)
+	      chunk = indigits_per_limb;
+	    int res_digit = str[str_pos++];
+	    int big_base = base;
+
+	    while (--chunk > 0)
+	      {
+		res_digit = res_digit * base + str[str_pos++];
+		big_base *= base;
+	      }
+
+	    int cy_limb;
+	    if (size == 0)
+	      cy_limb = res_digit;
+	    else
+	      {
+		cy_limb = MPN.mul_1 (dest, dest, size, big_base);
+		cy_limb += MPN.add_1 (dest, dest, size, res_digit);
+	      }
+	    if (cy_limb != 0)
+	      dest[size++] = cy_limb;
+	  }
+       }
+    return size;
+  }
+
+  /** Compare x[0:size-1] with y[0:size-1], treating them as unsigned integers.
+   * @result -1, 0, or 1 depending on if x<y, x==y, or x>y.
+   * This is basically the same as gmp's mpn_cmp function.
+   */
+  public static int cmp (int[] x, int[] y, int size)
+  {
+    while (--size >= 0)
+      {
+	int x_word = x[size];
+	int y_word = y[size];
+	if (x_word != y_word)
+	  {
+	    // Invert the high-order bit, because:
+	    // (unsigned) X > (unsigned) Y iff
+	    // (int) (X^0x80000000) > (int) (Y^0x80000000).
+	    return (x_word ^ 0x80000000) > (y_word ^0x80000000) ? 1 : -1;
+	  }
+      }
+    return 0;
+  }
+
+  /** Compare x[0:xlen-1] with y[0:ylen-1], treating them as unsigned integers.
+   * @result -1, 0, or 1 depending on if x<y, x==y, or x>y.
+   */
+  public static int cmp (int[] x, int xlen, int[] y, int ylen)
+  {
+    return xlen > ylen ? 1 : xlen < ylen ? -1 : cmp (x, y, xlen);
+  }
+
+  /* Shift x[x_start:x_start+len-1]count bits to the "right"
+   * (i.e. divide by 2**count).
+   * Store the len least significant words of the result at dest.
+   * The bits shifted out to the right are returned.
+   * OK if dest==x.
+   * Assumes: 0 < count < 32
+   */
+
+  public static int rshift (int[] dest, int[] x, int x_start,
+			    int len, int count)
+  {
+    int count_2 = 32 - count;
+    int low_word = x[x_start];
+    int retval = low_word << count_2;
+    int i = 1;
+    for (; i < len;  i++)
+      {
+	int high_word = x[x_start+i];
+	dest[i-1] = (low_word >>> count) | (high_word << count_2);
+	low_word = high_word;
+      }
+    dest[i-1] = low_word >>> count;
+    return retval;
+  }
+
+  /** Return the long-truncated value of right shifting.
+  * @param x a two's-complement "bignum"
+  * @param len the number of significant words in x
+  * @param count the shift count
+  * @return (long)(x[0..len-1] >> count).
+  */
+  public static long rshift_long (int[] x, int len, int count)
+  {
+    int wordno = count >> 5;
+    count &= 31;
+    int sign = x[len-1] < 0 ? -1 : 0;
+    int w0 = wordno >= len ? sign : x[wordno];
+    wordno++;
+    int w1 = wordno >= len ? sign : x[wordno];
+    if (count != 0)
+      {
+	wordno++;
+	int w2 = wordno >= len ? sign : x[wordno];
+	w0 = (w0 >>> count) | (w1 << (32-count));
+	w1 = (w1 >>> count) | (w2 << (32-count));
+      }
+    return ((long)w1 << 32) | ((long)w0 & 0xffffffffL);
+  }
+
+  /* Shift x[0:len-1]count bits to the "right" (i.e. divide by 2**count).
+   * Store the len least significant words of the result at dest.
+   * OK if dest==x.
+   * OK if count > 32 (but must be >= 0).
+   */
+  public static void rshift (int[] dest, int[] x, int len, int count)
+  {
+    int word_count = count >> 5;
+    count &= 31;
+    rshift (dest, x, word_count, len, count);
+    while (word_count < len)
+      dest[word_count++] = 0;
+  }
+
+  /* Shift x[0:len-1] left by count bits, and store the len least
+   * significant words of the result in dest[d_offset:d_offset+len-1].
+   * Return the bits shifted out from the most significant digit.
+   * Assumes 0 < count < 32.
+   * OK if dest==x.
+   */
+
+  public static int lshift (int[] dest, int d_offset,
+			    int[] x, int len, int count)
+  {
+    int count_2 = 32 - count;
+    int i = len - 1;
+    int high_word = x[i];
+    int retval = high_word >>> count_2;
+    d_offset++;
+    while (--i >= 0)
+      {
+	int low_word = x[i];
+	dest[d_offset+i] = (high_word << count) | (low_word >>> count_2);
+	high_word = low_word;
+      }
+    dest[d_offset+i] = high_word << count;
+    return retval;
+  }
+
+  /** Return least i such that word&(1<<i). Assumes word!=0. */
+
+  static int findLowestBit (int word)
+  {
+    int i = 0;
+    while ((word & 0xF) == 0)
+      {
+	word >>= 4;
+	i += 4;
+      }
+    if ((word & 3) == 0)
+      {
+	word >>= 2;
+	i += 2;
+      }
+    if ((word & 1) == 0)
+      i += 1;
+    return i;
+  }
+
+  /** Return least i such that words & (1<<i). Assumes there is such an i. */
+
+  static int findLowestBit (int[] words)
+  {
+    for (int i = 0;  ; i++)
+      {
+	if (words[i] != 0)
+	  return 32 * i + findLowestBit (words[i]);
+      }
+  }
+
+  /** Calculate Greatest Common Divisior of x[0:len-1] and y[0:len-1].
+    * Assumes both arguments are non-zero.
+    * Leaves result in x, and returns len of result.
+    * Also destroys y (actually sets it to a copy of the result). */
+
+  public static int gcd (int[] x, int[] y, int len)
+  {
+    int i, word;
+    // Find sh such that both x and y are divisible by 2**sh.
+    for (i = 0; ; i++)
+      {
+	word = x[i] | y[i];
+	if (word != 0)
+	  {
+	    // Must terminate, since x and y are non-zero.
+	    break;
+	  }
+      }
+    int initShiftWords = i;
+    int initShiftBits = findLowestBit (word);
+    // Logically: sh = initShiftWords * 32 + initShiftBits
+
+    // Temporarily devide both x and y by 2**sh.
+    len -= initShiftWords;
+    MPN.rshift (x, x, initShiftWords, len, initShiftBits);
+    MPN.rshift (y, y, initShiftWords, len, initShiftBits);
+
+    int[] odd_arg; /* One of x or y which is odd. */
+    int[] other_arg; /* The other one can be even or odd. */
+    if ((x[0] & 1) != 0)
+      {
+	odd_arg = x;
+	other_arg = y;
+      }
+    else
+      {
+	odd_arg = y;
+	other_arg = x;
+      }
+
+    for (;;)
+      {
+	// Shift other_arg until it is odd; this doesn't
+	// affect the gcd, since we divide by 2**k, which does not
+	// divide odd_arg.
+	for (i = 0; other_arg[i] == 0; ) i++;
+	if (i > 0)
+	  {
+	    int j;
+	    for (j = 0; j < len-i; j++)
+		other_arg[j] = other_arg[j+i];
+	    for ( ; j < len; j++)
+	      other_arg[j] = 0;
+	  }
+	i = findLowestBit(other_arg[0]);
+	if (i > 0)
+	  MPN.rshift (other_arg, other_arg, 0, len, i);
+
+	// Now both odd_arg and other_arg are odd.
+
+	// Subtract the smaller from the larger.
+	// This does not change the result, since gcd(a-b,b)==gcd(a,b).
+	i = MPN.cmp(odd_arg, other_arg, len);
+	if (i == 0)
+	    break;
+	if (i > 0)
+	  { // odd_arg > other_arg
+	    MPN.sub_n (odd_arg, odd_arg, other_arg, len);
+	    // Now odd_arg is even, so swap with other_arg;
+	    int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp;
+	  }
+	else
+	  { // other_arg > odd_arg
+	    MPN.sub_n (other_arg, other_arg, odd_arg, len);
+	}
+	while (odd_arg[len-1] == 0 && other_arg[len-1] == 0)
+	  len--;
+    }
+    if (initShiftWords + initShiftBits > 0)
+      {
+	if (initShiftBits > 0)
+	  {
+	    int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits);
+	    if (sh_out != 0)
+	      x[(len++)+initShiftWords] = sh_out;
+	  }
+	else
+	  {
+	    for (i = len; --i >= 0;)
+	      x[i+initShiftWords] = x[i];
+	  }
+	for (i = initShiftWords;  --i >= 0; )
+	  x[i] = 0;
+	len += initShiftWords;
+      }
+    return len;
+  }
+
+  public static int intLength (int i)
+  {
+    return 32 - count_leading_zeros (i < 0 ? ~i : i);
+  }
+
+  /** Calcaulte the Common Lisp "integer-length" function.
+   * Assumes input is canonicalized:  len==IntNum.wordsNeeded(words,len) */
+  public static int intLength (int[] words, int len)
+  {
+    len--;
+    return intLength (words[len]) + 32 * len;
+  }
+
+  /* DEBUGGING:
+  public static void dprint (IntNum x)
+  {
+    if (x.words == null)
+      System.err.print(Long.toString((long) x.ival & 0xffffffffL, 16));
+    else
+      dprint (System.err, x.words, x.ival);
+  }
+  public static void dprint (int[] x) { dprint (System.err, x, x.length); }
+  public static void dprint (int[] x, int len) { dprint (System.err, x, len); }
+  public static void dprint (java.io.PrintStream ps, int[] x, int len)
+  {
+    ps.print('(');
+    for (int i = 0;  i < len; i++)
+      {
+	if (i > 0)
+	  ps.print (' ');
+	ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16));
+      }
+    ps.print(')');
+  }
+  */
+}

libjava/java/math/BigInteger.java

+// BigInteger.java -- an arbitrary-precision integer
+
+/* Copyright (C) 1999, 2000  Red Hat, Inc.
+
+   This file is part of libgcj.
+
+This software is copyrighted work licensed under the terms of the
+Libgcj License.  Please consult the file "LIBGCJ_LICENSE" for
+details.  */
+
+package java.math;
+import gnu.gcj.math.*;
+import java.util.Random;
+
+/**
+ * @author Warren Levy <warrenl@cygnus.com>
+ * @date December 20, 1999.
+ */
+
+/**
+ * Written using on-line Java Platform 1.2 API Specification, as well
+ * as "The Java Class Libraries", 2nd edition (Addison-Wesley, 1998).
+ * 
+ * Based primarily on IntNum.java by Per Bothner <per@bothner.com>
+ * (found in Kawa 1.6.62).
+ *
+ * Status:  Believed complete and correct.
+ */
+
+public class BigInteger extends Number implements Comparable
+{
+  /** All integers are stored in 2's-complement form.
+   * If words == null, the ival is the value of this BigInteger.
+   * Otherwise, the first ival elements of words make the value
+   * of this BigInteger, stored in little-endian order, 2's-complement form. */
+  public int ival;
+  public int[] words;
+
+
+  /** We pre-allocate integers in the range minFixNum..maxFixNum. */
+  private static final int minFixNum = -100;
+  private static final int maxFixNum = 1024;
+  private static final int numFixNum = maxFixNum-minFixNum+1;
+  private static final BigInteger[] smallFixNums = new BigInteger[numFixNum];
+
+  static {
+    for (int i = numFixNum;  --i >= 0; )
+      smallFixNums[i] = new BigInteger(i + minFixNum);
+  }
+
+  // JDK1.2
+  public static final BigInteger ZERO = smallFixNums[-minFixNum];
+
+  // JDK1.2
+  public static final BigInteger ONE = smallFixNums[1 - minFixNum];
+
+  /* Rounding modes: */
+  private static final int FLOOR = 1;
+  private static final int CEILING = 2;
+  private static final int TRUNCATE = 3;
+  private static final int ROUND = 4;
+
+  private BigInteger()
+  {
+  }
+
+  /* Create a new (non-shared) BigInteger, and initialize to an int. */
+  private BigInteger(int value)
+  {
+    ival = value;
+  }
+
+  /* Create a new (non-shared) BigInteger, and initialize from a byte array. */
+  public BigInteger(byte[] val)
+  {
+    if (val == null || val.length < 1)
+      throw new NumberFormatException();
+
+    words = byteArrayToIntArray(val, val[0] < 0 ? -1 : 0);
+    BigInteger result = make(words, words.length);
+    this.ival = result.ival;
+    this.words = result.words;
+  }
+
+  public BigInteger(int signum, byte[] magnitude)
+  {
+    if (magnitude == null || signum > 1 || signum < -1)
+      throw new NumberFormatException();
+
+    if (signum == 0)
+      {
+	int i;
+	for (i = magnitude.length - 1; i >= 0 && magnitude[i] == 0; --i)
+	  ;
+	if (i >= 0)
+	  throw new NumberFormatException();
+        return;
+      }
+
+    // Magnitude is always positive, so don't ever pass a sign of -1.
+    words = byteArrayToIntArray(magnitude, 0);
+    BigInteger result = make(words, words.length);
+    this.ival = result.ival;
+    this.words = result.words;
+
+    if (signum < 0)
+      setNegative();
+  }
+
+  public BigInteger(int numBits, Random rnd)
+  {
+    if (numBits < 0)
+      throw new IllegalArgumentException();
+
+    // Result is always positive so tack on an extra zero word, it will be
+    // canonicalized out later if necessary.
+    int nwords = numBits / 32 + 2;
+    words = new int[nwords];
+    words[--nwords] = 0;
+    words[--nwords] = rnd.nextInt() >>> (numBits % 32);
+    while (--nwords >= 0)
+      words[nwords] = rnd.nextInt();
+
+    BigInteger result = make(words, words.length);
+    this.ival = result.ival;
+    this.words = result.words;
+  }
+
+
+  /** Return a (possibly-shared) BigInteger with a given long value. */
+  private static BigInteger make(long value)
+  {
+    if (value >= minFixNum && value <= maxFixNum)
+      return smallFixNums[(int)value - minFixNum];
+    int i = (int) value;
+    if ((long)i == value)
+      return new BigInteger(i);
+    BigInteger result = alloc(2);
+    result.ival = 2;
+    result.words[0] = i;
+    result.words[1] = (int) (value >> 32);
+    return result;
+  }
+
+  // FIXME: Could simply rename 'make' method above as valueOf while
+  // changing all instances of 'make'.  Don't do this until this class
+  // is done as the Kawa class this is based on has 'make' methods
+  // with other parameters; wait to see if they are used in BigInteger.
+  public static BigInteger valueOf(long val)
+  {
+    return make(val);
+  }
+
+  /** Make a canonicalized BigInteger from an array of words.
+   * The array may be reused (without copying). */
+  private static BigInteger make(int[] words, int len)
+  {
+    if (words == null)
+      return make(len);
+    len = BigInteger.wordsNeeded(words, len);
+    if (len <= 1)
+      return len == 0 ? ZERO : make(words[0]);
+    BigInteger num = new BigInteger();
+    num.words = words;
+    num.ival = len;
+    return num;
+  }
+
+  /** Convert a big-endian byte array to a little-endian array of words. */
+  private static int[] byteArrayToIntArray(byte[] bytes, int sign)
+  {
+    // Determine number of words needed.
+    int[] words = new int[(bytes.length + 3) / 4 + 1];
+    int nwords = words.length;
+
+    // For simplicity, tack on an extra word of sign at the front,
+    // it will be canonicalized out later. */
+    words[--nwords] = sign;
+
+    // Create a int out of modulo 4 high order bytes.
+    int bptr = 0;
+    int word = sign;
+    for (int i = bytes.length % 4; i > 0; --i, bptr++)
+      word = (word << 8) | (((int) bytes[bptr]) & 0xff);
+    words[--nwords] = word;
+
+    // Elements remaining in byte[] is a  multiple of 4.
+    while (nwords > 0)
+      words[--nwords] = bytes[bptr++] << 24 |
+			(((int) bytes[bptr++]) & 0xff) << 16 |
+			(((int) bytes[bptr++]) & 0xff) << 8 |
+			(((int) bytes[bptr++]) & 0xff);
+    return words;
+  }
+
+  /** Allocate a new non-shared BigInteger.
+   * @param nwords number of words to allocate
+   */
+  private static BigInteger alloc(int nwords)
+  {
+    if (nwords <= 1)
+      return new BigInteger();
+    BigInteger result = new BigInteger();
+    result.words = new int[nwords];
+    return result;
+  }
+
+  /** Change words.length to nwords.
+   * We allow words.length to be upto nwords+2 without reallocating.
+   */
+  private void realloc(int nwords)
+  {
+    if (nwords == 0)
+      {
+	if (words != null)
+	  {
+	    if (ival > 0)
+	      ival = words[0];
+	    words = null;
+	  }
+      }
+    else if (words == null
+	     || words.length < nwords
+	     || words.length > nwords + 2)
+      {
+	int[] new_words = new int [nwords];
+	if (words == null)
+	  {
+	    new_words[0] = ival;
+	    ival = 1;
+	  }
+	else
+	  {
+	    if (nwords < ival)
+	      ival = nwords;
+	    System.arraycopy(words, 0, new_words, 0, ival);
+	  }
+	words = new_words;
+      }
+  }
+
+  private final boolean isNegative()
+  {
+    return (words == null ? ival : words[ival - 1]) < 0;
+  }
+
+  public int signum()
+  {
+    int top = words == null ? ival : words[ival-1];
+    return top > 0 ? 1 : top < 0 ? -1 : 0;
+  }
+
+  private static int compareTo(BigInteger x, BigInteger y)
+  {
+    if (x.words == null && y.words == null)
+      return x.ival < y.ival ? -1 : x.ival > y.ival ? 1 : 0;
+    boolean x_negative = x.isNegative();
+    boolean y_negative = y.isNegative();
+    if (x_negative != y_negative)
+      return x_negative ? -1 : 1;
+    int x_len = x.words == null ? 1 : x.ival;
+    int y_len = y.words == null ? 1 : y.ival;
+    if (x_len != y_len)
+      return (x_len > y_len) != x_negative ? 1 : -1;
+    return MPN.cmp(x.words, y.words, x_len);
+  }
+
+  // JDK1.2
+  public int compareTo(Object obj)
+  {
+    if (obj instanceof BigInteger)
+      return compareTo(this, (BigInteger) obj);
+    throw new ClassCastException();
+  }
+
+  public int compareTo(BigInteger val)
+  {
+    return compareTo(this, val);
+  }
+
+  private final boolean isOdd()
+  {
+    int low = words == null ? ival : words[0];
+    return (low & 1) != 0;
+  }
+
+  private final boolean isZero()
+  {
+    return words == null && ival == 0;
+  }
+
+  private final boolean isOne()
+  {
+    return words == null && ival == 1;
+  }
+
+  private final boolean isMinusOne()
+  {
+    return words == null && ival == -1;
+  }
+
+  /** Calculate how many words are significant in words[0:len-1].
+   * Returns the least value x such that x>0 && words[0:x-1]==words[0:len-1],
+   * when words is viewed as a 2's complement integer.
+   */
+  private static int wordsNeeded(int[] words, int len)
+  {
+    int i = len;
+    if (i > 0)
+      {
+	int word = words[--i];
+	if (word == -1)
+	  {
+	    while (i > 0 && (word = words[i - 1]) < 0)
+	      {
+		i--;
+		if (word != -1) break;
+	      }
+	  }
+	else
+	  {
+	    while (word == 0 && i > 0 && (word = words[i - 1]) >= 0)  i--;
+	  }
+      }
+    return i + 1;
+  }
+
+  private BigInteger canonicalize()
+  {
+    if (words != null
+	&& (ival = BigInteger.wordsNeeded(words, ival)) <= 1)
+      {
+	if (ival == 1)
+	  ival = words[0];
+	words = null;
+      }
+    if (words == null && ival >= minFixNum && ival <= maxFixNum)
+      return smallFixNums[(int) ival - minFixNum];
+    return this;
+  }
+
+  /** Add two ints, yielding an BigInteger. */
+  private static final BigInteger add(int x, int y)
+  {
+    return BigInteger.make((long) x + (long) y);
+  }
+
+  /** Add an BigInteger and an int, yielding a new BigInteger. */
+  private static BigInteger add(BigInteger x, int y)
+  {
+    if (x.words == null)
+      return BigInteger.add(x.ival, y);
+    BigInteger result = new BigInteger(0);
+    result.setAdd(x, y);
+    return result.canonicalize();
+  }
+
+  /** Set this to the sum of x and y.
+   * OK if x==this. */
+  private void setAdd(BigInteger x, int y)
+  {
+    if (x.words == null)
+      {
+	set((long) x.ival + (long) y);
+	return;
+      }
+    int len = x.ival;
+    realloc(len + 1);
+    long carry = y;
+    for (int i = 0;  i < len;  i++)
+      {
+	carry += ((long) x.words[i] & 0xffffffffL);
+	words[i] = (int) carry;
+	carry >>= 32;
+      }
+    if (x.words[len - 1] < 0)
+      carry--;
+    words[len] = (int) carry;
+    ival = wordsNeeded(words, len + 1);
+  }
+
+  /** Destructively add an int to this. */
+  private final void setAdd(int y)
+  {
+    setAdd(this, y);
+  }
+
+  /** Destructively set the value of this to a long. */
+  private final void set(long y)
+  {
+    int i = (int) y;
+    if ((long) i == y)
+      {
+	ival = i;
+	words = null;
+      }
+    else
+      {
+	realloc(2);
+	words[0] = i;
+	words[1] = (int) (y >> 32);
+	ival = 2;
+      }
+  }
+
+  /** Destructively set the value of this to the given words.
+  * The words array is reused, not copied. */
+  private final void set(int[] words, int length)
+  {
+    this.ival = length;
+    this.words = words;
+  }
+
+  /** Destructively set the value of this to that of y. */
+  private final void set(BigInteger y)
+  {
+    if (y.words == null)
+      set(y.ival);
+    else if (this != y)
+      {
+	realloc(y.ival);
+	System.arraycopy(y.words, 0, words, 0, y.ival);
+	ival = y.ival;
+      }
+  }
+
+  /** Add two BigIntegers, yielding their sum as another BigInteger. */
+  private static BigInteger add(BigInteger x, BigInteger y, int k)
+  {
+    if (x.words == null && y.words == null)
+      return BigInteger.make((long) k * (long) y.ival + (long) x.ival);
+    if (k != 1)
+      {
+	if (k == -1)
+	  y = BigInteger.neg(y);
+	else
+	  y = BigInteger.times(y, BigInteger.make(k));
+      }
+    if (x.words == null)
+      return BigInteger.add(y, x.ival);
+    if (y.words == null)
+      return BigInteger.add(x, y.ival);
+    // Both are big
+    int len;
+    if (y.ival > x.ival)
+      { // Swap so x is longer then y.
+	BigInteger tmp = x;  x = y;  y = tmp;
+      }
+    BigInteger result = alloc(x.ival + 1);
+    int i = y.ival;
+    long carry = MPN.add_n(result.words, x.words, y.words, i);
+    long y_ext = y.words[i - 1] < 0 ? 0xffffffffL : 0;
+    for (; i < x.ival;  i++)
+      {
+	carry += ((long) x.words[i] & 0xffffffffL) + y_ext;;
+	result.words[i] = (int) carry;
+	carry >>>= 32;
+      }
+    if (x.words[i - 1] < 0)
+      y_ext--;
+    result.words[i] = (int) (carry + y_ext);
+    result.ival = i+1;
+    return result.canonicalize();
+  }
+
+  public BigInteger add(BigInteger val)
+  {
+    return add(this, val, 1);
+  }
+
+  public BigInteger subtract(BigInteger val)
+  {
+    return add(this, val, -1);
+  }
+
+  private static final BigInteger times(BigInteger x, int y)
+  {
+    if (y == 0)
+      return ZERO;
+    if (y == 1)
+      return x;
+    int[] xwords = x.words;
+    int xlen = x.ival;
+    if (xwords == null)
+      return BigInteger.make((long) xlen * (long) y);
+    boolean negative;
+    BigInteger result = BigInteger.alloc(xlen + 1);
+    if (xwords[xlen - 1] < 0)
+      {
+	negative = true;
+	negate(result.words, xwords, xlen);
+	xwords = result.words;
+      }
+    else
+      negative = false;
+    if (y < 0)
+      {
+	negative = !negative;
+	y = -y;
+      }
+    result.words[xlen] = MPN.mul_1(result.words, xwords, xlen, y);
+    result.ival = xlen + 1;
+    if (negative)
+      result.setNegative();
+    return result.canonicalize();
+  }
+
+  private static final BigInteger times(BigInteger x, BigInteger y)
+  {
+    if (y.words == null)
+      return times(x, y.ival);
+    if (x.words == null)
+      return times(y, x.ival);
+    boolean negative = false;
+    int[] xwords;
+    int[] ywords;
+    int xlen = x.ival;
+    int ylen = y.ival;
+    if (x.isNegative())
+      {
+	negative = true;
+	xwords = new int[xlen];
+	negate(xwords, x.words, xlen);
+      }
+    else
+      {
+	negative = false;
+	xwords = x.words;
+      }
+    if (y.isNegative())
+      {
+	negative = !negative;
+	ywords = new int[ylen];
+	negate(ywords, y.words, ylen);
+      }
+    else
+      ywords = y.words;
+    // Swap if x is shorter then y.
+    if (xlen < ylen)
+      {
+	int[] twords = xwords;  xwords = ywords;  ywords = twords;
+	int tlen = xlen;  xlen = ylen;  ylen = tlen;
+      }
+    BigInteger result = BigInteger.alloc(xlen+ylen);
+    MPN.mul(result.words, xwords, xlen, ywords, ylen);
+    result.ival = xlen+ylen;
+    if (negative)
+      result.setNegative();
+    return result.canonicalize();
+  }
+
+  public BigInteger multiply(BigInteger y)
+  {
+    return times(this, y);
+  }
+
+  private static void divide(long x, long y,
+			     BigInteger quotient, BigInteger remainder,
+			     int rounding_mode)
+  {
+    boolean xNegative, yNegative;
+    if (x < 0)
+      {
+	xNegative = true;
+	if (x == Long.MIN_VALUE)
+	  {
+	    divide(BigInteger.make(x), BigInteger.make(y),
+		   quotient, remainder, rounding_mode);
+	    return;
+	  }
+	x = -x;
+      }
+    else
+      xNegative = false;
+
+    if (y < 0)
+      {
+	yNegative = true;
+	if (y == Long.MIN_VALUE)
+	  {
+	    if (rounding_mode == TRUNCATE)
+	      { // x != Long.Min_VALUE implies abs(x) < abs(y)
+		if (quotient != null)
+		  quotient.set(0);
+		if (remainder != null)
+		  remainder.set(x);
+	      }
+	    else
+	      divide(BigInteger.make(x), BigInteger.make(y),
+		      quotient, remainder, rounding_mode);
+	    return;
+	  }
+	y = -y;
+      }
+    else
+      yNegative = false;
+
+    long q = x / y;
+    long r = x % y;
+    boolean qNegative = xNegative ^ yNegative;
+
+    boolean add_one = false;
+    if (r != 0)
+      {
+	switch (rounding_mode)
+	  {
+	  case TRUNCATE:
+	    break;
+	  case CEILING:
+	  case FLOOR:
+	    if (qNegative == (rounding_mode == FLOOR))
+	      add_one = true;
+	    break;
+	  case ROUND:
+	    add_one = r > ((y - (q & 1)) >> 1);
+	    break;
+	  }
+      }
+    if (quotient != null)
+      {
+	if (add_one)
+	  q++;
+	if (qNegative)
+	  q = -q;
+	quotient.set(q);
+      }
+    if (remainder != null)
+      {
+	// The remainder is by definition: X-Q*Y
+	if (add_one)
+	  {
+	    // Subtract the remainder from Y.
+	    r = y - r;
+	    // In this case, abs(Q*Y) > abs(X).
+	    // So sign(remainder) = -sign(X).
+	    xNegative = ! xNegative;
+	  }
+	else
+	  {
+	    // If !add_one, then: abs(Q*Y) <= abs(X).
+	    // So sign(remainder) = sign(X).
+	  }
+	if (xNegative)
+	  r = -r;
+	remainder.set(r);
+      }
+  }
+
+  /** Divide two integers, yielding quotient and remainder.
+   * @param x the numerator in the division
+   * @param y the denominator in the division
+   * @param quotient is set to the quotient of the result (iff quotient!=null)
+   * @param remainder is set to the remainder of the result
+   *  (iff remainder!=null)
+   * @param rounding_mode one of FLOOR, CEILING, TRUNCATE, or ROUND.
+   */
+  private static void divide(BigInteger x, BigInteger y,
+			     BigInteger quotient, BigInteger remainder,
+			     int rounding_mode)
+  {
+    if ((x.words == null || x.ival <= 2)
+	&& (y.words == null || y.ival <= 2))
+      {
+	long x_l = x.longValue();
+	long y_l = y.longValue();
+	if (x_l != Long.MIN_VALUE && y_l != Long.MIN_VALUE)
+	  {
+	    divide(x_l, y_l, quotient, remainder, rounding_mode);
+	    return;
+	  }
+      }
+
+    boolean xNegative = x.isNegative();
+    boolean yNegative = y.isNegative();
+    boolean qNegative = xNegative ^ yNegative;
+
+    int ylen = y.words == null ? 1 : y.ival;
+    int[] ywords = new int[ylen];
+    y.getAbsolute(ywords);
+    while (ylen > 1 && ywords[ylen - 1] == 0)  ylen--;
+
+    int xlen = x.words == null ? 1 : x.ival;
+    int[] xwords = new int[xlen+2];
+    x.getAbsolute(xwords);
+    while (xlen > 1 && xwords[xlen-1] == 0)  xlen--;
+
+    int qlen, rlen;
+
+    int cmpval = MPN.cmp(xwords, xlen, ywords, ylen);
+    if (cmpval < 0)  // abs(x) < abs(y)
+      { // quotient = 0;  remainder = num.
+	int[] rwords = xwords;  xwords = ywords;  ywords = rwords;
+	rlen = xlen;  qlen = 1;  xwords[0] = 0;
+      }
+    else if (cmpval == 0)  // abs(x) == abs(y)
+      {
+	xwords[0] = 1;  qlen = 1;  // quotient = 1
+	ywords[0] = 0;  rlen = 1;  // remainder = 0;
+      }
+    else if (ylen == 1)
+      {
+	qlen = xlen;
+	rlen = 1;
+	ywords[0] = MPN.divmod_1(xwords, xwords, xlen, ywords[0]);
+      }
+    else  // abs(x) > abs(y)
+      {
+	// Normalize the denominator, i.e. make its most significant bit set by
+	// shifting it normalization_steps bits to the left.  Also shift the
+	// numerator the same number of steps (to keep the quotient the same!).
+
+	int nshift = MPN.count_leading_zeros(ywords[ylen - 1]);
+	if (nshift != 0)
+	  {
+	    // Shift up the denominator setting the most significant bit of
+	    // the most significant word.
+	    MPN.lshift(ywords, 0, ywords, ylen, nshift);
+
+	    // Shift up the numerator, possibly introducing a new most
+	    // significant word.
+	    int x_high = MPN.lshift(xwords, 0, xwords, xlen, nshift);
+	    xwords[xlen++] = x_high;
+	}
+
+	if (xlen == ylen)
+	  xwords[xlen++] = 0;
+	MPN.divide(xwords, xlen, ywords, ylen);
+	rlen = ylen;
+	if (remainder != null || rounding_mode != TRUNCATE)
+	  {
+	    if (nshift == 0)
+	      System.arraycopy(xwords, 0, ywords, 0, rlen);
+	    else
+	      MPN.rshift(ywords, xwords, 0, rlen, nshift);
+	  }
+
+	qlen = xlen + 1 - ylen;
+	if (quotient != null)
+	  {
+	    for (int i = 0;  i < qlen;  i++)
+	      xwords[i] = xwords[i+ylen];
+	  }
+      }
+
+    // Now the quotient is in xwords, and the remainder is in ywords.
+
+    boolean add_one = false;
+    if (rlen > 1 || ywords[0] != 0)
+      { // Non-zero remainder i.e. in-exact quotient.
+	switch (rounding_mode)
+	  {
+	  case TRUNCATE:
+	    break;
+	  case CEILING:
+	  case FLOOR:
+	    if (qNegative == (rounding_mode == FLOOR))
+	      add_one = true;
+	    break;
+	  case ROUND:
+	    // int cmp = compareTo(remainder<<1, abs(y));
+	    BigInteger tmp = remainder == null ? new BigInteger() : remainder;
+	    tmp.set(ywords, rlen);
+	    tmp = shift(tmp, 1);
+	    if (yNegative)
+	      tmp.setNegative();
+	    int cmp = compareTo(tmp, y);
+	    // Now cmp == compareTo(sign(y)*(remainder<<1), y)
+	    if (yNegative)
+	      cmp = -cmp;
+	    add_one = (cmp == 1) || (cmp == 0 && (xwords[0]&1) != 0);
+	  }
+      }
+    if (quotient != null)
+      {
+	quotient.set(xwords, qlen);
+	if (qNegative)
+	  {
+	    if (add_one)  // -(quotient + 1) == ~(quotient)
+	      quotient.setInvert();
+	    else
+	      quotient.setNegative();
+	  }
+	else if (add_one)
+	  quotient.setAdd(1);
+      }
+    if (remainder != null)
+      {
+	// The remainder is by definition: X-Q*Y
+	remainder.set(ywords, rlen);
+	if (add_one)
+	  {
+	    // Subtract the remainder from Y:
+	    // abs(R) = abs(Y) - abs(orig_rem) = -(abs(orig_rem) - abs(Y)).
+	    BigInteger tmp;
+	    if (y.words == null)
+	      {
+		tmp = remainder;
+		tmp.set(yNegative ? ywords[0] + y.ival : ywords[0] - y.ival);
+	      }
+	    else
+	      tmp = BigInteger.add(remainder, y, yNegative ? 1 : -1);
+	    // Now tmp <= 0.
+	    // In this case, abs(Q) = 1 + floor(abs(X)/abs(Y)).
+	    // Hence, abs(Q*Y) > abs(X).
+	    // So sign(remainder) = -sign(X).
+	    if (xNegative)
+	      remainder.setNegative(tmp);
+	    else
+	      remainder.set(tmp);
+	  }
+	else
+	  {
+	    // If !add_one, then: abs(Q*Y) <= abs(X).
+	    // So sign(remainder) = sign(X).
+	    if (xNegative)
+	      remainder.setNegative();
+	  }
+      }
+  }
+
+  public BigInteger divide(BigInteger val)
+  {
+    if (val.isZero())
+      throw new ArithmeticException("divisor is zero");
+
+    BigInteger quot = new BigInteger();
+    divide(this, val, quot, null, TRUNCATE);
+    return quot.canonicalize();
+  }
+
+  public BigInteger remainder(BigInteger val)
+  {
+    if (val.isZero())
+      throw new ArithmeticException("divisor is zero");
+
+    BigInteger rem = new BigInteger();
+    divide(this, val, null, rem, TRUNCATE);
+    return rem.canonicalize();
+  }
+
+  public BigInteger[] divideAndRemainder(BigInteger val)
+  {
+    if (val.isZero())
+      throw new ArithmeticException("divisor is zero");
+
+    BigInteger[] result = new BigInteger[2];
+    result[0] = new BigInteger();
+    result[1] = new BigInteger();
+    divide(this, val, result[0], result[1], TRUNCATE);
+    result[0].canonicalize();
+    result[1].canonicalize();
+    return result;
+  }
+
+  public BigInteger mod(BigInteger m)
+  {
+    if (m.isNegative() || m.isZero())
+      throw new ArithmeticException("non-positive modulus");
+
+    BigInteger rem = new BigInteger();
+    divide(this, m, null, rem, FLOOR);
+    return rem.canonicalize();
+  }
+
+  /** Calculate power for BigInteger exponents.
+   * @param y exponent assumed to be non-negative. */
+  private BigInteger pow(BigInteger y)
+  {
+    if (isOne())
+      return this;
+    if (isMinusOne())
+      return y.isOdd () ? this : ONE;
+    if (y.words == null && y.ival >= 0)
+      return pow(y.ival);
+
+    // Assume exponent is non-negative.
+    if (isZero())
+      return this;
+
+    // Implemented by repeated squaring and multiplication.
+    BigInteger pow2 = this;
+    BigInteger r = null;
+    for (;;)  // for (i = 0;  ; i++)
+      {
+        // pow2 == x**(2**i)
+        // prod = x**(sum(j=0..i-1, (y>>j)&1))
+        if (y.isOdd())
+          r = r == null ? pow2 : times(r, pow2);  // r *= pow2
+        y = BigInteger.shift(y, -1);
+        if (y.isZero())
+          break;
+        // pow2 *= pow2;
+        pow2 = times(pow2, pow2);
+      }
+    return r == null ? ONE : r;
+  }
+
+  /** Calculate the integral power of a BigInteger.
+   * @param exponent the exponent (must be non-negative)
+   */
+  public BigInteger pow(int exponent)
+  {
+    if (exponent <= 0)
+      {
+	if (exponent == 0)
+	  return ONE;
+	else
+	  throw new ArithmeticException("negative exponent");
+      }
+    if (isZero())
+      return this;
+    int plen = words == null ? 1 : ival;  // Length of pow2.
+    int blen = ((bitLength() * exponent) >> 5) + 2 * plen;
+    boolean negative = isNegative() && (exponent & 1) != 0;
+    int[] pow2 = new int [blen];
+    int[] rwords = new int [blen];
+    int[] work = new int [blen];
+    getAbsolute(pow2);	// pow2 = abs(this);
+    int rlen = 1;
+    rwords[0] = 1; // rwords = 1;
+    for (;;)  // for (i = 0;  ; i++)
+      {
+	// pow2 == this**(2**i)
+	// prod = this**(sum(j=0..i-1, (exponent>>j)&1))
+	if ((exponent & 1) != 0)
+	  { // r *= pow2
+	    MPN.mul(work, pow2, plen, rwords, rlen);
+	    int[] temp = work;  work = rwords;  rwords = temp;
+	    rlen += plen;
+	    while (rwords[rlen - 1] == 0)  rlen--;
+	  }
+	exponent >>= 1;
+	if (exponent == 0)
+	  break;
+	// pow2 *= pow2;
+	MPN.mul(work, pow2, plen, pow2, plen);
+	int[] temp = work;  work = pow2;  pow2 = temp;  // swap to avoid a copy
+	plen *= 2;
+	while (pow2[plen - 1] == 0)  plen--;
+      }
+    if (rwords[rlen - 1] < 0)
+      rlen++;
+    if (negative)
+      negate(rwords, rwords, rlen);
+    return BigInteger.make(rwords, rlen);
+  }
+
+  private static final int[] euclidInv(int a, int b, int prevDiv)
+  {
+    // Storage for return values, plus one slot for a temp int (see below).
+    int[] xy;
+
+    if (b == 0)
+      throw new ArithmeticException("not invertible");
+    else if (b == 1)
+      {
+	// Success:  values are indeed invertible!
+	// Bottom of the recursion reached; start unwinding.
+        xy = new int[3];
+	xy[0] = -prevDiv;
+	xy[1] = 1;
+	return xy;
+      }
+
+    xy = euclidInv(b, a % b, a / b);	// Recursion happens here.
+
+    // xy[2] is just temp storage for intermediate results in the following
+    // calculation.  This saves us a bit of space over having an int
+    // allocated at every level of this recursive method.
+    xy[2] = xy[0];
+    xy[0] = xy[2] * -prevDiv + xy[1];
+    xy[1] = xy[2];
+    return xy;
+  }
+
+  private static final BigInteger[]
+    euclidInv(BigInteger a, BigInteger b, BigInteger prevDiv)
+  {
+    // FIXME: This method could be more efficient memory-wise and should be
+    // modified as such since it is recursive.
+
+    // Storage for return values, plus one slot for a temp int (see below).
+    BigInteger[] xy;
+
+    if (b.isZero())
+      throw new ArithmeticException("not invertible");
+    else if (b.isOne())
+      {
+	// Success:  values are indeed invertible!
+	// Bottom of the recursion reached; start unwinding.
+        xy = new BigInteger[3];
+	xy[0] = neg(prevDiv);
+	xy[1] = ONE;
+	return xy;
+      }
+
+    // Recursion happens in the following conditional!
+
+    // If a just contains an int, then use integer math for the rest.
+    if (a.words == null)
+      {
+        int[] xyInt = euclidInv(b.ival, a.ival % b.ival, a.ival / b.ival);
+        xy = new BigInteger[3];
+	xy[0] = new BigInteger(xyInt[0]);
+	xy[1] = new BigInteger(xyInt[1]);
+      }
+    else
+      {
+	BigInteger rem = new BigInteger();
+	BigInteger quot = new BigInteger();
+	divide(a, b, quot, rem, FLOOR);
+        xy = euclidInv(b, rem, quot);
+      }
+
+    // xy[2] is just temp storage for intermediate results in the following
+    // calculation.  This saves us a bit of space over having a BigInteger
+    // allocated at every level of this recursive method.
+    xy[2] = xy[0];
+    xy[0] = add(xy[1], times(xy[2], prevDiv), -1);
+    xy[1] = xy[2];
+    return xy;
+  }
+
+  public BigInteger modInverse(BigInteger y)
+  {
+    if (y.isNegative() || y.isZero())
+      throw new ArithmeticException("non-positive modulo");
+
+    // Degenerate cases.
+    if (y.isOne())
+      return ZERO;
+    else if (isOne())
+      return ONE;
+
+    // Use Euclid's algorithm as in gcd() but do this recursively
+    // rather than in a loop so we can use the intermediate results as we
+    // unwind from the recursion.
+    // Used http://www.math.nmsu.edu/~crypto/EuclideanAlgo.html as reference.
+    BigInteger result = new BigInteger();
+    int xval = ival;
+    int yval = y.ival;
+    boolean swapped = false;
+
+    if (y.words == null)
+      {
+	// The result is guaranteed to be less than the modulus, y (which is
+	// an int), so simplify this by working with the int result of this
+	// modulo y.  Also, if this is negative, make it positive via modulo
+	// math.  Note that BigInteger.mod() must be used even if this is
+	// already an int as the % operator would provide a negative result if
+	// this is negative, BigInteger.mod() never returns negative values.
+	if (words != null || isNegative())
+	  xval = mod(y).ival;
+
+	// Swap values so x > y.
+	if (yval > xval)
+	  {
+	    int tmp = xval; xval = yval; yval = tmp;
+	    swapped = true;
+	  }
+	// Normally, the result is in the 2nd element of the array, but
+	// if originally x < y, then x and y were swapped and the result
+	// is in the 1st element of the array.
+	result.ival =
+	  euclidInv(yval, xval % yval, xval / yval)[swapped ? 0 : 1];
+
+	// Result can't be negative, so make it positive by adding the
+	// original modulus, y.ival (not the possibly "swapped" yval).
+	if (result.ival < 0)
+	  result.ival += y.ival;
+      }
+    else
+      {
+	BigInteger x = this;
+
+	// As above, force this to be a positive value via modulo math.
+	if (isNegative())
+	  x = mod(y);
+
+	// Swap values so x > y.
+	if (x.compareTo(y) < 0)
+	  {
+	    BigInteger tmp = x; x = y; y = tmp;
+	    swapped = true;
+	  }
+	// As above (for ints), result will be in the 2nd element unless
+	// the original x and y were swapped.
+	BigInteger rem = new BigInteger();
+	BigInteger quot = new BigInteger();
+	divide(x, y, quot, rem, FLOOR);
+	result = euclidInv(y, rem, quot)[swapped ? 0 : 1];
+
+	// Result can't be negative, so make it positive by adding the
+	// original modulus, y (which is now x if they were swapped).
+	if (result.isNegative())
+	  result = add(result, swapped ? x : y, 1);
+      }
+    
+    return result;
+  }
+
+  public BigInteger modPow(BigInteger exponent, BigInteger m)
+  {
+    if (m.isNegative() || m.isZero())
+      throw new ArithmeticException("non-positive modulo");
+
+    if (exponent.isNegative())
+      return modInverse(m);
+    if (exponent.isOne())
+      return mod(m);
+
+    return pow(exponent).mod(m);
+  }
+
+  /** Calculate Greatest Common Divisor for non-negative ints. */
+  private static final int gcd(int a, int b)
+  {
+    // Euclid's algorithm, copied from libg++.
+    if (b > a)
+      {
+	int tmp = a; a = b; b = tmp;
+      }
+    for(;;)
+      {
+	if (b == 0)
+	  return a;
+	else if (b == 1)
+	  return b;
+	else
+	  {
+	    int tmp = b;
+	    b = a % b;
+	    a = tmp;
+	  }
+      }
+  }
+
+  public BigInteger gcd(BigInteger y)
+  {
+    int xval = ival;
+    int yval = y.ival;
+    if (words == null)
+      {
+	if (xval == 0)
+	  return BigInteger.abs(y);
+	if (y.words == null
+	    && xval != Integer.MIN_VALUE && yval != Integer.MIN_VALUE)
+	  {
+	    if (xval < 0)
+	      xval = -xval;
+	    if (yval < 0)
+	      yval = -yval;
+	    return BigInteger.make(BigInteger.gcd(xval, yval));
+	  }
+	xval = 1;
+      }
+    if (y.words == null)
+      {
+	if (yval == 0)
+	  return BigInteger.abs(this);
+	yval = 1;
+      }
+    int len = (xval > yval ? xval : yval) + 1;
+    int[] xwords = new int[len];
+    int[] ywords = new int[len];
+    getAbsolute(xwords);
+    y.getAbsolute(ywords);
+    len = MPN.gcd(xwords, ywords, len);
+    BigInteger result = new BigInteger(0);
+    result.ival = len;
+    result.words = xwords;
+    return result.canonicalize();
+  }
+
+  private void setInvert()
+  {
+    if (words == null)
+      ival = ~ival;
+    else
+      {
+	for (int i = ival;  --i >= 0; )
+	  words[i] = ~words[i];
+      }
+  }
+
+  public BigInteger not()
+  {
+    BigInteger result = new BigInteger();
+    result.ival = ival;
+    result.words = words;
+    result.setInvert();
+    return result;
+  }
+
+  private void setShiftLeft(BigInteger x, int count)
+  {
+    int[] xwords;
+    int xlen;
+    if (x.words == null)
+      {
+	if (count < 32)
+	  {
+	    set((long) x.ival << count);
+	    return;
+	  }
+	xwords = new int[1];
+	xwords[0] = x.ival;
+	xlen = 1;
+      }
+    else
+      {
+	xwords = x.words;
+	xlen = x.ival;
+      }
+    int word_count = count >> 5;
+    count &= 31;
+    int new_len = xlen + word_count;
+    if (count == 0)
+      {
+	realloc(new_len);
+	for (int i = xlen;  --i >= 0; )
+	  words[i+word_count] = xwords[i];
+      }
+    else
+      {
+	new_len++;
+	realloc(new_len);
+	int shift_out = MPN.lshift(words, word_count, xwords, xlen, count);
+	count = 32 - count;
+	words[new_len-1] = (shift_out << count) >> count;  // sign-extend.
+      }
+    ival = new_len;
+    for (int i = word_count;  --i >= 0; )
+      words[i] = 0;
+  }
+
+  private void setShiftRight(BigInteger x, int count)
+  {
+    if (x.words == null)
+      set(count < 32 ? x.ival >> count : x.ival < 0 ? -1 : 0);
+    else if (count == 0)
+      set(x);
+    else
+      {
+	boolean neg = x.isNegative();
+	int word_count = count >> 5;
+	count &= 31;
+	int d_len = x.ival - word_count;
+	if (d_len <= 0)
+	  set(neg ? -1 : 0);
+	else
+	  {
+	    if (words == null || words.length < d_len)
+	      realloc(d_len);
+	    MPN.rshift(words, x.words, word_count, d_len, count);
+	    ival = d_len;
+	    if (neg)
+	      words[ival-1] |= -1 << (32 - count);
+	  }
+      }
+  }
+
+  private void setShift(BigInteger x, int count)
+  {
+    if (count > 0)
+      setShiftLeft(x, count);
+    else
+      setShiftRight(x, -count);
+  }
+
+  private static BigInteger shift(BigInteger x, int count)
+  {
+    if (x.words == null)
+      {
+	if (count <= 0)
+	  return make(count > -32 ? x.ival >> (-count) : x.ival < 0 ? -1 : 0);
+	if (count < 32)
+	  return make((long) x.ival << count);
+      }
+    if (count == 0)
+      return x;
+    BigInteger result = new BigInteger(0);
+    result.setShift(x, count);
+    return result.canonicalize();
+  }
+
+  public BigInteger shiftLeft(int n)
+  {
+    return shift(this, n);
+  }
+
+  public BigInteger shiftRight(int n)
+  {
+    return shift(this, -n);
+  }
+
+  private void format(int radix, StringBuffer buffer)
+  {
+    if (words == null)
+      buffer.append(Integer.toString(ival, radix));
+    else if (ival <= 2)
+      buffer.append(Long.toString(longValue(), radix));
+    else
+      {
+	boolean neg = isNegative();
+	int[] work;
+	if (neg || radix != 16)
+	  {
+	    work = new int[ival];
+	    getAbsolute(work);
+	  }
+	else
+	  work = words;
+	int len = ival;
+
+	int buf_size = len * (MPN.chars_per_word(radix) + 1);
+	if (radix == 16)
+	  {
+	    if (neg)
+	      buffer.append('-');
+	    int buf_start = buffer.length();
+	    for (int i = len;  --i >= 0; )
+	      {
+		int word = work[i];
+		for (int j = 8;  --j >= 0; )
+		  {
+		    int hex_digit = (word >> (4 * j)) & 0xF;
+		    // Suppress leading zeros:
+		    if (hex_digit > 0 || buffer.length() > buf_start)
+		      buffer.append(Character.forDigit(hex_digit, 16));
+		  }
+	      }
+	  }
+	else
+	  {
+	    int i = buffer.length();
+	    for (;;)
+	      {
+		int digit = MPN.divmod_1(work, work, len, radix);
+		buffer.append(Character.forDigit(digit, radix));
+		while (len > 0 && work[len-1] == 0) len--;
+		if (len == 0)
+		  break;
+	      }
+	    if (neg)
+	      buffer.append('-');
+	    /* Reverse buffer. */
+	    int j = buffer.length() - 1;
+	    while (i < j)
+	      {
+		char tmp = buffer.charAt(i);
+		buffer.setCharAt(i, buffer.charAt(j));
+		buffer.setCharAt(j, tmp);
+		i++;  j--;
+	      }
+	  }
+      }
+  }
+
+  public String toString()
+  {
+    return toString(10);
+  }
+
+  public String toString(int radix)
+  {
+    if (words == null)
+      return Integer.toString(ival, radix);
+    else if (ival <= 2)
+      return Long.toString(longValue(), radix);
+    int buf_size = ival * (MPN.chars_per_word(radix) + 1);
+    StringBuffer buffer = new StringBuffer(buf_size);
+    format(radix, buffer);
+    return buffer.toString();
+  }
+
+  public int intValue()
+  {
+    if (words == null)
+      return ival;
+    return words[0];
+  }
+
+  public long longValue()
+  {
+    if (words == null)
+      return ival;
+    if (ival == 1)
+      return words[0];
+    return ((long)words[1] << 32) + ((long)words[0] & 0xffffffffL);
+  }
+
+  public int hashCode()
+  {
+    // FIXME: May not match hashcode of JDK.
+    return words == null ? ival : (words[0] + words[ival - 1]);
+  }
+
+  /* Assumes x and y are both canonicalized. */
+  private static boolean equals(BigInteger x, BigInteger y)
+  {
+    if (x.words == null && y.words == null)
+      return x.ival == y.ival;
+    if (x.words == null || y.words == null || x.ival != y.ival)
+      return false;
+    for (int i = x.ival; --i >= 0; )
+      {
+	if (x.words[i] != y.words[i])
+	  return false;
+      }
+    return true;
+  }
+
+  /* Assumes this and obj are both canonicalized. */
+  public boolean equals(Object obj)
+  {
+    if (obj == null || ! (obj instanceof BigInteger))
+      return false;
+    return BigInteger.equals(this, (BigInteger) obj);
+  }
+
+  public double doubleValue()
+  {
+    if (words == null)
+      return (double) ival;
+    if (ival <= 2)
+      return (double) longValue();
+    if (isNegative())
+      return BigInteger.neg(this).roundToDouble(0, true, false);
+    else
+      return roundToDouble(0, false, false);
+  }
+
+  public float floatValue()
+  {
+    return (float) doubleValue();
+  }
+
+  /** Return true if any of the lowest n bits are one.
+   * (false if n is negative).  */
+  private boolean checkBits(int n)
+  {
+    if (n <= 0)
+      return false;
+    if (words == null)
+      return n > 31 || ((ival & ((1 << n) - 1)) != 0);
+    int i;
+    for (i = 0; i < (n >> 5) ; i++)
+      if (words[i] != 0)
+	return true;
+    return (n & 31) != 0 && (words[i] & ((1 << (n & 31)) - 1)) != 0;
+  }
+
+  /** Convert a semi-processed BigInteger to double.
+   * Number must be non-negative.  Multiplies by a power of two, applies sign,
+   * and converts to double, with the usual java rounding.
+   * @param exp power of two, positive or negative, by which to multiply
+   * @param neg true if negative
+   * @param remainder true if the BigInteger is the result of a truncating
+   * division that had non-zero remainder.  To ensure proper rounding in
+   * this case, the BigInteger must have at least 54 bits.  */
+  private double roundToDouble(int exp, boolean neg, boolean remainder)
+  {
+    // Compute length.
+    int il = bitLength();
+
+    // Exponent when normalized to have decimal point directly after
+    // leading one.  This is stored excess 1023 in the exponent bit field.
+    exp += il - 1;
+
+    // Gross underflow.  If exp == -1075, we let the rounding
+    // computation determine whether it is minval or 0 (which are just
+    // 0x0000 0000 0000 0001 and 0x0000 0000 0000 0000 as bit
+    // patterns).
+    if (exp < -1075)
+      return neg ? -0.0 : 0.0;
+
+    // gross overflow
+    if (exp > 1023)
+      return neg ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
+
+    // number of bits in mantissa, including the leading one.
+    // 53 unless it's denormalized
+    int ml = (exp >= -1022 ? 53 : 53 + exp + 1022);
+
+    // Get top ml + 1 bits.  The extra one is for rounding.
+    long m;
+    int excess_bits = il - (ml + 1);
+    if (excess_bits > 0)
+      m = ((words == null) ? ival >> excess_bits
+	   : MPN.rshift_long(words, ival, excess_bits));
+    else
+      m = longValue() << (- excess_bits);
+
+    // Special rounding for maxval.  If the number exceeds maxval by
+    // any amount, even if it's less than half a step, it overflows.
+    if (exp == 1023 && ((m >> 1) == (1L << 53) - 1))
+      {
+	if (remainder || checkBits(il - ml))
+	  return neg ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
+	else
+	  return neg ? - Double.MAX_VALUE : Double.MAX_VALUE;
+      }
+
+    // Normal round-to-even rule: round up if the bit dropped is a one, and
+    // the bit above it or any of the bits below it is a one.
+    if ((m & 1) == 1
+	&& ((m & 2) == 2 || remainder || checkBits(excess_bits)))
+      {
+	m += 2;
+	// Check if we overflowed the mantissa
+	if ((m & (1L << 54)) != 0)
+	  {
+	    exp++;
+	    // renormalize
+	    m >>= 1;
+	  }
+	// Check if a denormalized mantissa was just rounded up to a
+	// normalized one.
+	else if (ml == 52 && (m & (1L << 53)) != 0)
+	  exp++;
+      }
+	
+    // Discard the rounding bit
+    m >>= 1;
+
+    long bits_sign = neg ? (1L << 63) : 0;
+    exp += 1023;
+    long bits_exp = (exp <= 0) ? 0 : ((long)exp) << 52;
+    long bits_mant = m & ~(1L << 52);
+    return Double.longBitsToDouble(bits_sign | bits_exp | bits_mant);
+  }
+
+  /** Copy the abolute value of this into an array of words.
+   * Assumes words.length >= (this.words == null ? 1 : this.ival).
+   * Result is zero-extended, but need not be a valid 2's complement number.
+   */
+    
+  private void getAbsolute(int[] words)
+  {
+    int len;
+    if (this.words == null)
+      {
+	len = 1;
+	words[0] = this.ival;
+      }
+    else
+      {
+	len = this.ival;
+	for (int i = len;  --i >= 0; )
+	  words[i] = this.words[i];
+      }
+    if (words[len - 1] < 0)
+      negate(words, words, len);
+    for (int i = words.length;  --i > len; )
+      words[i] = 0;
+  }
+
+  /** Set dest[0:len-1] to the negation of src[0:len-1].
+   * Return true if overflow (i.e. if src is -2**(32*len-1)).
+   * Ok for src==dest. */
+  private static boolean negate(int[] dest, int[] src, int len)
+  {
+    long carry = 1;
+    boolean negative = src[len-1] < 0;
+    for (int i = 0;  i < len;  i++)
+      {
+        carry += ((long) (~src[i]) & 0xffffffffL);
+        dest[i] = (int) carry;
+        carry >>= 32;
+      }
+    return (negative && dest[len-1] < 0);
+  }
+
+  /** Destructively set this to the negative of x.
+   * It is OK if x==this.*/
+  private void setNegative(BigInteger x)
+  {
+    int len = x.ival;
+    if (x.words == null)
+      {
+	if (len == Integer.MIN_VALUE)
+	  set(- (long) len);
+	else
+	  set(-len);
+	return;
+      }
+    realloc(len + 1);
+    if (BigInteger.negate(words, x.words, len))
+      words[len++] = 0;
+    ival = len;
+  }
+
+  /** Destructively negate this. */
+  private final void setNegative()
+  {
+    setNegative(this);
+  }
+
+  private static BigInteger abs(BigInteger x)
+  {
+    return x.isNegative() ? neg(x) : x;
+  }
+
+  public BigInteger abs()
+  {
+    return abs(this);
+  }
+
+  public static BigInteger neg(BigInteger x)
+  {
+    if (x.words == null && x.ival != Integer.MIN_VALUE)
+      return make(- x.ival);
+    BigInteger result = new BigInteger(0);
+    result.setNegative(x);
+    return result.canonicalize();
+  }
+
+  public BigInteger negate()
+  {
+    return BigInteger.neg(this);
+  }
+
+  /** Calculates ceiling(log2(this < 0 ? -this : this+1))
+   * See Common Lisp: the Language, 2nd ed, p. 361.
+   */
+  public int bitLength()
+  {
+    if (words == null)
+      return MPN.intLength(ival);
+    else
+      return MPN.intLength(words, ival);
+  }
+
+/* TODO:
+
+  public BigInteger(String val, int radix)
+
+  public BigInteger(String val)
+
+  public BigInteger(int bitLength, int certainty, Random rnd)
+
+  public BigInteger and(BigInteger val)
+
+  public BigInteger or(BigInteger val)
+
+  public BigInteger xor(BigInteger val
+
+  public BigInteger andNot(BigInteger val)
+
+  public BigInteger clearBit(int n)
+  {
+    if (n < 0)
+      throw new ArithmeticException();
+
+    return and(ONE.shiftLeft(n).not());
+  }
+
+  public BigInteger setBit(int n)
+  {
+    if (n < 0)
+      throw new ArithmeticException();
+
+    return or(ONE.shiftLeft(n));
+  }
+
+  public boolean testBit(int n)
+
+  public BigInteger flipBit(int n)
+
+  public int getLowestSetBit()
+
+  public int bitCount()
+
+  public boolean isProbablePrime(int certainty)
+
+  public BigInteger min(BigInteger val)
+
+  public BigInteger max(BigInteger val)
+
+  public byte[] toByteArray()
+*/
+}