FpUtils contains static utility methods for
* manipulating and inspecting float and
* double floating-point numbers. These methods include
* functionality recommended or required by the IEEE 754
* floating-point standard.
*
* @author Joseph D. Darcy
*/
public class FpUtils {
/*
* The methods in this class are reasonably implemented using
* direct or indirect bit-level manipulation of floating-point
* values. However, having access to the IEEE 754 recommended
* functions would obviate the need for most programmers to engage
* in floating-point bit-twiddling.
*
* An IEEE 754 number has three fields, from most significant bit
* to to least significant, sign, exponent, and significand.
*
* msb lsb
* [sign|exponent| fractional_significand]
*
* Using some encoding cleverness, explained below, the high order
* bit of the logical significand does not need to be explicitly
* stored, thus "fractional_significand" instead of simply
* "significand" in the figure above.
*
* For finite normal numbers, the numerical value encoded is
*
* (-1)^sign * 2^(exponent)*(1.fractional_significand)
*
* Most finite floating-point numbers are normalized; the exponent
* value is reduced until the leading significand bit is 1.
* Therefore, the leading 1 is redundant and is not explicitly
* stored. If a numerical value is so small it cannot be
* normalized, it has a subnormal representation. Subnormal
* numbers don't have a leading 1 in their significand; subnormals
* are encoding using a special exponent value. In other words,
* the high-order bit of the logical significand can be elided in
* from the representation in either case since the bit's value is
* implicit from the exponent value.
*
* The exponent field uses a biased representation; if the bits of
* the exponent are interpreted as a unsigned integer E, the
* exponent represented is E - E_bias where E_bias depends on the
* floating-point format. E can range between E_min and E_max,
* constants which depend on the floating-point format. E_min and
* E_max are -126 and +127 for float, -1022 and +1023 for double.
*
* The 32-bit float format has 1 sign bit, 8 exponent bits, and 23
* bits for the significand (which is logically 24 bits wide
* because of the implicit bit). The 64-bit double format has 1
* sign bit, 11 exponent bits, and 52 bits for the significand
* (logically 53 bits).
*
* Subnormal numbers and zero have the special exponent value
* E_min -1; the numerical value represented by a subnormal is:
*
* (-1)^sign * 2^(E_min)*(0.fractional_significand)
*
* Zero is represented by all zero bits in the exponent and all
* zero bits in the significand; zero can have either sign.
*
* Infinity and NaN are encoded using the exponent value E_max +
* 1. Signed infinities have all significand bits zero; NaNs have
* at least one non-zero significand bit.
*
* The details of IEEE 754 floating-point encoding will be used in
* the methods below without further comment. For further
* exposition on IEEE 754 numbers, see "IEEE Standard for Binary
* Floating-Point Arithmetic" ANSI/IEEE Std 754-1985 or William
* Kahan's "Lecture Notes on the Status of IEEE Standard 754 for
* Binary Floating-Point Arithmetic",
* http://www.cs.berkeley.edu/~wkahan/ieee754status/ieee754.ps.
*
* Many of this class's methods are members of the set of IEEE 754
* recommended functions or similar functions recommended or
* required by IEEE 754R. Discussion of various implementation
* techniques for these functions have occurred in:
*
* W.J. Cody and Jerome T. Coonen, "Algorithm 772 Functions to
* Support the IEEE Standard for Binary Floating-Point
* Arithmetic," ACM Transactions on Mathematical Software,
* vol. 19, no. 4, December 1993, pp. 443-451.
*
* Joseph D. Darcy, "Writing robust IEEE recommended functions in
* ``100% Pure Java''(TM)," University of California, Berkeley
* technical report UCB//CSD-98-1009.
*/
/**
* Don't let anyone instantiate this class.
*/
private FpUtils() {}
// Constants used in scalb
static double twoToTheDoubleScaleUp = powerOfTwoD(512);
static double twoToTheDoubleScaleDown = powerOfTwoD(-512);
// Helper Methods
// The following helper methods are used in the implementation of
// the public recommended functions; they generally omit certain
// tests for exception cases.
/**
* Returns unbiased exponent of a double.
*/
public static int getExponent(double 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 unbiased exponent of a float.
*/
public static int getExponent(float 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 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.
*/
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);
}
/**
* Returns the first floating-point argument with the sign of the
* second floating-point argument. Note that unlike the {@link
* FpUtils#copySign(double, double) copySign} method, this method
* does not require NaN sign arguments to be treated
* as positive values; implementations are permitted to treat some
* NaN arguments as positive and other NaN arguments as negative
* to allow greater performance.
*
* @param magnitude the parameter providing the magnitude of the result
* @param sign the parameter providing the sign of the result
* @return a value with the magnitude of magnitude
* and the sign of sign.
* @author Joseph D. Darcy
*/
public static double rawCopySign(double magnitude, double 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
* FpUtils#copySign(float, float) copySign} method, this method
* does not require NaN sign arguments to be treated
* as positive values; implementations are permitted to treat some
* NaN arguments as positive and other NaN arguments as negative
* to allow greater performance.
*
* @param magnitude the parameter providing the magnitude of the result
* @param sign the parameter providing the sign of the result
* @return a value with the magnitude of magnitude
* and the sign of sign.
* @author Joseph D. Darcy
*/
public static float rawCopySign(float magnitude, float sign) {
return Float.intBitsToFloat((Float.floatToRawIntBits(sign) &
(FloatConsts.SIGN_BIT_MASK)) |
(Float.floatToRawIntBits(magnitude) &
(FloatConsts.EXP_BIT_MASK |
FloatConsts.SIGNIF_BIT_MASK)));
}
/* ***************************************************************** */
/**
* Returns true if the argument is a finite
* floating-point value; returns false otherwise (for
* NaN and infinity arguments).
*
* @param d the double value to be tested
* @return true if the argument is a finite
* floating-point value, false otherwise.
*/
public static boolean isFinite(double d) {
return Math.abs(d) <= DoubleConsts.MAX_VALUE;
}
/**
* Returns true if the argument is a finite
* floating-point value; returns false otherwise (for
* NaN and infinity arguments).
*
* @param f the float value to be tested
* @return true if the argument is a finite
* floating-point value, false otherwise.
*/
public static boolean isFinite(float f) {
return Math.abs(f) <= FloatConsts.MAX_VALUE;
}
/**
* Returns true if the specified number is infinitely
* large in magnitude, false otherwise.
*
* Note that this method is equivalent to the {@link
* Double#isInfinite(double) Double.isInfinite} method; the
* functionality is included in this class for convenience.
*
* @param d the value to be tested.
* @return true if the value of the argument is positive
* infinity or negative infinity; false otherwise.
*/
public static boolean isInfinite(double d) {
return Double.isInfinite(d);
}
/**
* Returns true if the specified number is infinitely
* large in magnitude, false otherwise.
*
*
Note that this method is equivalent to the {@link
* Float#isInfinite(float) Float.isInfinite} method; the
* functionality is included in this class for convenience.
*
* @param f the value to be tested.
* @return true if the argument is positive infinity or
* negative infinity; false otherwise.
*/
public static boolean isInfinite(float f) {
return Float.isInfinite(f);
}
/**
* Returns true if the specified number is a
* Not-a-Number (NaN) value, false otherwise.
*
*
Note that this method is equivalent to the {@link
* Double#isNaN(double) Double.isNaN} method; the functionality is
* included in this class for convenience.
*
* @param d the value to be tested.
* @return true if the value of the argument is NaN;
* false otherwise.
*/
public static boolean isNaN(double d) {
return Double.isNaN(d);
}
/**
* Returns true if the specified number is a
* Not-a-Number (NaN) value, false otherwise.
*
*
Note that this method is equivalent to the {@link
* Float#isNaN(float) Float.isNaN} method; the functionality is
* included in this class for convenience.
*
* @param f the value to be tested.
* @return true if the argument is NaN;
* false otherwise.
*/
public static boolean isNaN(float f) {
return Float.isNaN(f);
}
/**
* Returns true if the unordered relation holds
* between the two arguments. When two floating-point values are
* unordered, one value is neither less than, equal to, nor
* greater than the other. For the unordered relation to be true,
* at least one argument must be a NaN.
*
* @param arg1 the first argument
* @param arg2 the second argument
* @return true if at least one argument is a NaN,
* false otherwise.
*/
public static boolean isUnordered(double arg1, double arg2) {
return isNaN(arg1) || isNaN(arg2);
}
/**
* Returns true if the unordered relation holds
* between the two arguments. When two floating-point values are
* unordered, one value is neither less than, equal to, nor
* greater than the other. For the unordered relation to be true,
* at least one argument must be a NaN.
*
* @param arg1 the first argument
* @param arg2 the second argument
* @return true if at least one argument is a NaN,
* false otherwise.
*/
public static boolean isUnordered(float arg1, float arg2) {
return isNaN(arg1) || isNaN(arg2);
}
/**
* Returns unbiased exponent of a double; for
* subnormal values, the number is treated as if it were
* normalized. That is for all finite, non-zero, positive numbers
* x, scalb(x, -ilogb(x)) is
* always in the range [1, 2).
*
* Special cases: *
float; for
* subnormal values, the number is treated as if it were
* normalized. That is for all finite, non-zero, positive numbers
* x, scalb(x, -ilogb(x)) is
* always in the range [1, 2).
* * Special cases: *
d ×
* 2scale_factor rounded as if performed
* by a single correctly rounded floating-point multiply to a
* member of the double value set. See §4.2.3
* of the Java
* Language Specification for a discussion of floating-point
* value sets. If the exponent of the result is between the
* double's minimum exponent and maximum exponent,
* the answer is calculated exactly. If the exponent of the
* result would be larger than doubles's maximum
* exponent, an infinity is returned. Note that if the result is
* subnormal, precision may be lost; that is, when scalb(x,
* n) is subnormal, scalb(scalb(x, n), -n) may
* not equal x. When the result is non-NaN, the result has
* the same sign as d.
*
** Special cases: *
d
* @return d * 2scale_factor
* @author Joseph D. Darcy
*/
public static double scalb(double d, int scale_factor) {
/*
* 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 ^ (scale_factor % 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 ^ (scale_factor %
* MIN_EXPONENT) since even in a strictfp program double
* rounding on underflow could occur; e.g. if the scale_factor
* 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(scale_factor < 0) {
scale_factor = Math.max(scale_factor, -MAX_SCALE);
scale_increment = -512;
exp_delta = twoToTheDoubleScaleDown;
}
else {
scale_factor = Math.min(scale_factor, MAX_SCALE);
scale_increment = 512;
exp_delta = twoToTheDoubleScaleUp;
}
// Calculate (scale_factor % +/-512), 512 = 2^9, using
// technique from "Hacker's Delight" section 10-2.
int t = (scale_factor >> 9-1) >>> 32 - 9;
exp_adjust = ((scale_factor + t) & (512 -1)) - t;
d *= powerOfTwoD(exp_adjust);
scale_factor -= exp_adjust;
while(scale_factor != 0) {
d *= exp_delta;
scale_factor -= scale_increment;
}
return d;
}
/**
* Return f ×
* 2scale_factor rounded as if performed
* by a single correctly rounded floating-point multiply to a
* member of the float value set. See §4.2.3
* of the Java
* Language Specification for a discussion of floating-point
* value set. If the exponent of the result is between the
* float's minimum exponent and maximum exponent, the
* answer is calculated exactly. If the exponent of the result
* would be larger than float's maximum exponent, an
* infinity is returned. Note that if the result is subnormal,
* precision may be lost; that is, when scalb(x, n)
* is subnormal, scalb(scalb(x, n), -n) may not equal
* x. When the result is non-NaN, the result has the same
* sign as f.
*
** Special cases: *
f
* @return f * 2scale_factor
* @author Joseph D. Darcy
*/
public static float scalb(float f, int scale_factor) {
// 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
scale_factor = Math.max(Math.min(scale_factor, 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(scale_factor));
}
/**
* Returns the floating-point number adjacent to the first
* argument in the direction of the second argument. If both
* arguments compare as equal the second argument is returned.
*
* * Special cases: *
direction
* is returned unchanged (as implied by the requirement of
* returning the second argument if the arguments compare as
* equal).
*
* start is
* ±Double.MIN_VALUE and direction
* has a value such that the result should have a smaller
* magnitude, then a zero with the same sign as start
* is returned.
*
* start is infinite and
* direction has a value such that the result should
* have a smaller magnitude, Double.MAX_VALUE with the
* same sign as start is returned.
*
* start is equal to ±
* Double.MAX_VALUE and direction has a
* value such that the result should have a larger magnitude, an
* infinity with same sign as start is returned.
* start's neighbors or start should
* be returned
* @return The floating-point number adjacent to start in the
* direction of direction.
* @author Joseph D. Darcy
*/
public static double nextAfter(double start, double 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 (isNaN(start) || 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
* arguments compare as equal, the second argument is returned.
*
* * Special cases: *
float
* zero with the same sign as direction is returned
* (as implied by the requirement of returning the second argument
* if the arguments compare as equal).
*
* start is
* ±Float.MIN_VALUE and direction
* has a value such that the result should have a smaller
* magnitude, then a zero with the same sign as start
* is returned.
*
* start is infinite and
* direction has a value such that the result should
* have a smaller magnitude, Float.MAX_VALUE with the
* same sign as start is returned.
*
* start is equal to ±
* Float.MAX_VALUE and direction has a
* value such that the result should have a larger magnitude, an
* infinity with same sign as start is returned.
* start's neighbors or start should
* be returned
* @return The floating-point number adjacent to start in the
* direction of direction.
* @author Joseph D. Darcy
*/
public static float nextAfter(float start, double 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 (isNaN(start) || 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 d in
* the direction of positive infinity. This method is
* semantically equivalent to nextAfter(d,
* Double.POSITIVE_INFINITY); however, a nextUp
* implementation may run faster than its equivalent
* nextAfter call.
*
* Special Cases: *
Double.MIN_VALUE
*
* f in
* the direction of positive infinity. This method is
* semantically equivalent to nextAfter(f,
* Double.POSITIVE_INFINITY); however, a nextUp
* implementation may run faster than its equivalent
* nextAfter call.
*
* Special Cases: *
Float.MIN_VALUE
*
* d in
* the direction of negative infinity. This method is
* semantically equivalent to nextAfter(d,
* Double.NEGATIVE_INFINITY); however, a
* nextDown implementation may run faster than its
* equivalent nextAfter call.
*
* Special Cases: *
-Double.MIN_VALUE
*
* f in
* the direction of negative infinity. This method is
* semantically equivalent to nextAfter(f,
* Float.NEGATIVE_INFINITY); however, a
* nextDown implementation may run faster than its
* equivalent nextAfter call.
*
* Special Cases: *
-Float.MIN_VALUE
*
* sign argument is always treated as if it were
* positive.
*
* @param magnitude the parameter providing the magnitude of the result
* @param sign the parameter providing the sign of the result
* @return a value with the magnitude of magnitude
* and the sign of sign.
* @author Joseph D. Darcy
* @since 1.5
*/
public static double copySign(double magnitude, double sign) {
return rawCopySign(magnitude, (isNaN(sign)?1.0d:sign));
}
/**
* Returns the first floating-point argument with the sign of the
* second floating-point argument. For this method, a NaN
* sign argument is always treated as if it were
* positive.
*
* @param magnitude the parameter providing the magnitude of the result
* @param sign the parameter providing the sign of the result
* @return a value with the magnitude of magnitude
* and the sign of sign.
* @author Joseph D. Darcy
*/
public static float copySign(float magnitude, float sign) {
return rawCopySign(magnitude, (isNaN(sign)?1.0f:sign));
}
/**
* Returns the size of an ulp of the argument. An ulp of a
* double value is the positive distance between this
* floating-point value and the double value next
* larger in magnitude. Note that for non-NaN x,
* ulp(-x) == ulp(x).
*
* Special Cases: *
Double.MIN_VALUE.
* Double.MAX_VALUE, then
* the result is equal to 2971.
* float value is the positive distance between this
* floating-point value and the float value next
* larger in magnitude. Note that for non-NaN x,
* ulp(-x) == ulp(x).
*
* Special Cases: *
Float.MIN_VALUE.
* Float.MAX_VALUE, then
* the result is equal to 2104.
* Special Cases: *
Special Cases: *