/* java.math.BigInteger -- Arbitary precision integers
Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.
The file was changed by Radek Polak to work as midlet in MIDP 1.0
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., 59 Temple Place, Suite 330, Boston, MA
02111-1307 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 ssh.v1;
/**
* @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) 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).
*
* Status: Believed complete and correct.
*/
public class BigInteger {
/**
* 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.
*/
transient private int ival;
transient private 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: RADEK: just floor */
// private static final int FLOOR = 1;
private BigInteger() {
}
/* Create a new (non-shared) BigInteger, and initialize to an int. */
private BigInteger( int value ) {
ival = value;
}
// RADEK: cunstruct new BigInteger assuming that signum = 1
public BigInteger( byte[] magnitude ) {
// 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;
}
/** Return a (possibly-shared) BigInteger with a given long value. */
public static BigInteger valueOf( long val ) {
if ( val >= minFixNum && val <= maxFixNum )
return smallFixNums[(int) val - minFixNum];
int i = (int) val;
if ( (long) i == val )
return new BigInteger( i );
BigInteger result = alloc( 2 );
result.ival = 2;
result.words[0] = i;
result.words[1] = (int) ( val >> 32 );
return result;
}
/**
* 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 valueOf( len );
len = BigInteger.wordsNeeded( words, len );
if ( len <= 1 )
return len == 0 ? ZERO : valueOf( 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 / 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 BigInteger.
*
* @param nwords
* number of words to allocate
*/
private static BigInteger alloc( int nwords ) {
BigInteger result = new BigInteger();
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 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;
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 ) ? 1 : -1;
return cmp( x.words, y.words, x_len );
}
public int compareTo( BigInteger val ) {
return compareTo( this, val );
}
private final boolean isZero() {
return words == null && ival == 0;
}
private final 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 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[ival - minFixNum];
return this;
}
// /** Add two ints, yielding a BigInteger. */
// private static final BigInteger add( int x, int y ) {
// return valueOf( (long) x + (long) y );
// }
// /** Add a 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 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 valueOf( (long) k * (long) y.ival + (long) x.ival );
// if ( k != 1 )
// y = BigInteger.times( y, valueOf( k ) );
// if ( x.words == null )
// return BigInteger.add( y, x.ival );
// if ( y.words == null )
// return BigInteger.add( x, y.ival );
// // Both are big
// 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 = 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();
// }
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 valueOf( (long) xlen * (long) y );
BigInteger result = BigInteger.alloc( xlen + 1 );
result.words[xlen] = mul_1( result.words, xwords, xlen, y );
result.ival = xlen + 1;
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 );
int[] xwords;
int[] ywords;
int xlen = x.ival;
int ylen = y.ival;
xwords = x.words;
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 );
mul( result.words, xwords, xlen, ywords, ylen );
result.ival = xlen + ylen;
return result.canonicalize();
}
private static void divide( long x, long y, BigInteger quotient, BigInteger remainder ) {
boolean xNegative, yNegative;
if ( x < 0 ) {
xNegative = true;
if ( x == Long.MIN_VALUE ) {
divide( valueOf( x ), valueOf( y ), quotient, remainder );
return;
}
x = -x;
}
else
xNegative = false;
if ( y < 0 ) {
yNegative = true;
if ( y == Long.MIN_VALUE ) {
divide( valueOf( x ), valueOf( y ), quotient, remainder );
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 ) {
if ( qNegative )
add_one = true;
}
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 ) {
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 );
return;
}
}
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 = 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] = 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 = count_leading_zeros( ywords[ylen - 1] );
if ( nshift != 0 ) {
// Shift up the denominator setting the most significant bit of
// the most significant word.
lshift( ywords, 0, ywords, ylen, nshift );
// Shift up the numerator, possibly introducing a new most
// significant word.
int x_high = lshift( xwords, 0, xwords, xlen, nshift );
xwords[xlen++] = x_high;
}
if ( xlen == ylen )
xwords[xlen++] = 0;
divide( xwords, xlen, ywords, ylen );
rlen = ylen;
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.
if ( quotient != null )
quotient.set( xwords, qlen );
if ( remainder != null )
// The remainder is by definition: X-Q*Y
remainder.set( ywords, rlen );
}
public BigInteger mod( BigInteger m ) {
BigInteger rem = new BigInteger();
divide( this, m, null, rem );
return rem.canonicalize();
}
public BigInteger modPow( BigInteger exponent, BigInteger 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);
u = valueOf( u.ival >> 1 );
t = times( t, t ).mod( m );
}
return s;
}
public long longValue() {
if ( words == null )
return ival;
if ( ival == 1 )
return words[0];
return ( (long) words[1] << 32 ) + ( (long) words[0] & 0xffffffffL );
}
/**
* 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];
}
for ( int i = words.length; --i > len; )
words[i] = 0;
}
/**
* Calculates ceiling(log2(this < 0 ? -this : this+1)) See Common Lisp: the
* Language, 2nd ed, p. 361.
*/
public int bitLength() {
if ( words == null )
return intLength( ival );
return intLength( words, ival );
}
public byte[] toByteArray() {
// Determine number of bytes needed. The method bitlength returns
// the size without the sign bit, so add one bit for that and then
// add 7 more to emulate the ceil function using integer math.
byte[] bytes = new byte[( bitLength() + 1 + 7 ) / 8];
int nbytes = bytes.length;
int wptr = 0;
int word;
// Deal with words array until one word or less is left to process.
// If BigInteger is an int, then it is in ival and nbytes will be <= 4.
while ( nbytes > 4 ) {
word = words[wptr++];
for ( int i = 4; i > 0; --i, word >>= 8 )
bytes[--nbytes] = (byte) word;
}
// Deal with the last few bytes. If BigInteger is an int, use ival.
word = ( words == null ) ? ival : words[wptr];
for ( ; nbytes > 0; word >>= 8 )
bytes[--nbytes] = (byte) word;
return bytes;
}
/** Return the logical (bit-wise) "and" of a BigInteger and an int. */
private static BigInteger and( BigInteger 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 BigIntegers. */
public BigInteger and( BigInteger y ) {
if ( y.words == null )
return and( this, y.ival );
else if ( words == null )
return and( y, ival );
BigInteger x = this;
if ( ival < y.ival ) {
BigInteger temp = this;
x = y;
y = temp;
}
int i;
int len = y.ival;
int[] words = new int[len];
for ( i = 0; i < y.ival; i++ )
words[i] = x.words[i] & y.words[i];
for ( ; i < len; i++ )
words[i] = x.words[i];
return make( words, len );
}
/**
* 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] = 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 = 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 = mul_1( dest, dest, size, big_base );
cy_limb += 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;
}
/*
* 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.
* OK if dest==x. Assumes: 0 <= count < 32 Same as rshift, but handles
* count==0 (and has no return value).
*/
public static void rshift0( int[] dest, int[] x, int x_start, int len, int count ) {
if ( count > 0 )
rshift( dest, x, x_start, len, count );
else
for ( int i = 0; i < len; i++ )
dest[i] = x[i + x_start];
}
/**
* 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] 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. */
public 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;
}
/**
* 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;
rshift0( x, x, initShiftWords, len, initShiftBits );
rshift0( 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 )
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 = cmp( odd_arg, other_arg, len );
if ( i == 0 )
break;
if ( i > 0 ) { // odd_arg > other_arg
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
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 = 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==BigInteger.wordsNeeded(words,len)
*/
public static int intLength( int[] words, int len ) {
len--;
return intLength( words[len] ) + 32 * len;
}
}