in commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java [1393:1472]
private Complex log(DoubleUnaryOperator log, double logOfeOver2, double logOf2, ComplexConstructor constructor) {
// Handle NaN
if (Double.isNaN(real) || Double.isNaN(imaginary)) {
// Return NaN unless infinite
if (isInfinite()) {
return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
}
// No-use of the input constructor
return NAN;
}
// Returns the real part:
// log(sqrt(x^2 + y^2))
// log(x^2 + y^2) / 2
// Compute with positive values
double x = Math.abs(real);
double y = Math.abs(imaginary);
// Find the larger magnitude.
if (x < y) {
final double tmp = x;
x = y;
y = tmp;
}
if (x == 0) {
// Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception.
return constructor.create(Double.NEGATIVE_INFINITY,
negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
}
double re;
// This alters the implementation of Hull et al (1994) which used a standard
// precision representation of |z|: sqrt(x*x + y*y).
// This formula should use the same definition of the magnitude returned
// by Complex.abs() which is a high precision computation with scaling.
// The checks for overflow thus only require ensuring the output of |z|
// will not overflow or underflow.
if (x > HALF && x < ROOT2) {
// x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2.
re = Math.log1p(x2y2m1(x, y)) * logOfeOver2;
} else {
// Check for over/underflow in |z|
// When scaling:
// log(a / b) = log(a) - log(b)
// So initialize the result with the log of the scale factor.
re = 0;
if (x > Double.MAX_VALUE / 2) {
// Potential overflow.
if (isPosInfinite(x)) {
// Handle infinity
return constructor.create(x, arg());
}
// Scale down.
x /= 2;
y /= 2;
// log(2)
re = logOf2;
} else if (y < Double.MIN_NORMAL) {
// Potential underflow.
if (y == 0) {
// Handle real only number
return constructor.create(log.applyAsDouble(x), arg());
}
// Scale up sub-normal numbers to make them normal by scaling by 2^54,
// i.e. more than the mantissa digits.
x *= 0x1.0p54;
y *= 0x1.0p54;
// log(2^-54) = -54 * log(2)
re = -54 * logOf2;
}
re += log.applyAsDouble(abs(x, y));
}
// All ISO C99 edge cases for the imaginary are satisfied by the Math library.
return constructor.create(re, arg());
}