--- old/src/share/classes/java/lang/Math.java 2011-09-16 18:10:37.000000000 -0700 +++ new/src/share/classes/java/lang/Math.java 2011-09-16 18:10:37.000000000 -0700 @@ -26,6 +26,8 @@ package java.lang; import java.util.Random; +import sun.misc.FloatConsts; +import sun.misc.DoubleConsts; /** * The class {@code Math} contains methods for performing basic @@ -963,7 +965,31 @@ * @since 1.5 */ public static double ulp(double d) { - return sun.misc.FpUtils.ulp(d); + int exp = getExponent(d); + + switch(exp) { + case DoubleConsts.MAX_EXPONENT+1: // NaN or infinity + return Math.abs(d); + + case DoubleConsts.MIN_EXPONENT-1: // zero or subnormal + return Double.MIN_VALUE; + + default: + assert exp <= DoubleConsts.MAX_EXPONENT && exp >= DoubleConsts.MIN_EXPONENT; + + // ulp(x) is usually 2^(SIGNIFICAND_WIDTH-1)*(2^ilogb(x)) + exp = exp - (DoubleConsts.SIGNIFICAND_WIDTH-1); + if (exp >= DoubleConsts.MIN_EXPONENT) { + return powerOfTwoD(exp); + } + else { + // return a subnormal result; left shift integer + // representation of Double.MIN_VALUE appropriate + // number of positions + return Double.longBitsToDouble(1L << + (exp - (DoubleConsts.MIN_EXPONENT - (DoubleConsts.SIGNIFICAND_WIDTH-1)) )); + } + } } /** @@ -990,7 +1016,31 @@ * @since 1.5 */ public static float ulp(float f) { - return sun.misc.FpUtils.ulp(f); + int exp = getExponent(f); + + switch(exp) { + case FloatConsts.MAX_EXPONENT+1: // NaN or infinity + return Math.abs(f); + + case FloatConsts.MIN_EXPONENT-1: // zero or subnormal + return FloatConsts.MIN_VALUE; + + default: + assert exp <= FloatConsts.MAX_EXPONENT && exp >= FloatConsts.MIN_EXPONENT; + + // ulp(x) is usually 2^(SIGNIFICAND_WIDTH-1)*(2^ilogb(x)) + exp = exp - (FloatConsts.SIGNIFICAND_WIDTH-1); + if (exp >= FloatConsts.MIN_EXPONENT) { + return powerOfTwoF(exp); + } + else { + // return a subnormal result; left shift integer + // representation of FloatConsts.MIN_VALUE appropriate + // number of positions + return Float.intBitsToFloat(1 << + (exp - (FloatConsts.MIN_EXPONENT - (FloatConsts.SIGNIFICAND_WIDTH-1)) )); + } + } } /** @@ -1011,7 +1061,7 @@ * @since 1.5 */ public static double signum(double d) { - return sun.misc.FpUtils.signum(d); + return (d == 0.0 || Double.isNaN(d))?d:copySign(1.0, d); } /** @@ -1032,7 +1082,7 @@ * @since 1.5 */ public static float signum(float f) { - return sun.misc.FpUtils.signum(f); + return (f == 0.0f || Float.isNaN(f))?f:copySign(1.0f, f); } /** @@ -1252,7 +1302,11 @@ * @since 1.6 */ public static double copySign(double magnitude, double sign) { - return sun.misc.FpUtils.rawCopySign(magnitude, sign); + return Double.longBitsToDouble((Double.doubleToRawLongBits(sign) & + (DoubleConsts.SIGN_BIT_MASK)) | + (Double.doubleToRawLongBits(magnitude) & + (DoubleConsts.EXP_BIT_MASK | + DoubleConsts.SIGNIF_BIT_MASK))); } /** @@ -1271,7 +1325,11 @@ * @since 1.6 */ public static float copySign(float magnitude, float sign) { - return sun.misc.FpUtils.rawCopySign(magnitude, sign); + return Float.intBitsToFloat((Float.floatToRawIntBits(sign) & + (FloatConsts.SIGN_BIT_MASK)) | + (Float.floatToRawIntBits(magnitude) & + (FloatConsts.EXP_BIT_MASK | + FloatConsts.SIGNIF_BIT_MASK))); } /** @@ -1289,7 +1347,13 @@ * @since 1.6 */ public static int getExponent(float f) { - return sun.misc.FpUtils.getExponent(f); + /* + * Bitwise convert f to integer, mask out exponent bits, shift + * to the right and then subtract out float's bias adjust to + * get true exponent value + */ + return ((Float.floatToRawIntBits(f) & FloatConsts.EXP_BIT_MASK) >> + (FloatConsts.SIGNIFICAND_WIDTH - 1)) - FloatConsts.EXP_BIAS; } /** @@ -1307,7 +1371,13 @@ * @since 1.6 */ public static int getExponent(double d) { - return sun.misc.FpUtils.getExponent(d); + /* + * Bitwise convert d to long, mask out exponent bits, shift + * to the right and then subtract out double's bias adjust to + * get true exponent value. + */ + return (int)(((Double.doubleToRawLongBits(d) & DoubleConsts.EXP_BIT_MASK) >> + (DoubleConsts.SIGNIFICAND_WIDTH - 1)) - DoubleConsts.EXP_BIAS); } /** @@ -1351,7 +1421,63 @@ * @since 1.6 */ public static double nextAfter(double start, double direction) { - return sun.misc.FpUtils.nextAfter(start, direction); + /* + * The cases: + * + * nextAfter(+infinity, 0) == MAX_VALUE + * nextAfter(+infinity, +infinity) == +infinity + * nextAfter(-infinity, 0) == -MAX_VALUE + * nextAfter(-infinity, -infinity) == -infinity + * + * are naturally handled without any additional testing + */ + + // First check for NaN values + if (Double.isNaN(start) || Double.isNaN(direction)) { + // return a NaN derived from the input NaN(s) + return start + direction; + } else if (start == direction) { + return direction; + } else { // start > direction or start < direction + // Add +0.0 to get rid of a -0.0 (+0.0 + -0.0 => +0.0) + // then bitwise convert start to integer. + long transducer = Double.doubleToRawLongBits(start + 0.0d); + + /* + * IEEE 754 floating-point numbers are lexicographically + * ordered if treated as signed- magnitude integers . + * Since Java's integers are two's complement, + * incrementing" the two's complement representation of a + * logically negative floating-point value *decrements* + * the signed-magnitude representation. Therefore, when + * the integer representation of a floating-point values + * is less than zero, the adjustment to the representation + * is in the opposite direction than would be expected at + * first . + */ + if (direction > start) { // Calculate next greater value + transducer = transducer + (transducer >= 0L ? 1L:-1L); + } else { // Calculate next lesser value + assert direction < start; + if (transducer > 0L) + --transducer; + else + if (transducer < 0L ) + ++transducer; + /* + * transducer==0, the result is -MIN_VALUE + * + * The transition from zero (implicitly + * positive) to the smallest negative + * signed magnitude value must be done + * explicitly. + */ + else + transducer = DoubleConsts.SIGN_BIT_MASK | 1L; + } + + return Double.longBitsToDouble(transducer); + } } /** @@ -1394,7 +1520,63 @@ * @since 1.6 */ public static float nextAfter(float start, double direction) { - return sun.misc.FpUtils.nextAfter(start, direction); + /* + * The cases: + * + * nextAfter(+infinity, 0) == MAX_VALUE + * nextAfter(+infinity, +infinity) == +infinity + * nextAfter(-infinity, 0) == -MAX_VALUE + * nextAfter(-infinity, -infinity) == -infinity + * + * are naturally handled without any additional testing + */ + + // First check for NaN values + if (Float.isNaN(start) || Double.isNaN(direction)) { + // return a NaN derived from the input NaN(s) + return start + (float)direction; + } else if (start == direction) { + return (float)direction; + } else { // start > direction or start < direction + // Add +0.0 to get rid of a -0.0 (+0.0 + -0.0 => +0.0) + // then bitwise convert start to integer. + int transducer = Float.floatToRawIntBits(start + 0.0f); + + /* + * IEEE 754 floating-point numbers are lexicographically + * ordered if treated as signed- magnitude integers . + * Since Java's integers are two's complement, + * incrementing" the two's complement representation of a + * logically negative floating-point value *decrements* + * the signed-magnitude representation. Therefore, when + * the integer representation of a floating-point values + * is less than zero, the adjustment to the representation + * is in the opposite direction than would be expected at + * first. + */ + if (direction > start) {// Calculate next greater value + transducer = transducer + (transducer >= 0 ? 1:-1); + } else { // Calculate next lesser value + assert direction < start; + if (transducer > 0) + --transducer; + else + if (transducer < 0 ) + ++transducer; + /* + * transducer==0, the result is -MIN_VALUE + * + * The transition from zero (implicitly + * positive) to the smallest negative + * signed magnitude value must be done + * explicitly. + */ + else + transducer = FloatConsts.SIGN_BIT_MASK | 1; + } + + return Float.intBitsToFloat(transducer); + } } /** @@ -1423,7 +1605,13 @@ * @since 1.6 */ public static double nextUp(double d) { - return sun.misc.FpUtils.nextUp(d); + if( Double.isNaN(d) || d == Double.POSITIVE_INFINITY) + return d; + else { + d += 0.0d; + return Double.longBitsToDouble(Double.doubleToRawLongBits(d) + + ((d >= 0.0d)?+1L:-1L)); + } } /** @@ -1452,7 +1640,13 @@ * @since 1.6 */ public static float nextUp(float f) { - return sun.misc.FpUtils.nextUp(f); + if( Float.isNaN(f) || f == FloatConsts.POSITIVE_INFINITY) + return f; + else { + f += 0.0f; + return Float.intBitsToFloat(Float.floatToRawIntBits(f) + + ((f >= 0.0f)?+1:-1)); + } } @@ -1487,7 +1681,80 @@ * @since 1.6 */ public static double scalb(double d, int scaleFactor) { - return sun.misc.FpUtils.scalb(d, scaleFactor); + /* + * This method does not need to be declared strictfp to + * compute the same correct result on all platforms. When + * scaling up, it does not matter what order the + * multiply-store operations are done; the result will be + * finite or overflow regardless of the operation ordering. + * However, to get the correct result when scaling down, a + * particular ordering must be used. + * + * When scaling down, the multiply-store operations are + * sequenced so that it is not possible for two consecutive + * multiply-stores to return subnormal results. If one + * multiply-store result is subnormal, the next multiply will + * round it away to zero. This is done by first multiplying + * by 2 ^ (scaleFactor % n) and then multiplying several + * times by by 2^n as needed where n is the exponent of number + * that is a covenient power of two. In this way, at most one + * real rounding error occurs. If the double value set is + * being used exclusively, the rounding will occur on a + * multiply. If the double-extended-exponent value set is + * being used, the products will (perhaps) be exact but the + * stores to d are guaranteed to round to the double value + * set. + * + * It is _not_ a valid implementation to first multiply d by + * 2^MIN_EXPONENT and then by 2 ^ (scaleFactor % + * MIN_EXPONENT) since even in a strictfp program double + * rounding on underflow could occur; e.g. if the scaleFactor + * argument was (MIN_EXPONENT - n) and the exponent of d was a + * little less than -(MIN_EXPONENT - n), meaning the final + * result would be subnormal. + * + * Since exact reproducibility of this method can be achieved + * without any undue performance burden, there is no + * compelling reason to allow double rounding on underflow in + * scalb. + */ + + // magnitude of a power of two so large that scaling a finite + // nonzero value by it would be guaranteed to over or + // underflow; due to rounding, scaling down takes takes an + // additional power of two which is reflected here + final int MAX_SCALE = DoubleConsts.MAX_EXPONENT + -DoubleConsts.MIN_EXPONENT + + DoubleConsts.SIGNIFICAND_WIDTH + 1; + int exp_adjust = 0; + int scale_increment = 0; + double exp_delta = Double.NaN; + + // Make sure scaling factor is in a reasonable range + + if(scaleFactor < 0) { + scaleFactor = Math.max(scaleFactor, -MAX_SCALE); + scale_increment = -512; + exp_delta = twoToTheDoubleScaleDown; + } + else { + scaleFactor = Math.min(scaleFactor, MAX_SCALE); + scale_increment = 512; + exp_delta = twoToTheDoubleScaleUp; + } + + // Calculate (scaleFactor % +/-512), 512 = 2^9, using + // technique from "Hacker's Delight" section 10-2. + int t = (scaleFactor >> 9-1) >>> 32 - 9; + exp_adjust = ((scaleFactor + t) & (512 -1)) - t; + + d *= powerOfTwoD(exp_adjust); + scaleFactor -= exp_adjust; + + while(scaleFactor != 0) { + d *= exp_delta; + scaleFactor -= scale_increment; + } + return d; } /** @@ -1521,6 +1788,49 @@ * @since 1.6 */ public static float scalb(float f, int scaleFactor) { - return sun.misc.FpUtils.scalb(f, scaleFactor); + // magnitude of a power of two so large that scaling a finite + // nonzero value by it would be guaranteed to over or + // underflow; due to rounding, scaling down takes takes an + // additional power of two which is reflected here + final int MAX_SCALE = FloatConsts.MAX_EXPONENT + -FloatConsts.MIN_EXPONENT + + FloatConsts.SIGNIFICAND_WIDTH + 1; + + // Make sure scaling factor is in a reasonable range + scaleFactor = Math.max(Math.min(scaleFactor, MAX_SCALE), -MAX_SCALE); + + /* + * Since + MAX_SCALE for float fits well within the double + * exponent range and + float -> double conversion is exact + * the multiplication below will be exact. Therefore, the + * rounding that occurs when the double product is cast to + * float will be the correctly rounded float result. Since + * all operations other than the final multiply will be exact, + * it is not necessary to declare this method strictfp. + */ + return (float)((double)f*powerOfTwoD(scaleFactor)); + } + + // Constants used in scalb + static double twoToTheDoubleScaleUp = powerOfTwoD(512); + static double twoToTheDoubleScaleDown = powerOfTwoD(-512); + + /** + * Returns a floating-point power of two in the normal range. + */ + static double powerOfTwoD(int n) { + assert(n >= DoubleConsts.MIN_EXPONENT && n <= DoubleConsts.MAX_EXPONENT); + return Double.longBitsToDouble((((long)n + (long)DoubleConsts.EXP_BIAS) << + (DoubleConsts.SIGNIFICAND_WIDTH-1)) + & DoubleConsts.EXP_BIT_MASK); + } + + /** + * Returns a floating-point power of two in the normal range. + */ + public static float powerOfTwoF(int n) { + assert(n >= FloatConsts.MIN_EXPONENT && n <= FloatConsts.MAX_EXPONENT); + return Float.intBitsToFloat(((n + FloatConsts.EXP_BIAS) << + (FloatConsts.SIGNIFICAND_WIDTH-1)) + & FloatConsts.EXP_BIT_MASK); } }