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