[Mpc-discuss] asin, atan, etc
Stephen Montgomery-Smith
stephen at missouri.edu
Sun Aug 12 19:58:06 CEST 2012
I was using mpc to check a casinh that I am writing for FreeBSD. I
noticed that your asinh functions can take a very long time to compute
when for certain data. I was wondering if you guys would be interested
in using the same algorithm that I used.
You guys effectively use asinh(z) = log(z + sqrt(z^2+1)), and if z=x+I*y
where x is very small and |y| <= 1, then this computation may have to be
performed in double the precision you actually want. And my sample data
was something like x=1e-3000, y=1-1e-3000.
I used the algorithm that can be found in "Implementing the complex
arcsine and arccosine functions using exception handling" by T. E. Hull,
Thomas F. Fairgrieve, and Ping Tak Peter Tang, published in ACM
Transactions on Mathematical Software, Volume 23 Issue 3, 1997, Pages
299-335, http://dl.acm.org/citation.cfm?id=275324.
It really is a very nice algorithm, and because you guys don't have
problems with underflow or overflow, you don't need to mess with the
exception handling.
I have described the algorithm at
http://people.freebsd.org/~stephen/catrig.c, but the part you need is
completely described in the opening comments, because underflow and
overflow is not a problem for you guys.
Also, regarding atan, instead of using the formula:
Im(atan(x+I*y)) = 1/4 * [log(x^2+(1+y)^2) - log (x^2 +(1-y)^2)]
you might like to use the formula
Im(atan(x+I*y)) = 1/4 * [log1p (4y / (x^2 +(1-y)^2))]
when y > 0, and then use the fact that atan is odd for the case y<0.
Stephen
More information about the Mpc-discuss
mailing list