From 1a1de0dfd7ba4b4fa2be1e464280755e9df917ab Mon Sep 17 00:00:00 2001 From: Ulf Adams Date: Mon, 28 May 2018 21:20:13 +0200 Subject: [PATCH 1/2] Rewrite Java implementation to use 32-bit factors This puts it in line with the C implementation fallback code. This may be required in order to half the lookup table size, which may require higher bit counts, and the previous implementation was limited to 4*31 = 124 bits. This may also be a deadend. If the code becomes significantly slower in a way that we cannot recover from, then we may have to undo this. --- src/main/java/info/adams/ryu/RyuDouble.java | 106 ++++++++++---------- 1 file changed, 55 insertions(+), 51 deletions(-) diff --git a/src/main/java/info/adams/ryu/RyuDouble.java b/src/main/java/info/adams/ryu/RyuDouble.java index 210f3999..b96cf8bd 100644 --- a/src/main/java/info/adams/ryu/RyuDouble.java +++ b/src/main/java/info/adams/ryu/RyuDouble.java @@ -46,16 +46,17 @@ public final class RyuDouble { private static final BigInteger[] POW5_INV = new BigInteger[NEG_TABLE_SIZE]; private static final int POW5_BITCOUNT = 121; // max 3*31 = 124 - private static final int POW5_QUARTER_BITCOUNT = 31; + private static final int POW5_QUARTER_BITCOUNT = 32; private static final int[][] POW5_SPLIT = new int[POS_TABLE_SIZE][4]; private static final int POW5_INV_BITCOUNT = 122; // max 3*31 = 124 - private static final int POW5_INV_QUARTER_BITCOUNT = 31; + private static final int POW5_INV_QUARTER_BITCOUNT = 32; private static final int[][] POW5_INV_SPLIT = new int[NEG_TABLE_SIZE][4]; static { BigInteger mask = BigInteger.valueOf(1).shiftLeft(POW5_QUARTER_BITCOUNT).subtract(BigInteger.ONE); BigInteger invMask = BigInteger.valueOf(1).shiftLeft(POW5_INV_QUARTER_BITCOUNT).subtract(BigInteger.ONE); + for (int i = 0; i < Math.max(POW5.length, POW5_INV.length); i++) { BigInteger pow = BigInteger.valueOf(5).pow(i); int pow5len = pow.bitLength(); @@ -68,10 +69,8 @@ public final class RyuDouble { } if (i < POW5_SPLIT.length) { for (int j = 0; j < 4; j++) { - POW5_SPLIT[i][j] = pow - .shiftRight(pow5len - POW5_BITCOUNT + (3 - j) * POW5_QUARTER_BITCOUNT) - .and(mask) - .intValueExact(); + POW5_SPLIT[i][j] = toUnsignedIntExact( + pow.shiftRight(pow5len - POW5_BITCOUNT + (3 - j) * POW5_QUARTER_BITCOUNT).and(mask)); } } @@ -81,16 +80,21 @@ public final class RyuDouble { BigInteger inv = BigInteger.ONE.shiftLeft(j).divide(pow).add(BigInteger.ONE); POW5_INV[i] = inv; for (int k = 0; k < 4; k++) { - if (k == 0) { - POW5_INV_SPLIT[i][k] = inv.shiftRight((3 - k) * POW5_INV_QUARTER_BITCOUNT).intValueExact(); - } else { - POW5_INV_SPLIT[i][k] = inv.shiftRight((3 - k) * POW5_INV_QUARTER_BITCOUNT).and(invMask).intValueExact(); - } + POW5_INV_SPLIT[i][k] = toUnsignedIntExact( + inv.shiftRight((3 - k) * POW5_INV_QUARTER_BITCOUNT).and(invMask)); } } } } + private static int toUnsignedIntExact(BigInteger v) { + long value = v.longValueExact(); + if (value > 0xffffffffL) { + throw new ArithmeticException("value too large for int"); + } + return (int) value; + } + public static void main(String[] args) { DEBUG = true; double value = Double.longBitsToDouble(0x7fefffffffffffffL); @@ -181,9 +185,9 @@ public static String doubleToString(double value, RoundingMode roundingMode) { e10 = q; if (DEBUG) { System.out.println(mv + " * 2^" + e2); - System.out.println("V+=" + dp); - System.out.println("V =" + dv); - System.out.println("V-=" + dm); + System.out.println("V+=" + Long.toUnsignedString(dp)); + System.out.println("V =" + Long.toUnsignedString(dv)); + System.out.println("V-=" + Long.toUnsignedString(dm)); } if (DEBUG) { long exact = POW5_INV[q] @@ -435,26 +439,7 @@ private static int pow5Factor(long value) { * point number. */ private static long mulPow5divPow2(long m, int i, int j) { - // m has at most 55 bits. - long mHigh = m >>> 31; - long mLow = m & 0x7fffffff; - long bits13 = mHigh * POW5_SPLIT[i][0]; // 124 - long bits03 = mLow * POW5_SPLIT[i][0]; // 93 - long bits12 = mHigh * POW5_SPLIT[i][1]; // 93 - long bits02 = mLow * POW5_SPLIT[i][1]; // 62 - long bits11 = mHigh * POW5_SPLIT[i][2]; // 62 - long bits01 = mLow * POW5_SPLIT[i][2]; // 31 - long bits10 = mHigh * POW5_SPLIT[i][3]; // 31 - long bits00 = mLow * POW5_SPLIT[i][3]; // 0 - int actualShift = j - 3 * 31 - 21; - if (actualShift < 0) { - throw new IllegalArgumentException("" + actualShift); - } - return (((((( - ((bits00 >>> 31) + bits01 + bits10) >>> 31) - + bits02 + bits11) >>> 31) - + bits03 + bits12) >>> 21) - + (bits13 << 10)) >>> actualShift; + return mulAndShiftRight(m, POW5_SPLIT[i], j); } /** @@ -462,26 +447,45 @@ private static long mulPow5divPow2(long m, int i, int j) { * decimal digits. i and j are already chosen appropriately. */ private static long mulPow5InvDivPow2(long m, int i, int j) { + return mulAndShiftRight(m, POW5_INV_SPLIT[i], j); + } + + private static long mulAndShiftRight(long m, int[] mul, int shift) { // m has at most 55 bits. - long mHigh = m >>> 31; - long mLow = m & 0x7fffffff; - long bits13 = mHigh * POW5_INV_SPLIT[i][0]; - long bits03 = mLow * POW5_INV_SPLIT[i][0]; - long bits12 = mHigh * POW5_INV_SPLIT[i][1]; - long bits02 = mLow * POW5_INV_SPLIT[i][1]; - long bits11 = mHigh * POW5_INV_SPLIT[i][2]; - long bits01 = mLow * POW5_INV_SPLIT[i][2]; - long bits10 = mHigh * POW5_INV_SPLIT[i][3]; - long bits00 = mLow * POW5_INV_SPLIT[i][3]; - - int actualShift = j - 3 * 31 - 21; + long mHi = m >>> 32; // 23 bit + long mLo = m & 0xffffffffL; + long bits11 = mHi * Integer.toUnsignedLong(mul[2]); // 64 + long bits01 = mLo * Integer.toUnsignedLong(mul[2]); // 32 + long bits10 = mHi * Integer.toUnsignedLong(mul[3]); // 32 + long bits00 = mLo * Integer.toUnsignedLong(mul[3]); // 0 + + long mLoMid = bits01 + bits10; + long carry96 = Long.compareUnsigned(mLoMid, bits01) < 0 ? 1 : 0; + long mLoLo = bits00 + (mLoMid << 32); + long carry64 = Long.compareUnsigned(mLoLo, bits00) < 0 ? 1 : 0; + long mLoHi = bits11 + (mLoMid >>> 32) + (carry96 << 32) + carry64; // Can't overflow, 64 + + long bits13 = mHi * Integer.toUnsignedLong(mul[0]); // 128 + long bits03 = mLo * Integer.toUnsignedLong(mul[0]); // 96 + long bits12 = mHi * Integer.toUnsignedLong(mul[1]); // 96 + long bits02 = mLo * Integer.toUnsignedLong(mul[1]); // 64 + + long mHiMid = bits03 + bits12; + long carry160 = Long.compareUnsigned(mHiMid, bits03) < 0 ? 1 : 0; + long mHiLo = bits02 + (mHiMid << 32); // 64 + long carry128 = Long.compareUnsigned(mHiLo, bits02) < 0 ? 1 : 0; + long mHiHi = bits13 + (mHiMid >>> 32) + (carry160 << 32) + carry128; // Can't overflow, 128 + + long lo = mHiLo + mLoHi; // 64 + long carry = Long.compareUnsigned(lo, mHiLo) < 0 ? 1 : 0; + long hi = mHiHi + carry; // 128 + + int actualShift = shift - 64; if (actualShift < 0) { throw new IllegalArgumentException("" + actualShift); } - return (((((( - ((bits00 >>> 31) + bits01 + bits10) >>> 31) - + bits02 + bits11) >>> 31) - + bits03 + bits12) >>> 21) - + (bits13 << 10)) >>> actualShift; + return actualShift >= 64 + ? hi >>> (actualShift - 64) + : (hi << (64 - actualShift)) + (lo >>> actualShift); } } From 20b070797c35db5a0837e566e675c0bde6a20163 Mon Sep 17 00:00:00 2001 From: Ulf Adams Date: Mon, 28 May 2018 21:49:26 +0200 Subject: [PATCH 2/2] Only access even-numbered entries in the lookup table Naively, this requires 65 bit. However, we can immediately divide by 10 since we know that the last digit is unnecessary, and that is the same as dividing by 2 and then by 5. We do this by increasing the shift distance by one, so it's essentially for free. However, this is somewhat costly in terms of performance. --- src/main/java/info/adams/ryu/RyuDouble.java | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/main/java/info/adams/ryu/RyuDouble.java b/src/main/java/info/adams/ryu/RyuDouble.java index b96cf8bd..b267d5ef 100644 --- a/src/main/java/info/adams/ryu/RyuDouble.java +++ b/src/main/java/info/adams/ryu/RyuDouble.java @@ -45,11 +45,11 @@ public final class RyuDouble { private static final BigInteger[] POW5 = new BigInteger[POS_TABLE_SIZE]; private static final BigInteger[] POW5_INV = new BigInteger[NEG_TABLE_SIZE]; - private static final int POW5_BITCOUNT = 121; // max 3*31 = 124 + private static final int POW5_BITCOUNT = 124; // max 128 private static final int POW5_QUARTER_BITCOUNT = 32; private static final int[][] POW5_SPLIT = new int[POS_TABLE_SIZE][4]; - private static final int POW5_INV_BITCOUNT = 122; // max 3*31 = 124 + private static final int POW5_INV_BITCOUNT = 124; // max 128 private static final int POW5_INV_QUARTER_BITCOUNT = 32; private static final int[][] POW5_INV_SPLIT = new int[NEG_TABLE_SIZE][4]; @@ -97,7 +97,7 @@ private static int toUnsignedIntExact(BigInteger v) { public static void main(String[] args) { DEBUG = true; - double value = Double.longBitsToDouble(0x7fefffffffffffffL); + double value = 1.8531501765868567E21; //Double.longBitsToDouble(0x7fefffffffffffffL); String result = doubleToString(value); System.out.println(result + " " + value); } @@ -176,13 +176,20 @@ public static String doubleToString(double value, RoundingMode roundingMode) { boolean dpIsTrailingZeros = false, dmIsTrailingZeros = false; if (e2 >= 0) { int q = Math.max(0, (int) (e2 * LOG10_2_NUMERATOR / LOG10_2_DENOMINATOR) - 1); + boolean cutOne = (q & 1) != 0; + e10 = q; + q &= ~1; // k = constant + floor(log_2(5^q)) int k = POW5_INV_BITCOUNT + pow5bits(q) - 1; - int i = -e2 + q + k; + int i = -e2 + q + k + (cutOne ? 1 : 0); dv = mulPow5InvDivPow2(mv, q, i); dp = mulPow5InvDivPow2(mp, q, i); dm = mulPow5InvDivPow2(mm, q, i); - e10 = q; + if (cutOne) { + dv = Long.divideUnsigned(dv, 5); + dp = Long.divideUnsigned(dp, 5); + dm = Long.divideUnsigned(dm, 5); + } if (DEBUG) { System.out.println(mv + " * 2^" + e2); System.out.println("V+=" + Long.toUnsignedString(dp)); @@ -192,7 +199,9 @@ public static String doubleToString(double value, RoundingMode roundingMode) { if (DEBUG) { long exact = POW5_INV[q] .multiply(BigInteger.valueOf(mv)) - .shiftRight(-e2 + q + k).longValueExact(); + .shiftRight(-e2 + q + k) + .divide(cutOne ? BigInteger.TEN : BigInteger.ONE) + .longValueExact(); System.out.println(exact + " " + POW5_INV[q].bitCount()); if (dv != exact) { throw new IllegalStateException();