You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@harmony.apache.org by ml...@apache.org on 2006/11/13 07:03:28 UTC

svn commit: r474164 - in /incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math: BigInteger.java BitLevel.java Division.java Elementary.java

Author: mloenko
Date: Sun Nov 12 22:03:27 2006
New Revision: 474164

URL: http://svn.apache.org/viewvc?view=rev&rev=474164
Log:
applied newest patch for HARMONY-2091
[classlib][java.math] optimization of BigInteger.modInverse

Modified:
    incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BigInteger.java
    incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BitLevel.java
    incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Division.java
    incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Elementary.java

Modified: incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BigInteger.java
URL: http://svn.apache.org/viewvc/incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BigInteger.java?view=diff&rev=474164&r1=474163&r2=474164
==============================================================================
--- incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BigInteger.java (original)
+++ incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BigInteger.java Sun Nov 12 22:03:27 2006
@@ -569,34 +569,23 @@
         BigInteger val1 = this.abs();
         BigInteger val2 = val.abs();
         // To avoid a possible division by zero
-        if (sign == 0) {
+        if (val1.signum() == 0) {
             return val2;
-        }
-        if (val.sign == 0) {
+        } else if (val2.signum() == 0) {
             return val1;
         }
-        /* Implements one step of the Euclidean algorithm
-         * To reduce some operand if it's much smaller than the other one */
-        if (val1.numberLength > val2.numberLength) {
-            val1 = val1.remainder(val2);
-            if (val1.sign == 0) {
-                return val2;
-            }
-        } else if (val2.numberLength > val1.numberLength) {
-            val2 = val2.remainder(val1);
-            if (val2.sign == 0) {
-                return val1;
-            }
-        }
-        /* Optimization for small operands
-         * (val1.bitLength() < 64)   and   (val2.bitLength() < 64) */
-        if (((val1.numberLength == 1) || ((val1.numberLength == 2) && (val1.digits[1] > 0)))
-                && ((val2.numberLength == 1) || ((val2.numberLength == 2) && (val2.digits[1] > 0)))) {
-            return BigInteger.valueOf(Division.gcdBinary(val1.longValue(), val2
-                    .longValue()));
+        	
+        // Optimization for small operands
+        // (op2.bitLength() < 64) and (op1.bitLength() < 64)
+    	if ((( val1.numberLength == 1 )
+    			|| ( ( val1.numberLength == 2 ) && ( val1.digits[1] > 0 ) ))
+            && (val2.numberLength == 1 || (val2.numberLength == 2 && val2.digits[1] > 0) )) {
+                return BigInteger.valueOf(Division.gcdBinary(val1.longValue(),
+                        val2.longValue()));
         }
-        // Now 'val1' and 'val2' will be muttable
+    	
         return Division.gcdBinary(val1.copy(), val2.copy());
+    
     }
 
     /** @ar.org.fitc.spec_ref */
@@ -740,20 +729,24 @@
             throw new ArithmeticException(Messages.getString("math.18")); //$NON-NLS-1$
         }
         // If both are even, no inverse exists
-        if (!(testBit(0) || m.testBit(0))) {
+        if (!( testBit(0) || m.testBit(0) )) {
             // math.19=BigInteger not invertible.
             throw new ArithmeticException(Messages.getString("math.19")); //$NON-NLS-1$
         }
         if (m.isOne()) {
             return ZERO;
         }
+        
         // From now on: (m > 1)
-        BigInteger res = Division.modInverse(abs(), m);
+        BigInteger res = Division.modInverseMontgomery(abs().mod(m), m);
         if (res.sign == 0) {
             // math.19=BigInteger not invertible.
             throw new ArithmeticException(Messages.getString("math.19")); //$NON-NLS-1$
         }
-        return ((sign < 0) ? m.subtract(res) : res);
+        
+        res = ( ( sign < 0 ) ? m.subtract(res) : res );
+        return res;
+        
     }
 
     /** @ar.org.fitc.spec_ref */
@@ -941,4 +934,8 @@
         out.defaultWriteObject();
     }
 
+    
+    void unCache(){
+    	firstNonzeroDigit = -2;
+    }
 }

Modified: incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BitLevel.java
URL: http://svn.apache.org/viewvc/incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BitLevel.java?view=diff&rev=474164&r1=474163&r2=474164
==============================================================================
--- incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BitLevel.java (original)
+++ incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/BitLevel.java Sun Nov 12 22:03:27 2006
@@ -115,8 +115,9 @@
     /** @see BigInteger#shiftLeft(int) */
     static BigInteger shiftLeft(BigInteger source, int count) {
         int intCount = count >> 5;
-        count &= 31;
-        int resLength = source.numberLength + intCount + ((count == 0) ? 0 : 1);
+        count &= 31; // %= 32
+        int resLength = source.numberLength + intCount
+                + ( ( count == 0 ) ? 0 : 1 );
         int resDigits[] = new int[resLength];
 
         shiftLeft(resDigits, source.digits, intCount, count);
@@ -126,8 +127,23 @@
     }
 
     /**
-     * Shifts left an array of integers. Total shift distance in bits is
-     * intCount * 32 + count
+     * Performs {@code val <<= count}.
+     */
+    // val should have enough place (and one digit more)
+    static void inplaceShiftLeft(BigInteger val, int count) {
+        int intCount = count >> 5; // count of integers
+        val.numberLength += intCount
+                + ( Integer
+                .numberOfLeadingZeros(val.digits[val.numberLength - 1])
+                - ( count & 31 ) >= 0 ? 0 : 1 );
+        shiftLeft(val.digits, val.digits, intCount, count & 31);
+        val.cutOffLeadingZeroes();
+        val.unCache();
+    }
+    
+    /**
+     * Abstractly shifts left an array of integers in little endian (i.e. shift
+     * it right). Total shift distance in bits is intCount * 32 + count
      * 
      * @param result the destination array
      * @param source the source array
@@ -147,6 +163,10 @@
                 result[i - 1] = source[i - intCount - 1] << count;
             }
         }
+        
+        for (int i = 0; i < intCount; i++) {
+            result[i] = 0;
+        }
     }
 
     /** @see BigInteger#shiftRight(int) */
@@ -188,40 +208,69 @@
      * Performs {@code val >>= count} where {@code val} is a positive number.
      */
     static void inplaceShiftRight(BigInteger val, int count) {
-        // PRE: val > 0
+        int sign = val.signum();
+        if (count == 0 || val.signum() == 0)
+            return;
         int intCount = count >> 5; // count of integers
         val.numberLength -= intCount;
-        shiftRight(val.digits, val.numberLength, val.digits, intCount,
-                count & 31);
+        if (!shiftRight(val.digits, val.numberLength, val.digits, intCount,
+                count & 31)
+                && sign < 0) {
+            // remainder not zero: add one to the result
+            int i;
+            for (i = 0; ( i < val.numberLength ) && ( val.digits[i] == -1 ); i++) {
+                val.digits[i] = 0;
+            }
+            if (i == val.numberLength) {
+                val.numberLength++;
+            }
+            val.digits[i]++;
+        }
         val.cutOffLeadingZeroes();
+        val.unCache();
     }
 
     /**
      * Shifts right an array of integers. Total shift distance in bits is
      * intCount * 32 + count.
      * 
-     * @param result the destination array
-     * @param resultLen the destination array's length
-     * @param source the source array
-     * @param intCount the number of elements to be shifted
-     * @param count the number of bits to be shifted
+     * @param result
+     *            the destination array
+     * @param resultLen
+     *            the destination array's length
+     * @param source
+     *            the source array
+     * @param intCount
+     *            the number of elements to be shifted
+     * @param count
+     *            the number of bits to be shifted
+     * @return dropped bit's are all zero (i.e. remaider is zero)
      */
-    static void shiftRight(int result[], int resultLen, int source[],
+    static boolean shiftRight(int result[], int resultLen, int source[],
             int intCount, int count) {
+        int i;
+        boolean allZero = true;
+        for (i = 0; i < intCount; i++)
+            allZero &= source[i] == 0;
         if (count == 0) {
             System.arraycopy(source, intCount, result, 0, resultLen);
+            i = resultLen;
         } else {
-            int i;
             int leftShiftCount = 32 - count;
 
+            allZero &= ( source[i] << leftShiftCount ) == 0;
             for (i = 0; i < resultLen - 1; i++) {
-                result[i] = (source[i + intCount] >>> count)
-                        | (source[i + intCount + 1] << leftShiftCount);
+                result[i] = ( source[i + intCount] >>> count )
+                | ( source[i + intCount + 1] << leftShiftCount );
             }
-            result[i] = (source[i + intCount] >>> count);
+            result[i] = ( source[i + intCount] >>> count );
+            i++;
         }
+        
+        return allZero;
     }
 
+    
     /**
      * Performs a flipBit on the BigInteger, returning a BigInteger with the the
      * specified bit flipped.

Modified: incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Division.java
URL: http://svn.apache.org/viewvc/incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Division.java?view=diff&rev=474164&r1=474163&r2=474164
==============================================================================
--- incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Division.java (original)
+++ incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Division.java Sun Nov 12 22:03:27 2006
@@ -17,6 +17,8 @@
 
 package java.math;
 
+import org.apache.harmony.math.internal.nls.Messages;
+
 /**
  * Static library that provides all operations related with division and modular
  * arithmetic to {@link BigInteger}. Some methos are provided in both mutable
@@ -377,47 +379,68 @@
     }
 
     /**
-     * Return the greatest common divisor of op1 and op2, based in The Binary
-     * GCD Algorithm of D. Knuth.
+     * @param m a positive modulus
+     * Return the greatest common divisor of op1 and op2,
      * 
-     * @param op1 must be greater than zero
-     * @param op2 must be greater than zero
+     * @param op1
+     *            must be greater than zero
+     * @param op2
+     *            must be greater than zero
      * @see BigInteger#gcd(BigInteger)
-     * @ar.org.fitc.ref "The Art of Computer Programming Vo.2, Section 4.5.2,
-     *                  Algorithm B by D. Knuth."
      * @return {@code GCD(op1, op2)}
      */
     static BigInteger gcdBinary(BigInteger op1, BigInteger op2) {
         // PRE: (op1 > 0) and (op2 > 0)
+        
         /*
          * Divide both number the maximal possible times by 2 without rounding
-         * gcd(2*a, 2*b) = 2 * gcd(a,n)
+                 * gcd(2*a, 2*b) = 2 * gcd(a,b)
          */
         int lsb1 = op1.getLowestSetBit();
         int lsb2 = op2.getLowestSetBit();
         int pow2Count = Math.min(lsb1, lsb2);
 
-        if (lsb1 != 0) {
             BitLevel.inplaceShiftRight(op1, lsb1);
-        }
-        if (lsb2 != 0) {
             BitLevel.inplaceShiftRight(op2, lsb2);
-        }
-        do {// gcd(2*k+1,2*n) = gcd(2*k+1, n)
-            lsb1 = op1.getLowestSetBit();
-            if (lsb1 != 0) {
-                BitLevel.inplaceShiftRight(op1, lsb1);
-            }
-            // gcd(2*n,2*k+1) = gcd(n,2*k+1)
-            lsb2 = op2.getLowestSetBit();
-            if (lsb2 != 0) {
-                BitLevel.inplaceShiftRight(op2, lsb2);
+        
+        BigInteger swap;
+        // I want op2 > op1
+        if (op1.compareTo(op2) == BigInteger.GREATER) {
+            swap = op1;
+            op1 = op2;
+            op2 = swap;
+        } 
+        
+        do { // INV: op2 >= op1 && both are odd unless op1 = 0
+            
+            // Optimization for small operands
+            // (op2.bitLength() < 64) implies by INV (op1.bitLength() < 64)
+            if (( op2.numberLength == 1 )
+            || ( ( op2.numberLength == 2 ) && ( op2.digits[1] > 0 ) )) {
+                op2 = BigInteger.valueOf(Division.gcdBinary(op1.longValue(),
+                        op2.longValue()));
+                break;
+        }
+            
+            // Implements one step of the Euclidean algorithm
+            // To reduce one operand if it's much smaller than the other one
+            if (op2.numberLength > op1.numberLength * 1.2) {
+                op2 = op2.remainder(op1);
+                if (op2.signum() != 0) {
+                    BitLevel.inplaceShiftRight(op2, op2.getLowestSetBit());
             }
-            if (op1.compareTo(op2) >= BigInteger.EQUALS) {
-                Elementary.inplaceSubtract(op1, op2);
             } else {
-                Elementary.inplaceSubtract(op2, op1);
-            }
+                
+                // Use Knuth's algorithm of sucessive subtract and shifting
+                do {
+                    Elementary.inplaceSubtract(op2, op1); // both are odd
+                    BitLevel.inplaceShiftRight(op2, op2.getLowestSetBit()); // op2 is even
+                } while (op2.compareTo(op1) >= BigInteger.EQUALS);
+            }
+            // now op1 >= op2
+            swap = op2;
+            op2 = op1;
+            op1 = swap;
         } while (op1.sign != 0);
         return op2.shiftLeft(pow2Count);
     }
@@ -427,8 +450,10 @@
      * with numbers of 63 bits, represented in positives values of {@code long}
      * type.
      * 
-     * @param op1 a postive number
-     * @param op2 a postive number
+     * @param op1
+     *            a postive number
+     * @param op2
+     *            a postive number
      * @see #gcdBinary(BigInteger, BigInteger)
      * @return <code>GCD(op1, op2)</code>
      */
@@ -445,99 +470,290 @@
             op2 >>>= lsb2;
         }
         do {
-            lsb1 = Long.numberOfTrailingZeros(op1);
-            if (lsb1 != 0) {
-                op1 >>>= lsb1;
-            }
-            lsb2 = Long.numberOfTrailingZeros(op2);
-            if (lsb2 != 0) {
-                op2 >>>= lsb2;
-            }
             if (op1 >= op2) {
                 op1 -= op2;
+                op1 >>>= Long.numberOfTrailingZeros(op1);
             } else {
                 op2 -= op1;
+                op2 >>>= Long.numberOfTrailingZeros(op2);
             }
         } while (op1 != 0);
-        return (op2 << pow2Count);
+        return ( op2 << pow2Count );
     }
 
+    
+    
+    
     /**
-     * Implements the "Shifting Euclidean modular inverse algorithm".
-     * 
-     * @see BigInteger#modInverse(BigInteger)
-     * @param a a positive number
-     * @param m a positive modulus
-     * @ar.org.fitc.ref "Laszlo Hars - Modular Inverse Algorithms Without
-     *                  Multiplications for Cryptographic Applications"
+     * Calculates a.modInverse(p) Based on: Savas, E; Koc, C "The Montgomery Modular
+     * Inverse - Revisted"
+     */
+    static BigInteger modInverseMontgomery(BigInteger a, BigInteger p) {
+
+        if (a.sign == 0){
+            // ZERO hasn't inverse
+            // math.19: BigInteger not invertible
+            throw new ArithmeticException(Messages.getString("math.19"));
+            }
+        
+        
+        if (!p.testBit(0)){
+            // montgomery inverse require even modulo
+            return modInverseLorencz(a, p);
+        }
+        
+        int m = p.numberLength * 32;
+        
+        Object[] save = almostMonInv(a, p);
+        BigInteger r = ( (BigInteger) save[0] );
+        int k = ( (Integer) save[1] ).intValue();
+//		assert ( k >= n && k <= m + n );
+        
+        long n1 = calcN(p);
+        
+        if (k > m) {
+            r = monPro(r, BigInteger.ONE, p, n1);
+            k = k - m;
+//			assert k < m;
+        }
+        r = monPro(r, BigInteger.ONE.shiftLeft(m - k), p, n1);
+        return r;
+        }
+    
+    /**
+     * Calculate the first digit of the inverse
      */
-    static BigInteger modInverse(BigInteger a, BigInteger m) {
-        // PRE: (a > 0) and (m > 0)
-        BigInteger u, v, r, s, temp;
-        // u = MAX(a,m), v = MIN(a,m)
-        if (a.compareTo(m) == BigInteger.LESS) {
-            u = m;
-            v = a;
-            r = BigInteger.ZERO;
-            s = BigInteger.ONE;
+    private static long calcN(BigInteger a) {
+        long m0 = a.digits[0] & 0xFFFFFFFFL;
+        long n2 = 1L; // this is a'[0]
+        long powerOfTwo = 2L;
+        
+        do {
+            if (( ( m0 * n2 ) & powerOfTwo ) != 0) {
+                n2 |= powerOfTwo;
+        }
+            powerOfTwo <<= 1;
+        } while (powerOfTwo < 0x100000000L);
+        n2 = -n2;
+        return n2;
+        }
+    
+    /**
+     * Used for an intermediate result of the modInverse algorithm
+     * @return the pair: ((BigInteger)r, (Integer)k) where r == a^(-1) * 2^k mod (module)
+     */
+    private static Object[] almostMonInv(BigInteger a, BigInteger module) {
+        // PRE: a \in [1, p - 1]
+        BigInteger u, v, r, s;
+        // make copy to use inplace method
+        u = module.copy();
+        v = a.copy();
+        
+        int max = Math.max(v.numberLength, u.numberLength);
+        r = new BigInteger(1, 1, new int[max + 1]);
+        s = new BigInteger(1, 1, new int[max + 1]);
+        s.digits[0] = 1;
+        // s == 1 && v == 0
+        
+        int k = 0;
+        
+        int lsbu = u.getLowestSetBit();
+        int lsbv = v.getLowestSetBit();
+        int toShift;
+        
+        if (lsbu > lsbv) {
+            BitLevel.inplaceShiftRight(u, lsbu);
+            BitLevel.inplaceShiftRight(v, lsbv);
+            BitLevel.inplaceShiftLeft(r, lsbv);
+            k += lsbu - lsbv;
         } else {
-            v = m;
-            u = a;
-            s = BigInteger.ZERO;
-            r = BigInteger.ONE;
-        }
-        int uLen = u.bitLength();
-        int vLen = v.bitLength();
-        int f = uLen - vLen;
-
-        while (vLen > 1) {
-            if (u.sign == v.sign) {
-                u = u.subtract(v.shiftLeft(f));
-                r = r.subtract(s.shiftLeft(f));
+            BitLevel.inplaceShiftRight(u, lsbu);
+            BitLevel.inplaceShiftRight(v, lsbv);
+            BitLevel.inplaceShiftLeft(s, lsbu);
+            k += lsbv - lsbu;
+            
+    }
+
+        r.sign = 1;
+        while (v.signum() > 0) {
+            // INV v >= 0, u >= 0, v odd, u odd (excepto last iteration when v is even (0))
+    
+            while (u.compareTo(v) > BigInteger.EQUALS) {
+                Elementary.inplaceSubtract(u, v);
+                toShift = u.getLowestSetBit();
+                BitLevel.inplaceShiftRight(u, toShift);
+                Elementary.inplaceAdd(r, s);
+                BitLevel.inplaceShiftLeft(s, toShift);
+                k += toShift;
+                
+            }
+            
+            while (u.compareTo(v) <= BigInteger.EQUALS) {
+                Elementary.inplaceSubtract(v, u);
+                if (v.signum() == 0)
+                    break;
+                toShift = v.getLowestSetBit();
+                BitLevel.inplaceShiftRight(v, toShift);
+                Elementary.inplaceAdd(s, r);
+                BitLevel.inplaceShiftLeft(r, toShift);
+                k += toShift;
+                
+            }
+            
+        }
+        if (!u.isOne()){
+            // in u is stored the gcd
+            // math.19: BigInteger not invertible.
+            throw new ArithmeticException(Messages.getString("math.19"));
+        }
+        
+        if (r.compareTo(module) >= BigInteger.EQUALS) {
+            Elementary.inplaceSubtract(r, module);
+        }
+        r = module.subtract(r);
+        
+        return new Object[] { r, k };
+    }
+    
+    /**
+     * @return bi == abs(2^exp)
+     */
+    private static boolean isPowerOfTwo(BigInteger bi, int exp) {
+        boolean result = false;
+        result = ( exp >> 5 == bi.numberLength - 1 )
+        && ( bi.digits[bi.numberLength - 1] == 1 << ( exp & 31 ) );
+        if (result) {
+            for (int i = 0; result && i < bi.numberLength - 1; i++) {
+                result = bi.digits[i] == 0;
+            }
+        }
+        return result;
+    }
+    
+    /**
+     * Calculate how many iteration of Lorencz's algorithm would perform the
+     * same operation
+     *
+     * @param bi
+     * @param n
+     * @return
+     */
+    private static int howManyIterations(BigInteger bi, int n) {
+        int i = n - 1;
+        if (bi.sign > 0) {
+            while (!bi.testBit(i))
+                i--;
+            return n - 1 - i;
+        } else {
+            while (bi.testBit(i))
+                i--;
+            return n - 1 - Math.max(i, bi.getLowestSetBit());
+        }
+        
+    }
+    
+    /**
+     *
+     * Based on "New Algorithm for Classical Modular Inverse" Róbert Lórencz.
+     * LNCS 2523 (2002)
+     *
+     * @return a^(-1) mod m
+     */
+    static BigInteger modInverseLorencz(BigInteger a, BigInteger modulo) {
+        // PRE: a is coprime with modulo, a < modulo
+        
+        int max = Math.max(a.numberLength, modulo.numberLength);
+        int uDigits[] = new int[max + 1]; // enough place to make all the inplace operation
+        int vDigits[] = new int[max + 1];
+        System.arraycopy(modulo.digits, 0, uDigits, 0, modulo.numberLength);
+        System.arraycopy(a.digits, 0, vDigits, 0, a.numberLength);
+        BigInteger u = new BigInteger(modulo.sign, modulo.numberLength,
+                uDigits);
+        BigInteger v = new BigInteger(a.sign, a.numberLength, vDigits);
+        
+        BigInteger r = new BigInteger(0, 1, new int[max + 1]); // BigInteger.ZERO;
+        BigInteger s = new BigInteger(1, 1, new int[max + 1]);
+        s.digits[0] = 1;
+        // r == 0 && s == 1, but with enough place
+        
+        int coefU = 0, coefV = 0;
+        int n = modulo.bitLength();
+        int k;
+        while (!isPowerOfTwo(u, coefU) && !isPowerOfTwo(v, coefV)) {
+            
+            // modification of orignal algorithm: I calculate how many times the algorithm will enter in the same branch of if
+            k = howManyIterations(u, n);
+            
+            if (k != 0) {
+                BitLevel.inplaceShiftLeft(u, k);
+                if (coefU >= coefV) {
+                    BitLevel.inplaceShiftLeft(r, k);
+                } else {
+                    BitLevel.inplaceShiftRight(s, Math.min(coefV - coefU, k));
+                    if (k - ( coefV - coefU ) > 0) {
+                        BitLevel.inplaceShiftLeft(r, k - coefV + coefU);
+                    }
+                }
+                coefU += k;
+            }
+            
+            k = howManyIterations(v, n);
+            if (k != 0) {
+                BitLevel.inplaceShiftLeft(v, k);
+                if (coefV >= coefU) {
+                    BitLevel.inplaceShiftLeft(s, k);
+                } else {
+                    BitLevel.inplaceShiftRight(r, Math.min(coefU - coefV, k));
+                    if (k - ( coefU - coefV ) > 0) {
+                        BitLevel.inplaceShiftLeft(s, k - coefU + coefV);
+                    }
+                }
+                coefV += k;
+                
+            }
+            
+            if (u.signum() == v.signum()) {
+                if (coefU <= coefV) {
+                    Elementary.completeInPlaceSubtract(u, v);
+                    Elementary.completeInPlaceSubtract(r, s);
+                } else {
+                    Elementary.completeInPlaceSubtract(v, u);
+                    Elementary.completeInPlaceSubtract(s, r);
+                }
             } else {
-                u = u.add(v.shiftLeft(f));
-                r = r.add(s.shiftLeft(f));
+                if (coefU <= coefV) {
+                    Elementary.completeInPlaceAdd(u, v);
+                    Elementary.completeInPlaceAdd(r, s);
+                } else {
+                    Elementary.completeInPlaceAdd(v, u);
+                    Elementary.completeInPlaceAdd(s, r);
+                }
+            }
+            if (v.signum() == 0 || u.signum() == 0){
+                // math.19: BigInteger not invertible
+                throw new ArithmeticException(Messages.getString("math.19"));
             }
-            uLen = u.abs().bitLength();
-            vLen = v.abs().bitLength();
-            f = uLen - vLen;
-            if (f < 0) {
-                // SWAP(u,v)
-                temp = u;
-                u = v;
-                v = temp;
-                // SWAP(r,s)
-                temp = r;
-                r = s;
-                s = temp;
-
-                f = -f;
-                vLen = uLen;
-            }
-        }
-        if (v.sign == 0) {
-            return BigInteger.ZERO;
-        }
-        if (v.sign < 0) {
-            s = s.negate();
         }
-        if (s.compareTo(m) == BigInteger.GREATER) {
-            return s.subtract(m);
+        
+        if (isPowerOfTwo(v, coefV)) {
+            r = s;
+            if (v.signum() != u.signum())
+                u = u.negate();
+        }
+        if (u.testBit(n)) {
+            if (r.signum() < 0) {
+                r = r.negate();
+            } else {
+                r = modulo.subtract(r);
+            }
         }
-        if (s.sign < 0) {
-            return s.add(m);
+        if (r.signum() < 0) {
+            r = r.add(modulo);
         }
-        return s; // a^(-1) mod m
+        
+        return r;
     }
-
     
-    /*Implements the Montgomery modular exponentiation based in <i>The square and multiply algorithm and the
-     * Montgomery Reduction</i>.
-     *@ar.org.fitc.ref "C. K. Koc - Analyzing and Comparing Montgomery
-     *                  Multiplication Algorithms"
-     *@see #oddModPow(BigInteger, BigInteger,
-     *                           BigInteger)
-     */
     static BigInteger squareAndMultiply(BigInteger x2, BigInteger a2, BigInteger exponent,BigInteger modulus, long n2  ){
         BigInteger res = x2;
         for (int i = exponent.bitLength() - 1; i >= 0; i--) {

Modified: incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Elementary.java
URL: http://svn.apache.org/viewvc/incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Elementary.java?view=diff&rev=474164&r1=474163&r2=474164
==============================================================================
--- incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Elementary.java (original)
+++ incubator/harmony/enhanced/classlib/trunk/modules/math/src/main/java/java/math/Elementary.java Sun Nov 12 22:03:27 2006
@@ -33,7 +33,8 @@
 class Elementary {
 
     /** Just to denote that this class can't be instantied */
-    private Elementary() {}
+	private Elementary() {
+	}
 
     /**
      * Compares two arrays. All elements are treated as unsigned integers. The
@@ -114,18 +115,20 @@
     }
 
     /**
-     * Performs {@code res = a + b}. It is assumed the magnitude of a is not
-     * less than the magnitude of b.
+	 * Performs {@code res = a + b}. 
      */
     private static void add(int res[], int a[], int aSize, int b[], int bSize) {
-        // PRE: a[] >= b[]
+		// PRE: a.length < max(aSize, bSize)
+
         int i;
-        long carry = (a[0] & 0xFFFFFFFFL) + (b[0] & 0xFFFFFFFFL);
+		long carry = ( a[0] & 0xFFFFFFFFL ) + ( b[0] & 0xFFFFFFFFL );
 
         res[0] = (int) carry;
         carry >>= 32;
+
+		if (aSize >= bSize) {
         for (i = 1; i < bSize; i++) {
-            carry += (a[i] & 0xFFFFFFFFL) + (b[i] & 0xFFFFFFFFL);
+				carry += ( a[i] & 0xFFFFFFFFL ) + ( b[i] & 0xFFFFFFFFL );
             res[i] = (int) carry;
             carry >>= 32;
         }
@@ -134,7 +137,93 @@
             res[i] = (int) carry;
             carry >>= 32;
         }
+		} else {
+			for (i = 1; i < aSize; i++) {
+				carry += ( a[i] & 0xFFFFFFFFL ) + ( b[i] & 0xFFFFFFFFL );
         res[i] = (int) carry;
+				carry >>= 32;
+			}
+			for (; i < bSize; i++) {
+				carry += b[i] & 0xFFFFFFFFL;
+				res[i] = (int) carry;
+				carry >>= 32;
+			}
+		}
+		if (carry != 0) {
+			res[i] = (int) carry;
+		}
+	}
+
+	/** @see BigInteger#subtract(BigInteger) */
+	static BigInteger subtract(BigInteger op1, BigInteger op2) {
+		int resSign;
+		int resDigits[];
+		int op1Sign = op1.sign;
+		int op2Sign = op2.sign;
+
+		if (op2Sign == 0) {
+			return op1;
+		}
+		if (op1Sign == 0) {
+			return op2.negate ();
+		}
+		int op1Len = op1.numberLength;
+		int op2Len = op2.numberLength;
+		if (op1Len + op2Len == 2) {
+			long a = ( op1.digits[0] & 0xFFFFFFFFL );
+			long b = ( op2.digits[0] & 0xFFFFFFFFL );
+			if (op1Sign < 0) {
+				a = -a;
+			}
+			if (op2Sign < 0) {
+				b = -b;
+			}
+			return BigInteger.valueOf (a - b);
+		}
+		int cmp = ( ( op1Len != op2Len ) ? ( ( op1Len > op2Len ) ? 1 : -1 )
+				: Elementary.compareArrays (op1.digits, op2.digits, op1Len) );
+
+		if (cmp == BigInteger.LESS) {
+			resSign = -op2Sign;
+			resDigits = ( op1Sign == op2Sign ) ? subtract (op2.digits, op2Len,
+					op1.digits, op1Len) : add (op2.digits, op2Len, op1.digits,
+					op1Len);
+		} else {
+			resSign = op1Sign;
+			if (op1Sign == op2Sign) {
+				if (cmp == BigInteger.EQUALS) {
+					return BigInteger.ZERO;
+				}
+				resDigits = subtract (op1.digits, op1Len, op2.digits, op2Len);
+			} else {
+				resDigits = add (op1.digits, op1Len, op2.digits, op2Len);
+			}
+		}
+		BigInteger res = new BigInteger (resSign, resDigits.length, resDigits);
+		res.cutOffLeadingZeroes ();
+		return res;
+	}
+
+	/**
+	 * Performs {@code res = a - b}. It is assumed the magnitude of a is not
+	 * less than the magnitude of b.
+	 */
+	private static void subtract(int res[], int a[], int aSize, int b[],
+			int bSize) {
+		// PRE: a[] >= b[]
+		int i;
+		long borrow = 0;
+
+		for (i = 0; i < bSize; i++) {
+			borrow += ( a[i] & 0xFFFFFFFFL ) - ( b[i] & 0xFFFFFFFFL );
+			res[i] = (int) borrow;
+			borrow >>= 32; // -1 or 0
+		}
+		for (; i < aSize; i++) {
+			borrow += a[i] & 0xFFFFFFFFL;
+			res[i] = (int) borrow;
+			borrow >>= 32; // -1 or 0
+		}
     }
 
     /**
@@ -161,9 +250,12 @@
      */
     static void inplaceAdd(BigInteger op1, BigInteger op2) {
         // PRE: op1 >= op2 > 0
-        add(op1.digits, op1.digits, op1.numberLength, op2.digits,
+		add (op1.digits, op1.digits, op1.numberLength, op2.digits,
                 op2.numberLength);
-        op1.cutOffLeadingZeroes();
+		op1.numberLength = Math.min (Math.max (op1.numberLength,
+				op2.numberLength) + 1, op1.digits.length);
+		op1.cutOffLeadingZeroes ();
+		op1.unCache();
     }
 
     /**
@@ -192,80 +284,60 @@
             op1.digits[op1.numberLength] = 1;
             op1.numberLength++;
         }
+		op1.unCache();
     }
 
-    /** @see BigInteger#subtract(BigInteger) */
-    static BigInteger subtract(BigInteger op1, BigInteger op2) {
-        int resSign;
-        int resDigits[];
-        int op1Sign = op1.sign;
-        int op2Sign = op2.sign;
-
-        if (op2Sign == 0) {
-            return op1;
-        }
-        if (op1Sign == 0) {
-            return op2.negate();
-        }
-        int op1Len = op1.numberLength;
-        int op2Len = op2.numberLength;
-        if (op1Len + op2Len == 2) {
-            long a = (op1.digits[0] & 0xFFFFFFFFL);
-            long b = (op2.digits[0] & 0xFFFFFFFFL);
-            if (op1Sign < 0) {
-                a = -a;
-            }
-            if (op2Sign < 0) {
-                b = -b;
-            }
-            return BigInteger.valueOf(a - b);
-        }
-        int cmp = ((op1Len != op2Len) ? ((op1Len > op2Len) ? 1 : -1)
-                : Elementary.compareArrays(op1.digits, op2.digits, op1Len));
-
-        if (cmp == BigInteger.LESS) {
-            resSign = -op2Sign;
-            resDigits = (op1Sign == op2Sign) ? subtract(op2.digits, op2Len,
-                    op1.digits, op1Len) : add(op2.digits, op2Len, op1.digits,
-                    op1Len);
-        } else {
-            resSign = op1Sign;
-            if (op1Sign == op2Sign) {
-                if (cmp == BigInteger.EQUALS) {
-                    return BigInteger.ZERO;
-                }
-                resDigits = subtract(op1.digits, op1Len, op2.digits, op2Len);
-            } else {
-                resDigits = add(op1.digits, op1Len, op2.digits, op2Len);
-            }
-        }
-        BigInteger res = new BigInteger(resSign, resDigits.length, resDigits);
-        res.cutOffLeadingZeroes();
-        return res;
+	/**
+	 * Performs {@code op1 -= op2}. {@code op1} must have enough place to store
+	 * the result (i.e. {@code op1.bitLength() >= op2.bitLength()}). Both
+	 * should be positive (what implies that {@code op1 >= op2}).
+	 * 
+	 * @param op1
+	 *            the input minuend, and the ouput result.
+	 * @param op2
+	 *            the subtrahend
+	 */
+	static void inplaceSubtract(BigInteger op1, BigInteger op2) {
+		// PRE: op1 >= op2 > 0
+		subtract (op1.digits, op1.digits, op1.numberLength, op2.digits,
+				op2.numberLength);
+		op1.cutOffLeadingZeroes ();
+		op1.unCache();
     }
 
     /**
-     * Performs {@code res = a - b}. It is assumed the magnitude of a is not
-     * less than the magnitude of b.
+	 * Performs {@code res = b - a}
      */
-    private static void subtract(int res[], int a[], int aSize, int b[],
+	private static void inverseSubtract(int res[], int a[], int aSize, int b[],
             int bSize) {
-        // PRE: a[] >= b[]
         int i;
         long borrow = 0;
-
+		if (aSize < bSize) {
+			for (i = 0; i < aSize; i++) {
+				borrow += ( b[i] & 0xFFFFFFFFL ) - ( a[i] & 0xFFFFFFFFL );
+				res[i] = (int) borrow;
+				borrow >>= 32; // -1 or 0
+			}
+			for (; i < bSize; i++) {
+				borrow += b[i] & 0xFFFFFFFFL;
+				res[i] = (int) borrow;
+				borrow >>= 32; // -1 or 0
+			}
+		} else {
         for (i = 0; i < bSize; i++) {
-            borrow += (a[i] & 0xFFFFFFFFL) - (b[i] & 0xFFFFFFFFL);
+				borrow += ( b[i] & 0xFFFFFFFFL ) - ( a[i] & 0xFFFFFFFFL );
             res[i] = (int) borrow;
             borrow >>= 32; // -1 or 0
         }
         for (; i < aSize; i++) {
-            borrow += a[i] & 0xFFFFFFFFL;
+				borrow -= a[i] & 0xFFFFFFFFL;
             res[i] = (int) borrow;
             borrow >>= 32; // -1 or 0
         }
     }
 
+	}
+
     /**
      * Subtracts the value represented by {@code b} from the value represented
      * by {@code a}. It is assumed the magnitude of a is not less than the
@@ -281,18 +353,90 @@
     }
 
     /**
-     * Performs {@code op1 -= op2}. {@code op1} must have enough place to store
-     * the result (i.e. {@code op1.bitLength() >= op2.bitLength()}). Both
-     * should be positive (i.e. {@code op1 >= op2}).
+	 * Same as
      * 
-     * @param op1 the input minuend, and the ouput result.
-     * @param op2 the subtrahend
+	 * @link #inplaceSubtract(BigInteger, BigInteger), but without the
+	 *       restriction of non-positive values
+	 * @param op1
+	 *            should have enough space to save the result
+	 * @param op2
      */
-    static void inplaceSubtract(BigInteger op1, BigInteger op2) {
-        // PRE: op1 >= op2 > 0
-        subtract(op1.digits, op1.digits, op1.numberLength, op2.digits,
+	static void completeInPlaceSubtract(BigInteger op1, BigInteger op2) {
+		int resultSign = op1.compareTo (op2);
+		if (op1.sign == 0) {
+			System.arraycopy (op2.digits, 0, op1.digits, 0, op2.numberLength);
+			op1.sign = -op2.sign;
+		} else if (op1.sign != op2.sign) {
+			add (op1.digits, op1.digits, op1.numberLength, op2.digits,
                 op2.numberLength);
-        op1.cutOffLeadingZeroes();
-    }
+			op1.sign = resultSign;
+		} else {
+			int sign = unsignedArraysCompare (op1.digits,
+					op2.digits, op1.numberLength, op2.numberLength);
+			if (sign > 0) {
+				subtract (op1.digits, op1.digits, op1.numberLength, op2.digits,
+						op2.numberLength);	// op1 = op1 - op2
+				// op1.sign remains equal
+			} else {
+				inverseSubtract (op1.digits, op1.digits, op1.numberLength,
+						op2.digits, op2.numberLength);	// op1 = op2 - op1
+				op1.sign = -op1.sign;
+			}
+		}
+		op1.numberLength = Math.max (op1.numberLength, op2.numberLength) + 1;
+		op1.cutOffLeadingZeroes ();
+		op1.unCache();
+	}
+
+	/**
+	 * Same as @link #inplaceAdd(BigInteger, BigInteger), but without the restriction of
+	 *       non-positive values
+	 * @param op1 any number
+	 * @param op2 any number
+	 */
+	static void completeInPlaceAdd(BigInteger op1, BigInteger op2) {
+		if (op1.sign == 0)
+			System.arraycopy (op2.digits, 0, op1.digits, 0, op2.numberLength);
+		else if (op2.sign == 0)
+			return;
+		else if (op1.sign == op2.sign)
+			add (op1.digits, op1.digits, op1.numberLength, op2.digits,
+					op2.numberLength);
+		else {
+			int sign = unsignedArraysCompare(op1.digits,
+					op2.digits, op1.numberLength, op2.numberLength);
+			if (sign > 0)
+				subtract (op1.digits, op1.digits, op1.numberLength, op2.digits,
+						op2.numberLength);
+			else {
+				inverseSubtract (op1.digits, op1.digits, op1.numberLength,
+						op2.digits, op2.numberLength);
+				op1.sign = -op1.sign;
+			}
+    }
+		op1.numberLength = Math.max (op1.numberLength, op2.numberLength) + 1;
+		op1.cutOffLeadingZeroes ();
+		op1.unCache();
+	}
+
+	/**
+	 * Compares two arrays, representing unsigned integer in little-endian order.
+	 * Returns +1,0,-1 if a is - respective - greater, equal or lesser then b 
+	 */
+	private static int unsignedArraysCompare(int[] a, int[] b, int aSize, int bSize){
+		if (aSize > bSize)
+			return 1;
+		else if (aSize < bSize)
+			return -1;
+		
+		else {
+			int i;
+			for (i = aSize - 1; i >= 0 && a[i] == b[i]; i-- )
+				;
+			return i < 0 ? BigInteger.EQUALS : (( a[i] & 0xFFFFFFFFL ) < (b[i] & 0xFFFFFFFFL ) ? BigInteger.LESS
+						: BigInteger.GREATER) ;
+		}
+	}
+	
 
 }