in commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java [2029:2143]
private static Complex acos(final double real, final double imaginary,
final ComplexConstructor constructor) {
// Compute with positive values and determine sign at the end
final double x = Math.abs(real);
final double y = Math.abs(imaginary);
// The result (without sign correction)
final double re;
final double im;
// Handle C99 special cases
if (isPosInfinite(x)) {
if (isPosInfinite(y)) {
re = PI_OVER_4;
im = y;
} else if (Double.isNaN(y)) {
// sign of the imaginary part of the result is unspecified
return constructor.create(imaginary, real);
} else {
re = 0;
im = Double.POSITIVE_INFINITY;
}
} else if (Double.isNaN(x)) {
if (isPosInfinite(y)) {
return constructor.create(x, -imaginary);
}
// No-use of the input constructor
return NAN;
} else if (isPosInfinite(y)) {
re = PI_OVER_2;
im = y;
} else if (Double.isNaN(y)) {
return constructor.create(x == 0 ? PI_OVER_2 : y, y);
} else {
// Special case for real numbers:
if (y == 0 && x <= 1) {
return constructor.create(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
}
final double xp1 = x + 1;
final double xm1 = x - 1;
if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
final double yy = y * y;
final double r = Math.sqrt(xp1 * xp1 + yy);
final double s = Math.sqrt(xm1 * xm1 + yy);
final double a = 0.5 * (r + s);
final double b = x / a;
if (b <= B_CROSSOVER) {
re = Math.acos(b);
} else {
final double apx = a + x;
if (x <= 1) {
re = Math.atan(Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))) / x);
} else {
re = Math.atan((y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))) / x);
}
}
if (a <= A_CROSSOVER) {
final double am1;
if (x < 1) {
am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
} else {
am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
}
im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
} else {
im = Math.log(a + Math.sqrt(a * a - 1));
}
} else {
// Hull et al: Exception handling code from figure 6
if (y <= (EPSILON * Math.abs(xm1))) {
if (x < 1) {
re = Math.acos(x);
im = y / Math.sqrt(xp1 * (1 - x));
} else {
// This deviates from Hull et al's paper as per
// https://svn.boost.org/trac/boost/ticket/7290
if ((Double.MAX_VALUE / xp1) > xm1) {
// xp1 * xm1 won't overflow:
re = y / Math.sqrt(xm1 * xp1);
im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
} else {
re = y / x;
im = LN_2 + Math.log(x);
}
}
} else if (y <= SAFE_MIN) {
// Hull et al: Assume x == 1.
// True if:
// E^2 > 8*sqrt(u)
//
// E = Machine epsilon: (1 + epsilon) = 1
// u = Double.MIN_NORMAL
re = Math.sqrt(y);
im = Math.sqrt(y);
} else if (EPSILON * y - 1 >= x) {
re = PI_OVER_2;
im = LN_2 + Math.log(y);
} else if (x > 1) {
re = Math.atan(y / x);
final double xoy = x / y;
im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
} else {
re = PI_OVER_2;
final double a = Math.sqrt(1 + y * y);
im = 0.5 * Math.log1p(2 * y * (y + a));
}
}
}
return constructor.create(negative(real) ? Math.PI - re : re,
negative(imaginary) ? im : -im);
}