src/share/classes/java/lang/Math.java
Print this page
@@ -24,10 +24,12 @@
*/
package java.lang;
import java.util.Random;
+import sun.misc.FloatConsts;
+import sun.misc.DoubleConsts;
/**
* The class {@code Math} contains methods for performing basic
* numeric operations such as the elementary exponential, logarithm,
* square root, and trigonometric functions.
@@ -961,11 +963,35 @@
* @return the size of an ulp of the argument
* @author Joseph D. Darcy
* @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)) ));
+ }
+ }
}
/**
* Returns the size of an ulp of the argument. An ulp, unit in
* the last place, of a {@code float} value is the positive
@@ -988,11 +1014,35 @@
* @return the size of an ulp of the argument
* @author Joseph D. Darcy
* @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)) ));
+ }
+ }
}
/**
* Returns the signum function of the argument; zero if the argument
* is zero, 1.0 if the argument is greater than zero, -1.0 if the
@@ -1009,11 +1059,11 @@
* @return the signum function of the argument
* @author Joseph D. Darcy
* @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);
}
/**
* Returns the signum function of the argument; zero if the argument
* is zero, 1.0f if the argument is greater than zero, -1.0f if the
@@ -1030,11 +1080,11 @@
* @return the signum function of the argument
* @author Joseph D. Darcy
* @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);
}
/**
* Returns the hyperbolic sine of a {@code double} value.
* The hyperbolic sine of <i>x</i> is defined to be
@@ -1250,11 +1300,15 @@
* @return a value with the magnitude of {@code magnitude}
* and the sign of {@code sign}.
* @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)));
}
/**
* Returns the first floating-point argument with the sign of the
* second floating-point argument. Note that unlike the {@link
@@ -1269,11 +1323,15 @@
* @return a value with the magnitude of {@code magnitude}
* and the sign of {@code sign}.
* @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)));
}
/**
* Returns the unbiased exponent used in the representation of a
* {@code float}. Special cases:
@@ -1287,11 +1345,17 @@
* @param f a {@code float} value
* @return the unbiased exponent of the argument
* @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;
}
/**
* Returns the unbiased exponent used in the representation of a
* {@code double}. Special cases:
@@ -1305,11 +1369,17 @@
* @param d a {@code double} value
* @return the unbiased exponent of the argument
* @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);
}
/**
* Returns the floating-point number adjacent to the first
* argument in the direction of the second argument. If both
@@ -1349,11 +1419,67 @@
* @return The floating-point number adjacent to {@code start} in the
* direction of {@code direction}.
* @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);
+ }
}
/**
* Returns the floating-point number adjacent to the first
* argument in the direction of the second argument. If both
@@ -1392,11 +1518,67 @@
* @return The floating-point number adjacent to {@code start} in the
* direction of {@code direction}.
* @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);
+ }
}
/**
* Returns the floating-point value adjacent to {@code d} in
* the direction of positive infinity. This method is
@@ -1421,11 +1603,17 @@
* @return The adjacent floating-point value closer to positive
* infinity.
* @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));
+ }
}
/**
* Returns the floating-point value adjacent to {@code f} in
* the direction of positive infinity. This method is
@@ -1450,11 +1638,17 @@
* @return The adjacent floating-point value closer to positive
* infinity.
* @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));
+ }
}
/**
* Return {@code d} ×
@@ -1485,11 +1679,84 @@
* @param scaleFactor power of 2 used to scale {@code d}
* @return {@code d} × 2<sup>{@code scaleFactor}</sup>
* @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;
}
/**
* Return {@code f} ×
* 2<sup>{@code scaleFactor}</sup> rounded as if performed
@@ -1519,8 +1786,51 @@
* @param scaleFactor power of 2 used to scale {@code f}
* @return {@code f} × 2<sup>{@code scaleFactor}</sup>
* @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);
}
}