in commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java [1097:1153]
private static Complex divide(double re1, double im1, double re2, double im2) {
double a = re1;
double b = im1;
double c = re2;
double d = im2;
int ilogbw = 0;
// Get the exponent to scale the divisor parts to the range [1, 2).
final int exponent = getScale(c, d);
if (exponent <= Double.MAX_EXPONENT) {
ilogbw = exponent;
c = Math.scalb(c, -ilogbw);
d = Math.scalb(d, -ilogbw);
}
final double denom = c * c + d * d;
// Note: Modification from the listing in ISO C99 G.5.1 (8):
// Avoid overflow if a or b are very big.
// Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
// when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
// (bc - ad) with large negative values.
// Use the maximum exponent as an approximation to the magnitude.
if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
ilogbw -= 2;
a /= 4;
b /= 4;
}
double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
// Recover infinities and zeros that computed as NaN+iNaN
// the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
if (Double.isNaN(x) && Double.isNaN(y)) {
if (denom == 0.0 &&
(!Double.isNaN(a) || !Double.isNaN(b))) {
// nonzero/zero
// This case produces the same result as divide by a real-only zero
// using Complex.divide(+/-0.0)
x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
} else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
Double.isFinite(c) && Double.isFinite(d)) {
// infinite/finite
a = boxInfinity(a);
b = boxInfinity(b);
x = Double.POSITIVE_INFINITY * (a * c + b * d);
y = Double.POSITIVE_INFINITY * (b * c - a * d);
} else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
Double.isFinite(a) && Double.isFinite(b)) {
// finite/infinite
c = boxInfinity(c);
d = boxInfinity(d);
x = 0.0 * (a * c + b * d);
y = 0.0 * (b * c - a * d);
}
}
return new Complex(x, y);
}