diff options
Diffstat (limited to 'libjava/classpath/gnu/java/math/MPN.java')
| -rw-r--r-- | libjava/classpath/gnu/java/math/MPN.java | 596 |
1 files changed, 298 insertions, 298 deletions
diff --git a/libjava/classpath/gnu/java/math/MPN.java b/libjava/classpath/gnu/java/math/MPN.java index 05491bb..9143ef0 100644 --- a/libjava/classpath/gnu/java/math/MPN.java +++ b/libjava/classpath/gnu/java/math/MPN.java @@ -7,7 +7,7 @@ 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 @@ -59,9 +59,9 @@ public class MPN long carry = (long) y & 0xffffffffL; for (int i = 0; i < size; i++) { - carry += ((long) x[i] & 0xffffffffL); - dest[i] = (int) carry; - carry >>= 32; + carry += ((long) x[i] & 0xffffffffL); + dest[i] = (int) carry; + carry >>= 32; } return (int) carry; } @@ -77,10 +77,10 @@ public class MPN 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; + carry += ((long) x[i] & 0xffffffffL) + + ((long) y[i] & 0xffffffffL); + dest[i] = (int) carry; + carry >>>= 32; } return (int) carry; } @@ -96,15 +96,15 @@ public class MPN 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; + 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; } @@ -142,23 +142,23 @@ public class MPN * This function is basically the same gmp's mpn_mul. */ public static void mul (int[] dest, - int[] x, int xlen, - int[] y, int ylen) + 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; + 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; } } @@ -174,71 +174,71 @@ public class MPN 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; - } + 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)) { + 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 { + 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; - } - } + 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); @@ -251,7 +251,7 @@ public class MPN */ public static int divmod_1 (int[] quotient, int[] dividend, - int len, int divisor) + int len, int divisor) { int i = len - 1; long r = dividend[i]; @@ -259,16 +259,16 @@ public class MPN r = 0; else { - quotient[i--] = 0; - r <<= 32; + 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; + int n0 = dividend[i]; + r = (r & ~0xffffffffL) | (n0 & 0xffffffffL); + r = udiv_qrnnd (r, divisor); + quotient[i] = (int) r; } return (int)(r >> 32); } @@ -285,19 +285,19 @@ public class MPN 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; + 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; @@ -328,37 +328,37 @@ public class MPN 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); + // 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; + { + 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); } @@ -372,19 +372,19 @@ public class MPN { 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; + 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; @@ -410,9 +410,9 @@ public class MPN for (int k = 16; k > 0; k = k >> 1) { int j = i >>> k; if (j == 0) - count += k; + count += k; else - i = j; + i = j; } return count; } @@ -422,61 +422,61 @@ public class MPN 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; + // The base is a power of 2. Read the input string from + // least to most significant character/digit. */ + + int next_bitpos = 0; + int bits_per_indigit = 0; + for (int i = base; (i >>= 1) != 0; ) bits_per_indigit++; + int res_digit = 0; + + for (int i = str_len; --i >= 0; ) + { + int inp_digit = str[i]; + res_digit |= inp_digit << next_bitpos; + next_bitpos += bits_per_indigit; + if (next_bitpos >= 32) + { + dest[size++] = res_digit; + next_bitpos -= 32; + res_digit = inp_digit >> (bits_per_indigit - next_bitpos); + } + } + + if (res_digit != 0) + dest[size++] = res_digit; } else { - // General case. The base is not a power of 2. - int indigits_per_limb = MPN.chars_per_word (base); - int str_pos = 0; - - while (str_pos < str_len) - { - int chunk = str_len - str_pos; - if (chunk > indigits_per_limb) - chunk = indigits_per_limb; - int res_digit = str[str_pos++]; - int big_base = base; - - while (--chunk > 0) - { - res_digit = res_digit * base + str[str_pos++]; - big_base *= base; - } - - int cy_limb; - if (size == 0) - cy_limb = res_digit; - else - { - cy_limb = MPN.mul_1 (dest, dest, size, big_base); - cy_limb += MPN.add_1 (dest, dest, size, res_digit); - } - if (cy_limb != 0) - dest[size++] = cy_limb; - } + // General case. The base is not a power of 2. + int indigits_per_limb = MPN.chars_per_word (base); + int str_pos = 0; + + while (str_pos < str_len) + { + int chunk = str_len - str_pos; + if (chunk > indigits_per_limb) + chunk = indigits_per_limb; + int res_digit = str[str_pos++]; + int big_base = base; + + while (--chunk > 0) + { + res_digit = res_digit * base + str[str_pos++]; + big_base *= base; + } + + int cy_limb; + if (size == 0) + cy_limb = res_digit; + else + { + cy_limb = MPN.mul_1 (dest, dest, size, big_base); + cy_limb += MPN.add_1 (dest, dest, size, res_digit); + } + if (cy_limb != 0) + dest[size++] = cy_limb; + } } return size; } @@ -489,22 +489,22 @@ public class MPN { 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; - } + 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. - * + * * @return -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) @@ -521,7 +521,7 @@ public class MPN * Assumes: 0 < count < 32 */ public static int rshift (int[] dest, int[] x, int x_start, - int len, int count) + int len, int count) { int count_2 = 32 - count; int low_word = x[x_start]; @@ -529,9 +529,9 @@ public class MPN 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; + 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; @@ -546,13 +546,13 @@ public class MPN * 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) + 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]; + dest[i] = x[i + x_start]; } /** Return the long-truncated value of right shifting. @@ -571,10 +571,10 @@ public class MPN 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)); + 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); } @@ -587,7 +587,7 @@ public class MPN */ public static int lshift (int[] dest, int d_offset, - int[] x, int len, int count) + int[] x, int len, int count) { int count_2 = 32 - count; int i = len - 1; @@ -596,9 +596,9 @@ public class MPN 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; + 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; @@ -611,13 +611,13 @@ public class MPN int i = 0; while ((word & 0xF) == 0) { - word >>= 4; - i += 4; + word >>= 4; + i += 4; } if ((word & 3) == 0) { - word >>= 2; - i += 2; + word >>= 2; + i += 2; } if ((word & 1) == 0) i += 1; @@ -630,8 +630,8 @@ public class MPN { for (int i = 0; ; i++) { - if (words[i] != 0) - return 32 * i + findLowestBit (words[i]); + if (words[i] != 0) + return 32 * i + findLowestBit (words[i]); } } @@ -646,12 +646,12 @@ public class MPN // 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; - } + 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); @@ -666,69 +666,69 @@ public class MPN int[] other_arg; /* The other one can be even or odd. */ if ((x[0] & 1) != 0) { - odd_arg = x; - other_arg = y; + odd_arg = x; + other_arg = y; } else { - odd_arg = y; - other_arg = x; + odd_arg = y; + other_arg = x; } for (;;) { - // Shift other_arg until it is odd; this doesn't - // affect the gcd, since we divide by 2**k, which does not - // divide odd_arg. - for (i = 0; other_arg[i] == 0; ) i++; - if (i > 0) - { - int j; - for (j = 0; j < len-i; j++) - other_arg[j] = other_arg[j+i]; - for ( ; j < len; j++) - other_arg[j] = 0; - } - i = findLowestBit(other_arg[0]); - if (i > 0) - MPN.rshift (other_arg, other_arg, 0, len, i); - - // Now both odd_arg and other_arg are odd. - - // Subtract the smaller from the larger. - // This does not change the result, since gcd(a-b,b)==gcd(a,b). - i = MPN.cmp(odd_arg, other_arg, len); - if (i == 0) - break; - if (i > 0) - { // odd_arg > other_arg - MPN.sub_n (odd_arg, odd_arg, other_arg, len); - // Now odd_arg is even, so swap with other_arg; - int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp; - } - else - { // other_arg > odd_arg - MPN.sub_n (other_arg, other_arg, odd_arg, len); - } - while (odd_arg[len-1] == 0 && other_arg[len-1] == 0) - len--; + // Shift other_arg until it is odd; this doesn't + // affect the gcd, since we divide by 2**k, which does not + // divide odd_arg. + for (i = 0; other_arg[i] == 0; ) i++; + if (i > 0) + { + int j; + for (j = 0; j < len-i; j++) + other_arg[j] = other_arg[j+i]; + for ( ; j < len; j++) + other_arg[j] = 0; + } + i = findLowestBit(other_arg[0]); + if (i > 0) + MPN.rshift (other_arg, other_arg, 0, len, i); + + // Now both odd_arg and other_arg are odd. + + // Subtract the smaller from the larger. + // This does not change the result, since gcd(a-b,b)==gcd(a,b). + i = MPN.cmp(odd_arg, other_arg, len); + if (i == 0) + break; + if (i > 0) + { // odd_arg > other_arg + MPN.sub_n (odd_arg, odd_arg, other_arg, len); + // Now odd_arg is even, so swap with other_arg; + int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp; + } + else + { // other_arg > odd_arg + MPN.sub_n (other_arg, other_arg, odd_arg, len); + } + while (odd_arg[len-1] == 0 && other_arg[len-1] == 0) + len--; } if (initShiftWords + initShiftBits > 0) { - if (initShiftBits > 0) - { - int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits); - if (sh_out != 0) - x[(len++)+initShiftWords] = sh_out; - } - else - { - for (i = len; --i >= 0;) - x[i+initShiftWords] = x[i]; - } - for (i = initShiftWords; --i >= 0; ) - x[i] = 0; - len += initShiftWords; + if (initShiftBits > 0) + { + int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits); + if (sh_out != 0) + x[(len++)+initShiftWords] = sh_out; + } + else + { + for (i = len; --i >= 0;) + x[i+initShiftWords] = x[i]; + } + for (i = initShiftWords; --i >= 0; ) + x[i] = 0; + len += initShiftWords; } return len; } @@ -761,9 +761,9 @@ public class MPN ps.print('('); for (int i = 0; i < len; i++) { - if (i > 0) - ps.print (' '); - ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16)); + if (i > 0) + ps.print (' '); + ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16)); } ps.print(')'); } |