Revisiting inverse, complex, hyperbolic, floating-point trig functions
I thought I had learned everything worth knowing about trigonometry years ago. But I've just been reading a paper by W. Kahan, “Branch Cuts for Complex Elementary Functions” (in Iserles and Powell (eds.), State of the Art of Numerical Analysis, 1987). Kahan is a professor at Berkeley, the father of the IEEE floating-point standard, and a consultant in the design of algorithms used in Hewlett Packard calculators. His paper reminded me of some important facts about computing trig functions.
Let's start with the hyperbolic cosine,
cosh(t) = (et + e-t)/2
This is the only trig function I've ever seen in Newsweek, where it was in an article about the Gateway Arch in St. Louis. It's the equation of an upside down arch, as well as the catenary, the curve formed by a chain hanging between two supports. I am really interested in the inverse function, acosh(x). To find its formula, we need to solve the equation
x = (et + e-t)/2
for t. An intermediate step is the quadratic in et,
e2t - 2xet + 1 = 0
We solve the quadratic and then take a logarithm to get
acosh(x) = log(x ± (x2 - 1)1/2)
This formula for acosh(x) can be found in almost any calculus or math reference book. But it turns out to be unsuitable for actual computation, except for a limited range of x. As a general-purpose computational formula, it has at least three defects:
We will look at each of these problems in a moment, but first let's look at that ± sign. For real t, the graph of x = cosh(t) never drops below x = 1, so, for now, we can assume x ≥ 1 and get a real square root term. The ± sign gives two possible values for the argument of log. But the arguments turn out to be reciprocals of each other, so the resulting logs are negatives of each other. This corresponds to the two legs of the double-valued function obtained by turning the St. Louis arch on its side. In principle, we could use either sign in a computation.
Now, let's look at the various computational problems with the textbook formula. Overflow is the least important, but the easiest to see. Say x is greater than 10154, but less than 10308. This is pretty extreme, but we still should be able to compute acosh(x) because it's essentially equal to log(2x), which is between roughly 355 and 710. However, the formula squares x, producing an intermediate result larger than the largest number allowed in IEEE floating-point double precision format. This particular problem is easily fixed with some scaling, but that doesn't help the more serious problems.
The accuracy problem shows up in the real case if we choose the minus sign in the formula. Let
x = cosh(10)
which is about 1.1013 • 104. Using the minus sign, we expect acosh(x) to be very close to -10. The trouble comes when we compute x - (x2 - 1)1/2. The result is about 4.54 • 10-5, but it is obtained with serious subtractive cancellation. Taking the logarithm produces -10.00000001350353. We've lost almost half of the digits.
The cancellation problem gets worse quickly. For all values of x greater than about cosh(18.8), the computed value of (x2 - 1)1/2 is actually equal to x, the subtraction produces zero, and the logarithm becomes negative infinity. For real x between 1 and the square root of overflow, the cancellation difficulty can be avoided by using the plus sign in the textbook formula. But that just transfers the problem to negative x and doesn't make the formula reliable in general.
Now let's consider complex arithmetic and use z in place of x. As functions of a complex argument, both log(z) and z1/2 have discontinuous jumps when z crosses the branch cut, which is conventionally taken to be the negative real axis. The jump in log(z) is always 2i, while the jump in z1/2 is |z|1/2.
All of the inverse trig functions must have branch cuts somewhere in the complex plane, but I had forgotten (or never knew) where they were supposed to be. The branch cuts for most of the inverse trig functions are symmetric with respect to the origin. But the branch cut for acosh(z) is unusual. Near z = 1, acosh(z) behaves like (z - 1)1/2. For large negative real z, acosh(z) behaves like log(z). So the branch cut for acosh(z) consists of the portion of the real axis to the left of +1.
Since the branch cut for acosh(z) is a portion of the real axis, we expect the function to be perfectly well behaved near the imaginary axis. But the textbook formula involves (z2 - 1)1/2. This maps the imaginary axis onto the negative real axis and introduces spurious jumps. Try using the formula to compute acosh(z) for z = 2i ± .0001. Both values should be close to 1.4436 + 1.5709i. But we get the wrong sign when z is in the left half plane.
Kahan's solution to all these difficulties is to use a better formula.
acosh(z) = 2 log(((z + 1)/2)1/2 + ((z - 1)/2)1/2)
This involves two square roots, but it is still faster than the textbook formula modified to overcome all its deficiencies. We will incorporate this, and some of Kahan's algorithms for other trig functions, into MATLAB 5.3, available later this year.
As I was writing this column, a colleague told me about the Web page of John Richardson. (The Web page no longer exists, CBM, 2/17/04.) Richardson describes an interesting way to visualize complex functions by mapping complex numbers to color values. Zero is mapped to black, large values to white, real values to shades of red, and imaginary values to shades of blue and green. I've used his approach to produce graphs of the inverse trig functions and have added contour lines for the real and imaginary parts of the function values. The heavy black lines are the branch cuts.