in commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java [3528:3640]
private static double hypot(double x, double y) {
// Differences to the fdlibm reference:
//
// 1. fdlibm orders the two parts using the magnitude of the upper 32-bits.
// This incorrectly orders numbers which differ only in the lower 32-bits.
// This invalidates hypot(x, y) = hypot(y, x) for small sub-normal numbers and a minority
// of cases of normal numbers. This implementation forces the |x| >= |y| order
// using the entire 63-bits of the unsigned doubles to ensure the function
// is commutative.
//
// 2. fdlibm computed scaling by directly writing changes to the exponent bits
// and maintained the high part (ha) during scaling for use in the high
// precision sum x^2 + y^2. Since exponent scaling cannot be applied to sub-normals
// the original version dropped the split number representation for sub-normals
// and can produce maximum errors above 1 ULP for sub-normal numbers.
// This version uses Dekker's method to split the number. This can be applied to
// sub-normals and allows dropping the condition to check for sub-normal numbers
// since all small numbers are handled with a single scaling factor.
// The effect is increased precision for the majority of sub-normal cases where
// the implementations compute a different result.
//
// 3. An alteration is done here to add an 'else if' instead of a second
// 'if' statement. Thus you cannot scale down and up at the same time.
//
// 4. There is no use of the absolute double value. The magnitude comparison is
// performed using the long bit representation. The computation x^2+y^2 is
// insensitive to the sign bit. Thus use of Math.abs(double) is only in edge-case
// branches.
//
// 5. The exponent different to ignore the smaller component has changed from 60 to 54.
//
// Original comments from fdlibm are in c style: /* */
// Extra comments added for reference.
//
// Note that the high 32-bits are compared to constants.
// The lowest 20-bits are the upper bits of the 52-bit mantissa.
// The next 11-bits are the biased exponent. The sign bit has been cleared.
// Scaling factors are powers of two for exact scaling.
// For clarity the values have been refactored to named constants.
// The mask is used to remove the sign bit.
final long xbits = Double.doubleToRawLongBits(x) & UNSIGN_MASK;
final long ybits = Double.doubleToRawLongBits(y) & UNSIGN_MASK;
// Order by magnitude: |a| >= |b|
double a;
double b;
/* High word of x & y */
int ha;
int hb;
if (ybits > xbits) {
a = y;
b = x;
ha = (int) (ybits >>> 32);
hb = (int) (xbits >>> 32);
} else {
a = x;
b = y;
ha = (int) (xbits >>> 32);
hb = (int) (ybits >>> 32);
}
// Check if the smaller part is significant.
// a^2 is computed in extended precision for an effective mantissa of 106-bits.
// An exponent difference of 54 is where b^2 will not overlap a^2.
if ((ha - hb) > EXP_54) {
/* a/b > 2**54 */
// or a is Inf or NaN.
// No addition of a + b for sNaN.
return Math.abs(a);
}
double rescale = 1.0;
if (ha > EXP_500) {
/* a > 2^500 */
if (ha >= EXP_1024) {
/* Inf or NaN */
// Check b is infinite for the IEEE754 result.
// No addition of a + b for sNaN.
return Math.abs(b) == Double.POSITIVE_INFINITY ?
Double.POSITIVE_INFINITY :
Math.abs(a);
}
/* scale a and b by 2^-600 */
// Before scaling: a in [2^500, 2^1023].
// After scaling: a in [2^-100, 2^423].
// After scaling: b in [2^-154, 2^423].
a *= TWO_POW_NEG_600;
b *= TWO_POW_NEG_600;
rescale = TWO_POW_600;
} else if (hb < EXP_NEG_500) {
// No special handling of sub-normals.
// These do not matter when we do not manipulate the exponent bits
// for scaling the split representation.
// Intentional comparison with zero.
if (b == 0) {
return Math.abs(a);
}
/* scale a and b by 2^600 */
// Effective min exponent of a sub-normal = -1022 - 52 = -1074.
// Before scaling: b in [2^-1074, 2^-501].
// After scaling: b in [2^-474, 2^99].
// After scaling: a in [2^-474, 2^153].
a *= TWO_POW_600;
b *= TWO_POW_600;
rescale = TWO_POW_NEG_600;
}
// High precision x^2 + y^2
return Math.sqrt(x2y2(a, b)) * rescale;
}