Commits

didier deshommes committed 0e0fe32

adding MPN.java, wrap GMP in GMPInteger.java

Comments (0)

Files changed (3)

+/* java.math.BigInteger -- Arbitary precision integers
+   Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007  Free Software Foundation, Inc.
+
+This file is part of GNU Classpath.
+
+GNU Classpath is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+ 
+GNU Classpath is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU Classpath; see the file COPYING.  If not, write to the
+Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301 USA.
+
+Linking this library statically or dynamically with other modules is
+making a combined work based on this library.  Thus, the terms and
+conditions of the GNU General Public License cover the whole
+combination.
+
+As a special exception, the copyright holders of this library give you
+permission to link this library with independent modules to produce an
+executable, regardless of the license terms of these independent
+modules, and to copy and distribute the resulting executable under
+terms of your choice, provided that you also meet, for each linked
+independent module, the terms and conditions of the license of that
+module.  An independent module is a module which is not derived from
+or based on this library.  If you modify this library, you may extend
+this exception to your version of the library, but you are not
+obligated to do so.  If you do not wish to do so, delete this
+exception statement from your version. */
+
+
+//package java.math;
+
+//import gnu.classpath.Configuration;
+
+//import gnu.java.lang.CPStringBuilder;
+//import GMP;
+//import MPN;
+
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
+import java.util.Random;
+import java.util.logging.Logger;
+
+/**
+ * Written using on-line Java Platform 1.2 API Specification, as well
+ * as "The Java Class Libraries", 2nd edition (Addison-Wesley, 1998) and
+ * "Applied Cryptography, Second Edition" by Bruce Schneier (Wiley, 1996).
+ * 
+ * Based primarily on IntNum.java BitOps.java by Per Bothner (per@bothner.com)
+ * (found in Kawa 1.6.62).
+ *
+ * @author Warren Levy (warrenl@cygnus.com)
+ * @date December 20, 1999.
+ * @status believed complete and correct.
+ */
+public class GMPInteger extends Number implements Comparable<GMPInteger>
+{
+  private static final Logger log = Logger.getLogger(GMPInteger.class.getName());
+
+  /** All integers are stored in 2's-complement form.
+   * If words == null, the ival is the value of this GMPInteger.
+   * Otherwise, the first ival elements of words make the value
+   * of this GMPInteger, stored in little-endian order, 2's-complement form. */
+  private transient int ival;
+  private transient int[] words;
+
+  // Serialization fields.
+  // the first three, although not used in the code, are present for
+  // compatibility with older RI versions of this class. DO NOT REMOVE.
+  private int bitCount = -1;
+  private int bitLength = -1;
+  private int lowestSetBit = -2;
+  private byte[] magnitude;
+  private int signum;
+  private static final long serialVersionUID = -8287574255936472291L;
+
+
+  /** We pre-allocate integers in the range minFixNum..maxFixNum. 
+   * Note that we must at least preallocate 0, 1, and 10.  */
+  private static final int minFixNum = -100;
+  private static final int maxFixNum = 1024;
+  private static final int numFixNum = maxFixNum-minFixNum+1;
+  private static final GMPInteger[] smallFixNums;
+  
+  /** The alter-ego GMP instance for this. */
+  private transient GMP mpz;
+
+  static
+  {
+    
+      smallFixNums = null;
+      ZERO = valueOf(0L);
+      ONE = valueOf(1L);
+      TEN = valueOf(10L);
+        
+  }
+
+  /**
+   * The constant zero as a GMPInteger.
+   * @since 1.2
+   */
+  public static final GMPInteger ZERO;
+
+  /**
+   * The constant one as a GMPInteger.
+   * @since 1.2
+   */
+  public static final GMPInteger ONE;
+
+  /**
+   * The constant ten as a GMPInteger.
+   * @since 1.5
+   */
+  public static final GMPInteger TEN;
+
+  /* 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;
+
+  /** When checking the probability of primes, it is most efficient to
+   * first check the factoring of small primes, so we'll use this array.
+   */
+  private static final int[] primes =
+    {   2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,
+       47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97, 101, 103, 107,
+      109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
+      191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251 };
+
+  /** HAC (Handbook of Applied Cryptography), Alfred Menezes & al. Table 4.4. */
+  private static final int[] k =
+      {100,150,200,250,300,350,400,500,600,800,1250, Integer.MAX_VALUE};
+  private static final int[] t =
+      { 27, 18, 15, 12,  9,  8,  7,  6,  5,  4,   3, 2};
+
+  private GMPInteger()
+  {
+    super();
+    mpz = new GMP();
+  }
+
+  /* Create a new (non-shared) GMPInteger, and initialize to an int. */
+  private GMPInteger(int value)
+  {
+    super();
+
+    ival = value;
+  }
+
+  public GMPInteger(String s, int radix)
+  {
+    this();
+
+    int len = s.length();
+    int i, digit;
+    boolean negative;
+    byte[] bytes;
+    char ch = s.charAt(0);
+    if (ch == '-')
+      {
+        negative = true;
+        i = 1;
+        bytes = new byte[len - 1];
+      }
+    else
+      {
+        negative = false;
+        i = 0;
+        bytes = new byte[len];
+      }
+    int byte_len = 0;
+    for ( ; i < len;  i++)
+      {
+        ch = s.charAt(i);
+        digit = Character.digit(ch, radix);
+        if (digit < 0)
+          throw new NumberFormatException("Invalid character at position #" + i);
+        bytes[byte_len++] = (byte) digit;
+      }
+
+    bytes = null;
+    if (mpz.fromString(s, radix) != 0)
+	throw new NumberFormatException("String \"" + s
+					+ "\" is NOT a valid number in base "
+					+ radix);
+         
+  }
+
+  public GMPInteger(String val)
+  {
+    this(val, 10);
+  }
+
+  /* Create a new (non-shared) GMPInteger, and initialize from a byte array. */
+  public GMPInteger(byte[] val)
+  {
+    this();
+
+    if (val == null || val.length < 1)
+      throw new NumberFormatException();
+    
+    mpz.fromByteArray(val);
+    
+  }
+
+  public GMPInteger(int signum, byte[] magnitude)
+  {
+    this();
+
+    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;
+      }
+    
+    mpz.fromSignedMagnitude(magnitude, signum == -1);
+    
+  }
+
+  public GMPInteger(int numBits, Random rnd)
+  {
+    this();
+
+    if (numBits < 0)
+      throw new IllegalArgumentException();
+
+    init(numBits, rnd);
+  }
+
+  private void init(int numBits, Random rnd)
+  {
+
+      int length = (numBits + 7) / 8;
+      byte[] magnitude = new byte[length];
+      rnd.nextBytes(magnitude);
+      int discardedBitCount = numBits % 8;
+      if (discardedBitCount != 0)
+          {
+	      discardedBitCount = 8 - discardedBitCount;
+	      magnitude[0] = (byte)((magnitude[0] & 0xFF) >>> discardedBitCount);
+          }
+      mpz.fromSignedMagnitude(magnitude, false);
+      magnitude = null;
+      return;
+      
+  }
+
+  public GMPInteger(int bitLength, int certainty, Random rnd)
+  {
+    this();
+
+    GMPInteger result = new GMPInteger();
+    while (true)
+      {
+        result.init(bitLength, rnd);
+        result = result.setBit(bitLength - 1);
+        if (result.isProbablePrime(certainty))
+          break;
+      }
+
+    mpz.fromBI(result.mpz);
+    
+  }
+
+  /** 
+   *  Return a GMPInteger that is bitLength bits long with a
+   *  probability < 2^-100 of being composite.
+   *
+   *  @param bitLength length in bits of resulting number
+   *  @param rnd random number generator to use
+   *  @throws ArithmeticException if bitLength < 2
+   *  @since 1.4
+   */
+  public static GMPInteger probablePrime(int bitLength, Random rnd)
+  {
+    if (bitLength < 2)
+      throw new ArithmeticException();
+
+    return new GMPInteger(bitLength, 100, rnd);
+  }
+
+  /** Return a (possibly-shared) GMPInteger with a given long value. */
+  public static GMPInteger valueOf(long val)
+  {
+      GMPInteger result = new GMPInteger();
+      result.mpz.fromLong(val);
+      return result;
+          
+  }
+
+  /**
+   * @return <code>true</code> if the GMP-based native implementation library
+   *         was successfully loaded. Returns <code>false</code> otherwise.
+   */
+  private static boolean initializeLibrary()
+  {
+    boolean result;
+    try
+    {
+      System.loadLibrary("javamath");
+      GMP.natInitializeLibrary();
+      result = true;
+    }
+    catch (Throwable x)
+    {
+      result = false;
+
+    }
+    return result;
+  }
+
+  /** Make a canonicalized GMPInteger from an array of words.
+   * The array may be reused (without copying). */
+  private static GMPInteger make(int[] words, int len)
+  {
+    if (words == null)
+      return valueOf(len);
+    len = GMPInteger.wordsNeeded(words, len);
+    if (len <= 1)
+      return len == 0 ? ZERO : valueOf(words[0]);
+    GMPInteger num = new GMPInteger();
+    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/4 + 1];
+    int nwords = words.length;
+
+    // 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) | (bytes[bptr] & 0xff);
+    words[--nwords] = word;
+
+    // Elements remaining in byte[] are a multiple of 4.
+    while (nwords > 0)
+      words[--nwords] = bytes[bptr++] << 24 |
+			(bytes[bptr++] & 0xff) << 16 |
+			(bytes[bptr++] & 0xff) << 8 |
+			(bytes[bptr++] & 0xff);
+    return words;
+  }
+
+  /** Allocate a new non-shared GMPInteger.
+   * @param nwords number of words to allocate
+   */
+  private static GMPInteger alloc(int nwords)
+  {
+    GMPInteger result = new GMPInteger();
+    if (nwords > 1)
+    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 boolean isNegative()
+  {
+    return (words == null ? ival : words[ival - 1]) < 0;
+  }
+
+  public int signum()
+  {
+      return mpz.compare(ZERO.mpz);
+  }
+
+  private static int compareTo(GMPInteger x, GMPInteger y)
+  {
+      int dummy = y.signum; // force NPE check
+      return x.mpz.compare(y.mpz);
+  }
+
+  /** @since 1.2 */
+  public int compareTo(GMPInteger val)
+  {
+    return compareTo(this, val);
+  }
+
+  public GMPInteger min(GMPInteger val)
+  {
+    return compareTo(this, val) < 0 ? this : val;
+  }
+
+  public GMPInteger max(GMPInteger val)
+  {
+    return compareTo(this, val) > 0 ? this : val;
+  }
+
+  private boolean isZero()
+  {
+    return words == null && ival == 0;
+  }
+
+  private boolean isOne()
+  {
+    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 GMPInteger canonicalize()
+  {
+    if (words != null
+	&& (ival = GMPInteger.wordsNeeded(words, ival)) <= 1)
+      {
+	if (ival == 1)
+	  ival = words[0];
+	words = null;
+      }
+    if (words == null && ival >= minFixNum && ival <= maxFixNum)
+      return smallFixNums[ival - minFixNum];
+    return this;
+  }
+
+  /** Add two ints, yielding a GMPInteger. */
+  private static GMPInteger add(int x, int y)
+  {
+    return valueOf((long) x + (long) y);
+  }
+
+  /** Add a GMPInteger and an int, yielding a new GMPInteger. */
+  private static GMPInteger add(GMPInteger x, int y)
+  {
+    if (x.words == null)
+      return GMPInteger.add(x.ival, y);
+    GMPInteger result = new GMPInteger(0);
+    result.setAdd(x, y);
+    return result.canonicalize();
+  }
+
+  /** Set this to the sum of x and y.
+   * OK if x==this. */
+  private void setAdd(GMPInteger 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 void setAdd(int y)
+  {
+    setAdd(this, y);
+  }
+
+  /** Destructively set the value of this to a long. */
+  private 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 void set(int[] words, int length)
+  {
+    this.ival = length;
+    this.words = words;
+  }
+
+  /** Destructively set the value of this to that of y. */
+  private void set(GMPInteger 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 GMPIntegers, yielding their sum as another GMPInteger. */
+  private static GMPInteger add(GMPInteger x, GMPInteger y, int k)
+  {
+    if (x.words == null && y.words == null)
+      return valueOf((long) k * (long) y.ival + (long) x.ival);
+    if (k != 1)
+      {
+	if (k == -1)
+	  y = GMPInteger.neg(y);
+	else
+	  y = GMPInteger.times(y, valueOf(k));
+      }
+    if (x.words == null)
+      return GMPInteger.add(y, x.ival);
+    if (y.words == null)
+      return GMPInteger.add(x, y.ival);
+    // Both are big
+    if (y.ival > x.ival)
+      { // Swap so x is longer then y.
+	GMPInteger tmp = x;  x = y;  y = tmp;
+      }
+    GMPInteger 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 GMPInteger add(GMPInteger val)
+  {
+      int dummy = val.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.add(val.mpz, result.mpz);
+      return result;
+      
+  }
+
+  public GMPInteger subtract(GMPInteger val)
+  {
+      int dummy = val.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.subtract(val.mpz, result.mpz);
+      return result;
+      
+  }
+
+  private static GMPInteger times(GMPInteger 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 valueOf((long) xlen * (long) y);
+    boolean negative;
+    GMPInteger result = GMPInteger.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 GMPInteger times(GMPInteger x, GMPInteger 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;
+      }
+    GMPInteger result = GMPInteger.alloc(xlen+ylen);
+    MPN.mul(result.words, xwords, xlen, ywords, ylen);
+    result.ival = xlen+ylen;
+    if (negative)
+      result.setNegative();
+    return result.canonicalize();
+  }
+
+  public GMPInteger multiply(GMPInteger y)
+  {
+      int dummy = y.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.multiply(y.mpz, result.mpz);
+      return result;
+
+  }
+
+  private static void divide(long x, long y,
+			     GMPInteger quotient, GMPInteger remainder,
+			     int rounding_mode)
+  {
+    boolean xNegative, yNegative;
+    if (x < 0)
+      {
+	xNegative = true;
+	if (x == Long.MIN_VALUE)
+	  {
+	    divide(valueOf(x), valueOf(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(valueOf(x), valueOf(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(GMPInteger x, GMPInteger y,
+			     GMPInteger quotient, GMPInteger 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;
+	// Need to leave room for a word of leading zeros if dividing by 1
+	// and the dividend has the high bit set.  It might be safe to
+	// increment qlen in all cases, but it certainly is only necessary
+	// in the following case.
+	if (ywords[0] == 1 && xwords[xlen-1] < 0)
+	  qlen++;
+	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;
+	MPN.rshift0 (ywords, xwords, 0, rlen, nshift);
+
+	qlen = xlen + 1 - ylen;
+	if (quotient != null)
+	  {
+	    for (int i = 0;  i < qlen;  i++)
+	      xwords[i] = xwords[i+ylen];
+	  }
+      }
+
+    if (ywords[rlen-1] < 0)
+      {
+        ywords[rlen] = 0;
+        rlen++;
+      }
+
+    // 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));
+	    GMPInteger tmp = remainder == null ? new GMPInteger() : 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)).
+	    GMPInteger tmp;
+	    if (y.words == null)
+	      {
+		tmp = remainder;
+		tmp.set(yNegative ? ywords[0] + y.ival : ywords[0] - y.ival);
+	      }
+	    else
+	      tmp = GMPInteger.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 GMPInteger divide(GMPInteger val)
+  {
+      if (val.compareTo(ZERO) == 0)
+          throw new ArithmeticException("divisor is zero");
+
+      GMPInteger result = new GMPInteger();
+      mpz.quotient(val.mpz, result.mpz);
+      return result;
+  }
+
+  public GMPInteger remainder(GMPInteger val)
+  {
+      if (val.compareTo(ZERO) == 0)
+          throw new ArithmeticException("divisor is zero");
+      
+      GMPInteger result = new GMPInteger();
+      mpz.remainder(val.mpz, result.mpz);
+      return result;
+      
+  }
+
+  public GMPInteger[] divideAndRemainder(GMPInteger val)
+  {
+      if (val.compareTo(ZERO) == 0)
+          throw new ArithmeticException("divisor is zero");
+
+      GMPInteger q = new GMPInteger();
+      GMPInteger r = new GMPInteger();
+      mpz.quotientAndRemainder(val.mpz, q.mpz, r.mpz);
+      return new GMPInteger[] { q, r };
+  }
+
+  public GMPInteger mod(GMPInteger m)
+  {
+      int dummy = m.signum; // force NPE check
+      if (m.compareTo(ZERO) < 1)
+          throw new ArithmeticException("non-positive modulus");
+
+      GMPInteger result = new GMPInteger();
+      mpz.modulo(m.mpz, result.mpz);
+      return result;
+  
+  }
+
+  /** Calculate the integral power of a GMPInteger.
+   * @param exponent the exponent (must be non-negative)
+   */
+  public GMPInteger pow(int exponent)
+  {
+    if (exponent <= 0)
+      {
+	if (exponent == 0)
+	  return ONE;
+	  throw new ArithmeticException("negative exponent");
+      }
+
+   
+    GMPInteger result = new GMPInteger();
+    mpz.pow(exponent, result.mpz);
+    return result;
+    
+  }
+
+  private static int[] euclidInv(int a, int b, int prevDiv)
+  {
+    if (b == 0)
+      throw new ArithmeticException("not invertible");
+
+    if (b == 1)
+	// Success:  values are indeed invertible!
+	// Bottom of the recursion reached; start unwinding.
+	return new int[] { -prevDiv, 1 };
+
+    int[] xy = euclidInv(b, a % b, a / b);	// Recursion happens here.
+    a = xy[0]; // use our local copy of 'a' as a work var
+    xy[0] = a * -prevDiv + xy[1];
+    xy[1] = a;
+    return xy;
+  }
+
+  private static void euclidInv(GMPInteger a, GMPInteger b,
+                                GMPInteger prevDiv, GMPInteger[] xy)
+  {
+    if (b.isZero())
+      throw new ArithmeticException("not invertible");
+
+    if (b.isOne())
+      {
+	// Success:  values are indeed invertible!
+	// Bottom of the recursion reached; start unwinding.
+	xy[0] = neg(prevDiv);
+        xy[1] = ONE;
+	return;
+      }
+
+    // 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[0] = new GMPInteger(xyInt[0]);
+        xy[1] = new GMPInteger(xyInt[1]);
+      }
+    else
+      {
+	GMPInteger rem = new GMPInteger();
+	GMPInteger quot = new GMPInteger();
+	divide(a, b, quot, rem, FLOOR);
+        // quot and rem may not be in canonical form. ensure
+        rem.canonicalize();
+        quot.canonicalize();
+	euclidInv(b, rem, quot, xy);
+      }
+
+    GMPInteger t = xy[0];
+    xy[0] = add(xy[1], times(t, prevDiv), -1);
+    xy[1] = t;
+  }
+
+  public GMPInteger modInverse(GMPInteger y)
+  {
+      int dummy = y.signum; // force NPE check
+      if (mpz.compare(ZERO.mpz) < 1)
+          throw new ArithmeticException("non-positive modulo");
+
+      GMPInteger result = new GMPInteger();
+      mpz.modInverse(y.mpz, result.mpz);
+      return result;
+          
+  }
+
+  public GMPInteger modPow(GMPInteger exponent, GMPInteger m)
+  {
+      int dummy = exponent.signum; // force NPE check
+      if (m.mpz.compare(ZERO.mpz) < 1)
+          throw new ArithmeticException("non-positive modulo");
+      
+      GMPInteger result = new GMPInteger();
+      mpz.modPow(exponent.mpz, m.mpz, result.mpz);
+      return result;
+
+  }
+
+  /** Calculate Greatest Common Divisor for non-negative ints. */
+  private static int gcd(int a, int b)
+  {
+    // Euclid's algorithm, copied from libg++.
+    int tmp;
+    if (b > a)
+      {
+	tmp = a; a = b; b = tmp;
+      }
+    for(;;)
+      {
+	if (b == 0)
+	  return a;
+        if (b == 1)
+	  return b;
+        tmp = b;
+	    b = a % b;
+	    a = tmp;
+	  }
+      }
+
+  public GMPInteger gcd(GMPInteger y)
+  {
+      int dummy = y.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.gcd(y.mpz, result.mpz);
+      return result;
+      
+  }
+
+  /**
+   * <p>Returns <code>true</code> if this GMPInteger is probably prime,
+   * <code>false</code> if it's definitely composite. If <code>certainty</code>
+   * is <code><= 0</code>, <code>true</code> is returned.</p>
+   *
+   * @param certainty a measure of the uncertainty that the caller is willing
+   * to tolerate: if the call returns <code>true</code> the probability that
+   * this GMPInteger is prime exceeds <code>(1 - 1/2<sup>certainty</sup>)</code>.
+   * The execution time of this method is proportional to the value of this
+   * parameter.
+   * @return <code>true</code> if this GMPInteger is probably prime,
+   * <code>false</code> if it's definitely composite.
+   */
+  public boolean isProbablePrime(int certainty)
+  {
+    if (certainty < 1)
+      return true;
+    
+    return mpz.testPrimality(certainty) != 0;
+
+  }
+
+  private void setInvert()
+  {
+    if (words == null)
+      ival = ~ival;
+    else
+      {
+	for (int i = ival;  --i >= 0; )
+	  words[i] = ~words[i];
+      }
+  }
+
+  private void setShiftLeft(GMPInteger 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(GMPInteger 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.rshift0 (words, x.words, word_count, d_len, count);
+	    ival = d_len;
+	    if (neg)
+	      words[d_len-1] |= -2 << (31 - count);
+	  }
+      }
+  }
+
+  private void setShift(GMPInteger x, int count)
+  {
+    if (count > 0)
+      setShiftLeft(x, count);
+    else
+      setShiftRight(x, -count);
+  }
+
+  private static GMPInteger shift(GMPInteger x, int count)
+  {
+    if (x.words == null)
+      {
+	if (count <= 0)
+	  return valueOf(count > -32 ? x.ival >> (-count) : x.ival < 0 ? -1 : 0);
+	if (count < 32)
+	  return valueOf((long) x.ival << count);
+      }
+    if (count == 0)
+      return x;
+    GMPInteger result = new GMPInteger(0);
+    result.setShift(x, count);
+    return result.canonicalize();
+  }
+
+  public GMPInteger shiftLeft(int n)
+  {
+    if (n == 0)
+      return this;
+
+
+    GMPInteger result = new GMPInteger();
+    if (n < 0)
+	mpz.shiftRight(-n, result.mpz);
+    else
+	mpz.shiftLeft(n, result.mpz);
+    return result;
+    
+  }
+
+  public GMPInteger shiftRight(int n)
+  {
+    if (n == 0)
+      return this;
+
+    GMPInteger result = new GMPInteger();
+    if (n < 0)
+	mpz.shiftLeft(-n, result.mpz);
+    else
+	mpz.shiftRight(n, result.mpz);
+    return result;
+
+  }
+
+  private void format(int radix, StringBuilder 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;
+
+	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)
+  {
+      return mpz.toString(radix);
+  }
+
+  public int intValue()
+  {
+      int result = mpz.absIntValue();
+      return mpz.compare(ZERO.mpz) < 0 ? - result : result;
+
+  }
+
+  public long longValue()
+  {
+      long result;
+      result = (abs().shiftRight(32)).mpz.absIntValue();
+      result <<= 32;
+      result |= mpz.absIntValue() & 0xFFFFFFFFL;
+      return this.compareTo(ZERO) < 0 ? - result : result;
+      
+  }
+
+  public int hashCode()
+  {
+    // FIXME: May not match hashcode of JDK.
+    
+      // TODO: profile to decide whether to make it native
+      byte[] bytes = this.toByteArray();
+      int result = 0;
+      for (int i = 0; i < bytes.length; i++)
+          result ^= (bytes[i] & 0xFF) << (8 * (i % 4));
+      return result;
+      
+  }
+
+  /* Assumes x and y are both canonicalized. */
+  private static boolean equals(GMPInteger x, GMPInteger y)
+  {
+      return x.mpz.compare(y.mpz) == 0;
+  }
+
+  /* Assumes this and obj are both canonicalized. */
+  public boolean equals(Object obj)
+  {
+    if (! (obj instanceof GMPInteger))
+      return false;
+    return equals(this, (GMPInteger) obj);
+  }
+
+  private static GMPInteger valueOf(byte[] digits, int byte_len,
+				    boolean negative, int radix)
+  {
+    int chars_per_word = MPN.chars_per_word(radix);
+    int[] words = new int[byte_len / chars_per_word + 1];
+    int size = MPN.set_str(words, digits, byte_len, radix);
+    if (size == 0)
+      return ZERO;
+    if (words[size-1] < 0)
+      words[size++] = 0;
+    if (negative)
+      negate(words, words, size);
+    return make(words, size);
+  }
+
+  public double doubleValue()
+  {
+      return mpz.doubleValue();
+  }
+
+  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 GMPInteger 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 GMPInteger is the result of a truncating
+   * division that had non-zero remainder.  To ensure proper rounding in
+   * this case, the GMPInteger 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(GMPInteger 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 (negate(words, x.words, len))
+      words[len++] = 0;
+    ival = len;
+  }
+
+  /** Destructively negate this. */
+  private void setNegative()
+  {
+    setNegative(this);
+  }
+
+  private static GMPInteger abs(GMPInteger x)
+  {
+    return x.isNegative() ? neg(x) : x;
+  }
+
+  public GMPInteger abs()
+  {
+      GMPInteger result = new GMPInteger();
+      mpz.abs(result.mpz);
+      return result;
+  }
+
+  private static GMPInteger neg(GMPInteger x)
+  {
+    if (x.words == null && x.ival != Integer.MIN_VALUE)
+      return valueOf(- x.ival);
+    GMPInteger result = new GMPInteger(0);
+    result.setNegative(x);
+    return result.canonicalize();
+  }
+
+  public GMPInteger negate()
+  {
+      GMPInteger result = new GMPInteger();
+      mpz.negate(result.mpz);
+      return result;
+  }
+
+  /** Calculates ceiling(log2(this < 0 ? -this : this+1))
+   * See Common Lisp: the Language, 2nd ed, p. 361.
+   */
+  public int bitLength()
+  {
+      return mpz.bitLength();
+  }
+
+  public byte[] toByteArray()
+  {
+    if (signum() == 0)
+      return new byte[1];
+
+    // the minimal number of bytes required to represent the MPI is function
+    // of (a) its bit-length, and (b) its sign.  only when this MPI is both
+    // positive, and its bit-length is a multiple of 8 do we add one zero
+    // bit for its sign.  we do this so if we construct a new MPI from the
+    // resulting byte array, we wouldn't mistake a positive number, whose
+    // bit-length is a multiple of 8, for a similar-length negative one.
+    int bits = bitLength();
+    if (bits % 8 == 0 || this.signum() == 1)
+	bits++;
+    byte[] bytes = new byte[(bits + 7) / 8];
+    mpz.toByteArray(bytes);
+    return bytes;
+
+  }
+
+  /** Return the boolean opcode (for bitOp) for swapped operands.
+   * I.e. bitOp(swappedOp(op), x, y) == bitOp(op, y, x).
+   */
+  private static int swappedOp(int op)
+  {
+    return
+    "\000\001\004\005\002\003\006\007\010\011\014\015\012\013\016\017"
+    .charAt(op);
+  }
+
+  /** Do one the the 16 possible bit-wise operations of two GMPIntegers. */
+  private static GMPInteger bitOp(int op, GMPInteger x, GMPInteger y)
+  {
+    switch (op)
+      {
+        case 0:  return ZERO;
+        case 1:  return x.and(y);
+        case 3:  return x;
+        case 5:  return y;
+        case 15: return valueOf(-1);
+      }
+    GMPInteger result = new GMPInteger();
+    setBitOp(result, op, x, y);
+    return result.canonicalize();
+  }
+
+  /** Do one the the 16 possible bit-wise operations of two GMPIntegers. */
+  private static void setBitOp(GMPInteger result, int op,
+			       GMPInteger x, GMPInteger y)
+  {
+    if ((y.words != null) && (x.words == null || x.ival < y.ival))
+      {
+	GMPInteger temp = x;  x = y;  y = temp;
+	op = swappedOp(op);
+      }
+    int xi;
+    int yi;
+    int xlen, ylen;
+    if (y.words == null)
+      {
+	yi = y.ival;
+	ylen = 1;
+      }
+    else
+      {
+	yi = y.words[0];
+	ylen = y.ival;
+      }
+    if (x.words == null)
+      {
+	xi = x.ival;
+	xlen = 1;
+      }
+    else
+      {
+	xi = x.words[0];
+	xlen = x.ival;
+      }
+    if (xlen > 1)
+      result.realloc(xlen);
+    int[] w = result.words;
+    int i = 0;
+    // Code for how to handle the remainder of x.
+    // 0:  Truncate to length of y.
+    // 1:  Copy rest of x.
+    // 2:  Invert rest of x.
+    int finish = 0;
+    int ni;
+    switch (op)
+      {
+      case 0:  // clr
+	ni = 0;
+	break;
+      case 1: // and
+	for (;;)
+	  {
+	    ni = xi & yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi < 0) finish = 1;
+	break;
+      case 2: // andc2
+	for (;;)
+	  {
+	    ni = xi & ~yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi >= 0) finish = 1;
+	break;
+      case 3:  // copy x
+	ni = xi;
+	finish = 1;  // Copy rest
+	break;
+      case 4: // andc1
+	for (;;)
+	  {
+	    ni = ~xi & yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi < 0) finish = 2;
+	break;
+      case 5: // copy y
+	for (;;)
+	  {
+	    ni = yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	break;
+      case 6:  // xor
+	for (;;)
+	  {
+	    ni = xi ^ yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	finish = yi < 0 ? 2 : 1;
+	break;
+      case 7:  // ior
+	for (;;)
+	  {
+	    ni = xi | yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi >= 0) finish = 1;
+	break;
+      case 8:  // nor
+	for (;;)
+	  {
+	    ni = ~(xi | yi);
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi >= 0)  finish = 2;
+	break;
+      case 9:  // eqv [exclusive nor]
+	for (;;)
+	  {
+	    ni = ~(xi ^ yi);
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	finish = yi >= 0 ? 2 : 1;
+	break;
+      case 10:  // c2
+	for (;;)
+	  {
+	    ni = ~yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	break;
+      case 11:  // orc2
+	for (;;)
+	  {
+	    ni = xi | ~yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi < 0)  finish = 1;
+	break;
+      case 12:  // c1
+	ni = ~xi;
+	finish = 2;
+	break;
+      case 13:  // orc1
+	for (;;)
+	  {
+	    ni = ~xi | yi;
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi >= 0) finish = 2;
+	break;
+      case 14:  // nand
+	for (;;)
+	  {
+	    ni = ~(xi & yi);
+	    if (i+1 >= ylen) break;
+	    w[i++] = ni;  xi = x.words[i];  yi = y.words[i];
+	  }
+	if (yi < 0) finish = 2;
+	break;
+      default:
+      case 15:  // set
+	ni = -1;
+	break;
+      }
+    // Here i==ylen-1; w[0]..w[i-1] have the correct result;
+    // and ni contains the correct result for w[i+1].
+    if (i+1 == xlen)
+      finish = 0;
+    switch (finish)
+      {
+      case 0:
+	if (i == 0 && w == null)
+	  {
+	    result.ival = ni;
+	    return;
+	  }
+	w[i++] = ni;
+	break;
+      case 1:  w[i] = ni;  while (++i < xlen)  w[i] = x.words[i];  break;
+      case 2:  w[i] = ni;  while (++i < xlen)  w[i] = ~x.words[i];  break;
+      }
+    result.ival = i;
+  }
+
+  /** Return the logical (bit-wise) "and" of a GMPInteger and an int. */
+  private static GMPInteger and(GMPInteger x, int y)
+  {
+    if (x.words == null)
+      return valueOf(x.ival & y);
+    if (y >= 0)
+      return valueOf(x.words[0] & y);
+    int len = x.ival;
+    int[] words = new int[len];
+    words[0] = x.words[0] & y;
+    while (--len > 0)
+      words[len] = x.words[len];
+    return make(words, x.ival);
+  }
+
+  /** Return the logical (bit-wise) "and" of two GMPIntegers. */
+  public GMPInteger and(GMPInteger y)
+  {
+      int dummy = y.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.and(y.mpz, result.mpz);
+      return result;
+      
+  }
+
+  /** Return the logical (bit-wise) "(inclusive) or" of two GMPIntegers. */
+  public GMPInteger or(GMPInteger y)
+  {
+      int dummy = y.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.or(y.mpz, result.mpz);
+      return result;
+      
+  }
+
+  /** Return the logical (bit-wise) "exclusive or" of two GMPIntegers. */
+  public GMPInteger xor(GMPInteger y)
+  {
+
+      int dummy = y.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.xor(y.mpz, result.mpz);
+      return result;
+
+  }
+
+  /** Return the logical (bit-wise) negation of a GMPInteger. */
+  public GMPInteger not()
+  {
+      GMPInteger result = new GMPInteger();
+      mpz.not(result.mpz);
+      return result;
+
+  }
+
+  public GMPInteger andNot(GMPInteger val)
+  {
+      int dummy = val.signum; // force NPE check
+      GMPInteger result = new GMPInteger();
+      mpz.andNot(val.mpz, result.mpz);
+      return result;
+      
+  }
+
+  public GMPInteger clearBit(int n)
+  {
+    if (n < 0)
+      throw new ArithmeticException();
+
+    GMPInteger result = new GMPInteger();
+    mpz.setBit(n, false, result.mpz);
+    return result;
+	
+  }
+
+  public GMPInteger setBit(int n)
+  {
+    if (n < 0)
+      throw new ArithmeticException();
+    
+    GMPInteger result = new GMPInteger();
+    mpz.setBit(n, true, result.mpz);
+    return result;
+    
+  }
+
+  public boolean testBit(int n)
+  {
+    if (n < 0)
+      throw new ArithmeticException();
+
+    return mpz.testBit(n) != 0;
+
+  }
+
+  public GMPInteger flipBit(int n)
+  {
+    if (n < 0)
+      throw new ArithmeticException();
+    
+    GMPInteger result = new GMPInteger();
+    mpz.flipBit(n, result.mpz);
+    return result;
+    
+  }
+
+  public int getLowestSetBit()
+  {
+      return mpz.compare(ZERO.mpz) == 0 ? -1 : mpz.lowestSetBit();
+  }
+
+  // bit4count[I] is number of '1' bits in I.
+  private static final byte[] bit4_count = { 0, 1, 1, 2,  1, 2, 2, 3,
+					     1, 2, 2, 3,  2, 3, 3, 4};
+
+  private static int bitCount(int i)
+  {
+    int count = 0;
+    while (i != 0)
+      {
+	count += bit4_count[i & 15];
+	i >>>= 4;
+      }
+    return count;
+  }
+
+  private static int bitCount(int[] x, int len)
+  {
+    int count = 0;
+    while (--len >= 0)
+      count += bitCount(x[len]);
+    return count;
+  }
+
+  /** Count one bits in a GMPInteger.
+   * If argument is negative, count zero bits instead. */
+  public int bitCount()
+  {
+      return mpz.bitCount();
+  }
+
+  private void readObject(ObjectInputStream s)
+    throws IOException, ClassNotFoundException
+  {
+      mpz = new GMP();
+      s.defaultReadObject();
+      if (signum != 0)
+          mpz.fromByteArray(magnitude);
+        
+  }
+
+  private void writeObject(ObjectOutputStream s)
+    throws IOException, ClassNotFoundException
+  {
+    signum = signum();
+    magnitude = signum == 0 ? new byte[0] : toByteArray();
+    s.defaultWriteObject();
+    magnitude = null; // not needed anymore
+  }
+
+  // inner class(es) ..........................................................
+
+}
+/* gnu.java.math.MPN
+   Copyright (C) 1999, 2000, 2001, 2004  Free Software Foundation, Inc.
+
+This file is part of GNU Classpath.
+
+GNU Classpath is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+ 
+GNU Classpath is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU Classpath; see the file COPYING.  If not, write to the
+Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301 USA.
+
+Linking this library statically or dynamically with other modules is
+making a combined work based on this library.  Thus, the terms and
+conditions of the GNU General Public License cover the whole
+combination.
+
+As a special exception, the copyright holders of this library give you
+permission to link this library with independent modules to produce an
+executable, regardless of the license terms of these independent
+modules, and to copy and distribute the resulting executable under
+terms of your choice, provided that you also meet, for each linked
+independent module, the terms and conditions of the license of that
+module.  An independent module is a module which is not derived from
+or based on this library.  If you modify this library, you may extend
+this exception to your version of the library, but you are not
+obligated to do so.  If you do not wish to do so, delete this
+exception statement from your version. */
+
+// Included from Kawa 1.6.62 with permission of the author,
+// Per Bothner <per@bothner.com>.
+
+//package gnu.java.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));