in commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java [2732:2896]
private static Complex atanh(final double real, final double imaginary,
final ComplexConstructor constructor) {
// Compute with positive values and determine sign at the end
double x = Math.abs(real);
double y = Math.abs(imaginary);
// The result (without sign correction)
double re;
double im;
// Handle C99 special cases
if (Double.isNaN(x)) {
if (isPosInfinite(y)) {
// The sign of the real part of the result is unspecified
return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
}
// No-use of the input constructor.
// Optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
return NAN;
} else if (Double.isNaN(y)) {
if (isPosInfinite(x)) {
return constructor.create(Math.copySign(0, real), Double.NaN);
}
if (x == 0) {
return constructor.create(real, Double.NaN);
}
// No-use of the input constructor
return NAN;
} else {
// x && y are finite or infinite.
// Check the safe region.
// The lower and upper bounds have been copied from boost::math::atanh.
// They are different from the safe region for asin and acos.
// x >= SAFE_UPPER: (1-x) == -x
// x <= SAFE_LOWER: 1 - x^2 = 1
if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
// Normal computation within a safe region.
// minus x plus 1: (-x+1)
final double mxp1 = 1 - x;
final double yy = y * y;
// The definition of real component is:
// real = log( ((x+1)^2+y^2) / ((1-x)^2+y^2) ) / 4
// This simplifies by adding 1 and subtracting 1 as a fraction:
// = log(1 + ((x+1)^2+y^2) / ((1-x)^2+y^2) - ((1-x)^2+y^2)/((1-x)^2+y^2) ) / 4
//
// real(atanh(z)) == log(1 + 4*x / ((1-x)^2+y^2)) / 4
// imag(atanh(z)) == tan^-1 (2y, (1-x)(1+x) - y^2) / 2
// imag(atanh(z)) == tan^-1 (2y, (1 - x^2 - y^2) / 2
// The division is done at the end of the function.
re = Math.log1p(4 * x / (mxp1 * mxp1 + yy));
// Modified from boost which does not switch the magnitude of x and y.
// The denominator for atan2 is 1 - x^2 - y^2.
// This can be made more precise if |x| > |y|.
final double numerator = 2 * y;
final double denominator;
if (x < y) {
final double tmp = x;
x = y;
y = tmp;
}
// 1 - x is precise if |x| >= 1
if (x >= 1) {
denominator = (1 - x) * (1 + x) - y * y;
} else {
// |x| < 1: Use high precision if possible:
// 1 - x^2 - y^2 = -(x^2 + y^2 - 1)
// Modified from boost to use the custom high precision method.
denominator = -x2y2m1(x, y);
}
im = Math.atan2(numerator, denominator);
} else {
// This section handles exception cases that would normally cause
// underflow or overflow in the main formulas.
// C99. G.7: Special case for imaginary only numbers
if (x == 0) {
if (imaginary == 0) {
return constructor.create(real, imaginary);
}
// atanh(iy) = i atan(y)
return constructor.create(real, Math.atan(imaginary));
}
// Real part:
// real = Math.log1p(4x / ((1-x)^2 + y^2))
// real = Math.log1p(4x / (1 - 2x + x^2 + y^2))
// real = Math.log1p(4x / (1 + x(x-2) + y^2))
// without either overflow or underflow in the squared terms.
if (x >= SAFE_UPPER) {
// (1-x) = -x to machine precision:
// log1p(4x / (x^2 + y^2))
if (isPosInfinite(x) || isPosInfinite(y)) {
re = 0;
} else if (y >= SAFE_UPPER) {
// Big x and y: divide by x*y
re = Math.log1p((4 / y) / (x / y + y / x));
} else if (y > 1) {
// Big x: divide through by x:
re = Math.log1p(4 / (x + y * y / x));
} else {
// Big x small y, as above but neglect y^2/x:
re = Math.log1p(4 / x);
}
} else if (y >= SAFE_UPPER) {
if (x > 1) {
// Big y, medium x, divide through by y:
final double mxp1 = 1 - x;
re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
} else {
// Big y, small x, as above but neglect (1-x)^2/y:
// Note: log1p(v) == v - v^2/2 + v^3/3 ... Taylor series when v is small.
// Here v is so small only the first term matters.
re = 4 * x / y / y;
}
} else if (x == 1) {
// x = 1, small y:
// Special case when x == 1 as (1-x) is invalid.
// Simplify the following formula:
// real = log( sqrt((x+1)^2+y^2) ) / 2 - log( sqrt((1-x)^2+y^2) ) / 2
// = log( sqrt(4+y^2) ) / 2 - log(y) / 2
// if: 4+y^2 -> 4
// = log( 2 ) / 2 - log(y) / 2
// = (log(2) - log(y)) / 2
// Multiply by 2 as it will be divided by 4 at the end.
// C99: if y=0 raises the ‘‘divide-by-zero’’ floating-point exception.
re = 2 * (LN_2 - Math.log(y));
} else {
// Modified from boost which checks y > SAFE_LOWER.
// if y*y -> 0 it will be ignored so always include it.
final double mxp1 = 1 - x;
re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
}
// Imaginary part:
// imag = atan2(2y, (1-x)(1+x) - y^2)
// if x or y are large, then the formula:
// atan2(2y, (1-x)(1+x) - y^2)
// evaluates to +(PI - theta) where theta is negligible compared to PI.
if (x >= SAFE_UPPER || y >= SAFE_UPPER) {
im = Math.PI;
} else if (x <= SAFE_LOWER) {
// (1-x)^2 -> 1
if (y <= SAFE_LOWER) {
// 1 - y^2 -> 1
im = Math.atan2(2 * y, 1);
} else {
im = Math.atan2(2 * y, 1 - y * y);
}
} else {
// Medium x, small y.
// Modified from boost which checks (y == 0) && (x == 1) and sets re = 0.
// This is same as the result from calling atan2(0, 0) so exclude this case.
// 1 - y^2 = 1 so ignore subtracting y^2
im = Math.atan2(2 * y, (1 - x) * (1 + x));
}
}
}
re /= 4;
im /= 2;
return constructor.create(changeSign(re, real),
changeSign(im, imaginary));
}