Use of C99 extra long double math functions after r236148
Stephen Montgomery-Smith
stephen at missouri.edu
Sun Aug 12 23:11:33 UTC 2012
Three things.
1. I think I now understand the casinh algorithm by Hull et al rather
well. I implemented it in the attached code. Bruce, could you try it
out against pari? It did compare rather well to Mathematica's ArcSinh
function. I tested it mostly near the branch cut ends at I and -I.
(The branch cut's behavior is somewhat similar to that of sqrt.)
After this, I think getting casin, cacos and cacosh should be relatively
straightforward. The hardest thing about cacos and cacosh will be
figuring out where all the branch cuts and discontinuities are.
I also worked quite hard to avoid underflows and overflows from
happening. But I need to go through it again several times to check
that the logic really works.
2. I have thought more about the problem of computing clog(z) when |z|
is close to 1. I now think it might even require precision that is 3
times better than double precision. It is possible that you could, by
chance, pick x and y so that |z| = 1 + 1e30. (I wrote a test program,
and this did actually happen a few times.) So when you compute
x^2+y^2-1, if you do it using double double precision, you will get
2e30, but only the most significant digits will be accurate. The ULP
will be about 1e15. You need triple double precision to get ULP's close
to 1.
And I cannot even figure out how to do log(1+z) when z is close to 0!
The trouble is, (1+x)^2+y^2-1 = 2x+x^2+y^2, and if x is negative and
approximately -y^2, you are in the situation of subtracting nearly equal
numbers! (Obviously if z is very close to 0, I could use Taylor's series.)
3. I haven't learned proper style yet. (Is it what you call KNF?) I
have always had a distrust of consistency, especially when it comes to
people's coding styles. Sorry, but this is going to be a bit painful
for me.
Stephen
More information about the freebsd-numerics
mailing list