-
Notifications
You must be signed in to change notification settings - Fork 0
Inverse Trig & Hyp Functions
If you were looking for the page on regular Trigonometric/Hyperbolic functions, please click here.
The 6 main inverse trig functions and the 6 main inverse hyperbolic functions are implemented in this library. They're all public, static functions and are defined in the Complex class.
Unlike the regular trig and hyperbolic functions, their inverses technically have infinitely many solutions for any given input. For that reason, these implemented functions return the principal value of these functions, as defined by the National Standards of Science and Technology (NIST).
The arc sine and arc tangent (asin and atan) are defined in terms of the area hyperbolic sine and area hyperbolic tangent (asinh and atanh, the term area is used because hyp functions tie in to a hyperbola's area, but have nothing to do with it's arc length). The arc cosine, however, is not defined in terms of the area hyperbolic cosine (acosh), but is instead defined as π/2 minus the arc sine. This is how the arc cosine is officially defined by the NIST.
The asec, acsc, acot, asech, acsch, and acoth functions are defined as the acos, asin, atan, acosh, asinh, and atanh of the reciprocal of the input, respectively.
There is a page on NIST about inverse trig functions, and one on inverse hyperbolic functions. When writing this library, I've gone through the liberty of looking at the definitions provided and simplifying them in such a way that they can be expressed in only a few lines of code, without altering the branch cuts in any way. It should be noted that the definitions they provided do not contradict one another, otherwise I would be unable to do that. There are, however, two things that should be addressed:
Firstly, NIST has not officially defined these functions on the branch cuts. For instance, the inverse cosh of -1 could either be considered πi or -πi under these definitions. For this reason, you'll have to pardon me, but I've taken some liberties on these branch cuts, choosing whichever of the two definitions made sense in that specific context. The library defines acosh(-2) as ln(2+√3)+πi, when it could just as well be ln(2+√3)-πi. I'd explain my reasoning for each of the branch cuts, but honestly...for the most part, I just asked Wolfram Alpha wgat it thought the answer should be, and picked whichever option better explained the results I got.
Secondly, NIST does not actually define the asin and asinh in terms of one another, nor does it define the atan and atanh in terms of each other. However, after performing my aforementioned simplifications, I found that defining the asin in terms of the asinh and the atan in terms of the atanh yielded the exact same results as if I simply used the definitions provided (granted, I still had to take those liberties with the branch cuts).
Also, if any one happens to find a discrepancy in my formulas with respect to the NIST definitions, please don't hesitate to send me a message informing me of such. Thanks!
Complex acosh()
Returns the area hyperbolic cosine. The function first checks for two special cases, before using the default formula. In the default case, the result is simply ln(z+csgn(z)*√(z²-1)) (using the principal value of the logarithm).
The first special case is when z is a real number between -1 and 1. In this case, we return i times the Math.acos of the real part. The second special case is when the absolute square is larger than 1E18. As our input gets larger and larger, z²-1 will quickly approach the overflow limit, even when √(z²-1) isn't an overflow. However, once z.re reaches around 1E18, csgn(z)*√(z²-1) and z will be completely indistinguishable on the IEEE double floating point scale. As such, the calculation is just ln(z+z) = ln(2z) = ln(z)+ln(2).
I made sure to be careful with the branch cuts on this one, ensuring to double check that there wasn't a chance the small difference between csgn(z)*√(z²-1) and z could cause us to shift between quadrant III and quadrant IV (which would shift the logarithm by 2πi). Luckily, I found that it was impossible for such a thing to happen.
Complex asinh()
Returns the area hyperbolic sine. The function first checks for three special cases, before using the default formula. In the default case, the result is simply ln(z+√(z²+1)). However, that's not what the computer does. If this formula were applied to a large negative number, we'd run into a roundoff error, with z being equal and opposite to √(z²+1). As such, the computer instead makes use of the fact that the asinh is an odd function, and applies this formula to the secondary absolute value (abs2) of z, negating the result if z.re is negative. Thus, the return result is ln(abs2(z)+√(z²+1))*csgn(z).
The first special case is when z is an imaginary number between -i and i. In this case, we return i times the Math.asin of the imaginary part. The second special case is when the absolute square is larger than 1E18, and occurs for the same reason as it does with acosh. Once again, as z.re reaches around 1E18, √(z²+1) and abs2(z) become completely indistinguishable on the double floating point scale. As such, the calculation is just ln(abs2(z)+abs2(z))*csgn(z) = (ln(abs2(z))+ln(2))*csgn(z). Once again, I made sure that the logarithm doesn't shift over.
The third special case is when the input is too small. If we performed the formula above on 1E-100, we'd get ln(1E-100+√(1+1E-200))*csgn(1E-100). In theory, that should give us about ln(1E-100+1)*1 = ln(1+1E-100) ≈ 1E-100, which would be the correct answer. However, due to a lack of precision, the term inside the logarithm would actually be stored as just 1, and it's logarithm would come out to 0. To avoid this, whenever z has a lazy absolute value less than 1E-4, we return the first two terms of the Taylor's series expansion of asinh around 0, equal to z-z³/6.
Complex atanh()
Returns the area hyperbolic tangent. The function first checks for four special cases, before using the default formula, then checking for one more special case. In the default case, the result is simply ln((1+z)/(1-z))/2. After that, we have to check that one last special case, which is that if z is a real number greater than 1, we have to subtract πi (in implementation, it's written as ans.im=-HALFPI, since before that step the imaginary part will always be π/2). This last step is done to preserve the "odd function" property of the atanh, ensuring that atanh(-z) is always equal to -atanh(z).
The first special case that is checked is when z is either 1 or -1. In that case, we either return ∞ or -∞, respectively. The second special case is when z is infinite. In that case, we either return πi/2 if z has a positive imaginary part (or is a negative real), or we return -πi/2 if z has a negative imaginary part (or is a positive real). This is equivalent to adding πi/2*csgn(z/i), but the ternary expression uses a bit less swapping around.
The third special case is when z is an imaginary number, in which case we return i times the Math.atan of the imaginary part. The fourth special case is when the input is too small. If we perform the default case on 1E-100, we'd get ln((1+1E-100)/(1-1E-100))/2, which should in theory should give us the right answer, but in practice, a lack of precision gives us ln(1)/2 = 0. As such, we once again have to use the Taylor's series expansion around 0. Computing from the first 2 terms, this yields z+z³/3.
Complex acos()
Returns the arc cosine (in radians). There is one special case. If z is a real between -1 and 1, we just cut the **** and return the Math.acos of the real part. Otherwise, we return π/2-asin(z). Well, technically, we actually return π/2+i*asinh(zi), but that's really just to cut out the middleman.
Complex asin()
Returns the arc sine (in radians). Effectively just yields asinh(zi)/i.
Complex atan()
Returns the arc tangent (in radians). Effectively just yields atanh(zi)/i.
Complex asec()
Complex acsc()
Complex acot()
Returns the arc secant, arc cosecant, and arc cotangent, respectively. Effectively just returns the arccosine, arcsine, and arctangent (respectively) of 1/z. It should be noted than many authors disagree on the definition of the arc cotangent. Some define it as π/2 minus the arctangent, while others define it as the arctangent of 1/z. However, NIST defines the latter as the official definition.
Complex asech()
Complex acsch()
Complex acoth()
Returns the area hyperbolic secant, area hyperbolic cosecant, and area hyperbolic cotangent, respectively. Effectively just returns the acosh, asinh, and atanh (respectively) of 1/z.
(This goes into detail about the branch cuts for inverse trig/hyp functions. For information about the extension of trig/hyp functions to complex numbers, see the background section of Trig & Hyp Functions. If you want to see proof as to why these formulas work, see the next section)
All 6 main trig functions are periodic. If you plug in any of them for a given value of x, then plug the same function in for x+2π, you'll get the same output. In fact, the same holds if you add any integer multiple of 2π (0, 2π, -2π, 4π, 100π, 1000π, etc.) As such, if you wanted to find which value for x, when fed into a trig function, would yield a specific value y, you'd get infinitely many solutions. However, they're all basically the same solution, just plus some multiple of 2π. That is to say, if you know one of them, you know all of them. So we really only need to know one of them, and we decided pretty early on which value was the preferred one.
In general, whichever value was closest to 0 was the preferred choice. For odd functions such as sine, tangent, cosecant, and cotangent, this yielded some value between -π/2 and π/2 inclusive. However, for the even functions, cosine and secant, this wasn't enough, as there'd always be a positive answer, and an equal and opposite negative answer. As such, for those functions, we defined the result to be whichever result was positive and closest to 0. A consequence of this is that the arc cosine plus arc sine is always π/2.
As we discovered/invented the hyperbolic functions, we now also had a new set of inverse functions with infinitely many solutions. This time, though, there was usually only one or two real solutions (the rest would be those solutions plus some multiple of 2πi). And when there are two real solutions, we can just pick whichever one is positive, just like last time. Of course, there are some inputs which have no real solutions, like for the tanh. These are a bit trickier.
In fact, as we learned more and more about complex numbers, we found ourselves in a position where there was no longer an official definition for the inverse anything of anything. And, to be fair, we still haven't resolved that issue. Why's that? Well, we really only need official definitions for real inputs. If anybody runs into a situation where they need to find the inverse tangent of, say, 3+7i, they're probably in a situation where either all solutions are equally valid, or where they themselves are better suited to decide where the branch cuts are. There will never be one official definition that equally suits the needs of all mathematical problems. If anything, someone in that situation might actually be better off using the logarithmic definitions of the arc tangent than the actual function itself, since the logarithm has more predictable branch cuts that are easier to manage.
With that said, when I wrote this code, I had to use a definition that made sense. And even though the values along the branch cuts were tricky, I should probably be lucky that the NIST had an official definition for these functions at all, much less whether or not they were defined on the branch cuts.
cosh(t)=(e^t+e^-t)/2
sinh(t)=(e^t-e^-t)/2
cosh(t)+sinh(t)=e^t
t = ln(cosh(t)+sinh(t))
Let's say x=cosh(t), and we're trying to solve for t
If we assume that t>=0, then we know that t=acosh(x) (t must be positive, because the acosh, like the acos, has to return a positive value)
acosh(x) = ln(x + sinh(acosh(x)))
cosh²(t)-sinh²(t)=1
sinh²(t)=cosh²(t)-1
t>=0, so sinh(t)>=0
sinh(t)=√(cosh²(t)-1) = √(x²-1)
acosh(x) = ln(x+√(x²-1))
Now let's go back to our original formula, but this time use it to solve for asinh.
t=ln(cosh(t)+sinh(t))
Let's say y=sinh(t). t no longer strictly has to be positive.
t=asinh(y)
asinh(y)=ln(cosh(asinh(y))+y)
cosh²(t)-sinh²(t)=1
cosh²(t)=sinh²(t)+1
cosh(t) is always positive, no matter what, as long as t is real.
cosh(t)=√(sinh²(t)+1)=√(y²+1)
asinh(y)=ln(y+√(y²+1))
So now we have a definition for acosh and asinh. What about atanh?
Well, there's a simple solution. Let's look back at our definitions for cosh and sinh. We saw that if we add the two together, we get e^t. But, if we subtract the two, we get
cosh(t)-sinh(t) = e^-t.
Likewise, that means (cosh(t)+sinh(t))/(cosh(t)-sinh(t)) = e^t/e^-t = e^(2t).
If we multiple the numerator and denominator by sech(t), we get
(1+tanh(t))/(1-tanh(t)) = e^(2t)
Take the log of both sides
ln((1+tanh(t))/(1-tanh(t))) = 2t
t = ln((1+tanh(t))/(1-tanh(t)))/2
If we define x=tanh(t), and likewise that t=atanh(x), we get
atanh(x) = ln((1+x)/(1-x))/2
So the three main inverse hyperbolic functions can be written as:
acosh(x) = ln(x+√(x²-1))
asinh(x) = ln(x+√(x²+1))
atanh(x) = ln((1+x)/(1-x))/2
It should be noted that if each square root was replaced with the negative square root, while we would no longer have a principal solution, we would still get a perfectly valid solution. For the definition I implemented for acosh, I had it use the negative square root whenever the real part was negative, which helps keep the function continuous along the complex plane (and helps it to fall in line with NIST definitions).
Since cosine, sine, and tangent can be defined in terms of the hyperbolic functions, that also gives us an easy way to define the inverse trig functions. Since
cos(x) = cosh(xi)
sin(x) = sinh(xi)/i
tan(x) = tanh(xi)/i
That means
acos(x) = acosh(x)/i
asin(x) = asinh(xi)/i
atan(x) = atanh(xi)/i
And if we plug it all in, we get:
acos(x) = -i*ln(x+√(x²-1))
asin(x) = -i*ln(xi+√(1-x²))
atan(x) = -i*ln((1+xi)/(1-xi))/2
Which should make sense, if you think about what's actually happening here. The first two equations just boil down to assembling Euler's identity e^(Θi) = cos(Θ)+sin(Θ)i, then taking the natural log and dividing by i. The third one assembles Euler's identity (but multiplied by the secant), then divides by another instance of Euler's identity (still multiplied by the secant, but this time conjugated). What we get is e^(Θi)/e^(-Θi) = e^(2Θi). We then take the log, divide by i, split it in half, and voila.
It's also common to see the atan written as
atan(x) = i*ln((1-xi)/(1+xi))/2
The asec, acsc, acot, etc. are all defined as these same things, but with the input inverted. Of course, there are a few extra steps to these functions that were added upon implementation just to make the branch cuts line up with the official standards, but it's still just another one of the valid solutions.
So, I already mentioned that once you have a single solution to these functions, you have all the solutions. But I still have yet to tell you how to get the rest of the solutions. Well...here you go...
Note that N, from henceforth, will represent any real integer.
acos(x): ±acos(x)+2Nπ
asin(x): asin(x)+2Nπ, -asin(x)+(2N+1)π
atan(x): atan(x)+Nπ
asec(x): ±asec(x)+2Nπ
acsc(x): acsc(x)+2Nπ, -acsc(x)+(2N+1)π
acot(x): acot(x)+Nπ
acosh(x): ±acosh(x)+2Nπi
asinh(x): asinh(x)+2Nπi, -asinh(x)+(2N+1)πi
atanh(x): atanh(x)+Nπi
asech(x): ±asech(x)+2Nπi
acsch(x): acsch(x)+2Nπi, -acsch(x)+(2N+1)πi
acoth(x): acoth(x)+Nπi