biginteger.java

俄罗斯高人Mamaich的Pocket gcc编译器（运行在PocketPC上）的全部源代码。

	// original modulus, y.ival (not the possibly "swapped" yval).
	if (result.ival < 0)
	  result.ival += y.ival;
      }
    else
      {
	// As above, force this to be a positive value via modulo math.
	BigInteger x = isNegative() ? this.mod(y) : this;

	// Swap values so x > y.
	if (x.compareTo(y) < 0)
	  {
	    result = x; x = y; y = result; // use 'result' as a work var
	    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);
        // quot and rem may not be in canonical form. ensure
        rem.canonicalize();
        quot.canonicalize();
	BigInteger[] xy = new BigInteger[2];
	euclidInv(y, rem, quot, xy);
	result = swapped ? xy[0] : xy[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);

    // To do this naively by first raising this to the power of exponent
    // and then performing modulo m would be extremely expensive, especially
    // for very large numbers.  The solution is found in Number Theory
    // where a combination of partial powers and moduli can be done easily.
    //
    // We'll use the algorithm for Additive Chaining which can be found on
    // p. 244 of "Applied Cryptography, Second Edition" by Bruce Schneier.
    BigInteger s = ONE;
    BigInteger t = this;
    BigInteger u = exponent;

    while (!u.isZero())
      {
	if (u.and(ONE).isOne())
	  s = times(s, t).mod(m);
	u = u.shiftRight(1);
	t = times(t, t).mod(m);
      }

    return s;
  }

  /** Calculate Greatest Common Divisor for non-negative ints. */
  private static final 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 BigInteger gcd(BigInteger y)
  {
    int xval = ival;
    int yval = y.ival;
    if (words == null)
      {
	if (xval == 0)
	  return 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 valueOf(gcd(xval, yval));
	  }
	xval = 1;
      }
    if (y.words == null)
      {
	if (yval == 0)
	  return 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();
  }

  /**
   * <p>Returns <code>true</code> if this BigInteger 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 BigInteger 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 BigInteger is probably prime,
   * <code>false</code> if it's definitely composite.
   */
  public boolean isProbablePrime(int certainty)
  {
    if (certainty < 1)
      return true;

    /** We'll use the Rabin-Miller algorithm for doing a probabilistic
     * primality test.  It is fast, easy and has faster decreasing odds of a
     * composite passing than with other tests.  This means that this
     * method will actually have a probability much greater than the
     * 1 - .5^certainty specified in the JCL (p. 117), but I don't think
     * anyone will complain about better performance with greater certainty.
     *
     * The Rabin-Miller algorithm can be found on pp. 259-261 of "Applied
     * Cryptography, Second Edition" by Bruce Schneier.
     */

    // First rule out small prime factors
    BigInteger rem = new BigInteger();
    int i;
    for (i = 0; i < primes.length; i++)
      {
	if (words == null && ival == primes[i])
	  return true;

        divide(this, smallFixNums[primes[i] - minFixNum], null, rem, TRUNCATE);
        if (rem.canonicalize().isZero())
	  return false;
      }

    // Now perform the Rabin-Miller test.

    // Set b to the number of times 2 evenly divides (this - 1).
    // I.e. 2^b is the largest power of 2 that divides (this - 1).
    BigInteger pMinus1 = add(this, -1);
    int b = pMinus1.getLowestSetBit();

    // Set m such that this = 1 + 2^b * m.
    BigInteger m = pMinus1.divide(valueOf(2L << b - 1));

    // The HAC (Handbook of Applied Cryptography), Alfred Menezes & al. Note
    // 4.49 (controlling the error probability) gives the number of trials
    // for an error probability of 1/2**80, given the number of bits in the
    // number to test.  we shall use these numbers as is if/when 'certainty'
    // is less or equal to 80, and twice as much if it's greater.
    int bits = this.bitLength();
    for (i = 0; i < k.length; i++)
      if (bits <= k[i])
        break;
    int trials = t[i];
    if (certainty > 80)
      trials *= 2;
    BigInteger z;
    for (int t = 0; t < trials; t++)
      {
        // The HAC (Handbook of Applied Cryptography), Alfred Menezes & al.
        // Remark 4.28 states: "...A strategy that is sometimes employed
        // is to fix the bases a to be the first few primes instead of
        // choosing them at random.
	z = smallFixNums[primes[t] - minFixNum].modPow(m, this);
	if (z.isOne() || z.equals(pMinus1))
	  continue;			// Passes the test; may be prime.

	for (i = 0; i < b; )
	  {
	    if (z.isOne())
	      return false;
	    i++;
	    if (z.equals(pMinus1))
	      break;			// Passes the test; may be prime.

	    z = z.modPow(valueOf(2), this);
	  }

	if (i == b && !z.equals(pMinus1))
	  return false;
      }
    return true;
  }

  private void setInvert()
  {
    if (words == null)
      ival = ~ival;
    else
      {
	for (int i = ival;  --i >= 0; )
	  words[i] = ~words[i];
      }
  }

  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.rshift0 (words, x.words, word_count, d_len, count);
	    ival = d_len;
	    if (neg)
	      words[d_len-1] |= -2 << (31 - 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 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;
    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;

	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);
    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 equals(this, (BigInteger) obj);
  }

  private static BigInteger valueOf(String s, int radix)
       throws NumberFormatException
  {
    int len = s.length();
    // Testing (len < MPN.chars_per_word(radix)) would be more accurate,
    // but slightly more expensive, for little practical gain.
    if (len <= 15 && radix <= 16)
      return valueOf(Long.parseLong(s, radix));
    
    int byte_len = 0;
    byte[] bytes = new byte[len];
    boolean negative = false;
    for (int i = 0;  i < len;  i++)
      {
	char ch = s.charAt(i);
	if (ch == '-')
	  negative = true;
	else if (ch == '_' || (byte_len == 0 && (ch == ' ' || ch == '\t')))
	  continue;
	else
	  {
	    int digit = Character.digit(ch, radix);
	    if (digit < 0)
	      break;
	    bytes[byte_len++] = (byte) digit;
	  }
      }
    return valueOf(bytes, byte_len, negative, radix);
  }

  private static BigInteger 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()
  {
    if (words == null)
      return (double) ival;
    if (ival <= 2)
      return (double) longValue();
    if (isNegative())
      return neg(this).roundToDouble(0, true, false);
      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;