in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/DD.java [945:1010]
static DD ldexp(double x, double xx, int exp, DD r) {
// Handle scaling when 2^n can be represented with a single normal number
// n >= -1022 && n <= 1023
// Using unsigned compare => n + 1022 <= 1023 + 1022
if (exp + CMP_UNSIGNED_1022 < CMP_UNSIGNED_2046) {
final double s = twoPow(exp);
r.hi = x * s;
r.lo = xx * s;
return r;
}
// Scale by multiples of 2^512 (largest representable power of 2).
// Scaling requires max 5 multiplications to under/overflow any normal value.
// Break this down into e.g.: 2^512^(exp / 512) * 2^(exp % 512)
// Number of multiples n = exp / 512 : exp >>> 9
// Remainder m = exp % 512 : exp & 511 (exp must be positive)
int n;
int m;
double p;
if (exp < 0) {
// Downscaling
// (Note: Using an unsigned shift handles negation of min value: -2^31)
n = -exp >>> 9;
// m = exp % 512
m = -(-exp & 511);
p = TWO_POW_M512;
} else {
// Upscaling
n = exp >>> 9;
m = exp & 511;
p = TWO_POW_512;
}
// Multiply by the remainder scaling factor first. The remaining multiplications
// are either 2^512 or 2^-512.
// Down-scaling to sub-normal will use the final multiplication into a sub-normal result.
// Note here that n >= 1 as the n in [-1022, 1023] case has been handled.
// Handle n : 1, 2, 3, 4, 5
if (n >= 5) {
// n >= 5 will be over/underflow. Use an extreme scale factor.
// Do not use +/- infinity as this creates NaN if x = 0.
// p -> 2^1023 or 2^-1025
p *= p * 0.5;
r.hi = x * p * p * p;
r.lo = xx * p * p * p;
return r;
}
final double s = twoPow(m);
if (n == 4) {
r.hi = x * s * p * p * p * p;
r.lo = xx * s * p * p * p * p;
} else if (n == 3) {
r.hi = x * s * p * p * p;
r.lo = xx * s * p * p * p;
} else if (n == 2) {
r.hi = x * s * p * p;
r.lo = xx * s * p * p;
} else {
// n = 1. Occurs only if exp = -1023.
r.hi = x * s * p;
r.lo = xx * s * p;
}
return r;
}