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} &times;

@@ -1485,11 +1679,84 @@
      * @param scaleFactor power of 2 used to scale {@code d}
      * @return {@code d} &times; 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} &times;
      * 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} &times; 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);
     }
 }