Complex arg-trig functions

Stephen Montgomery-Smith stephen at missouri.edu
Fri Aug 3 20:35:47 UTC 2012


On 07/30/12 17:56, Stephen Montgomery-Smith wrote:
> I will be posting my updates on this web page:
>
> http://people.freebsd.org/~stephen/
>
> The file catrig.c contains implementations of cacos(h), casin(h),
> catan(h).  I have been working on it everyday - I seem unable to put it
> down.

I know you guys are busy and haven't had time to look this over.  But in 
many ways this is a good thing, as I am constantly finding problems with 
the code.

One of the hardest things is trying to avoid unwarranted underflows. 
For example, when |z| is large (bigger than 1e20) I use the approximation

catanh(z) = 1/z + cpack(0, copysign(M_PI_2, cimag(z)))

If the real part of 1/z underflows, then this is a legitimate underflow. 
  But if the imaginary part of 1/z underflows, then this is not 
legitimate, because then this really small number is being added to plus 
or minus I*PI/2, and it should just disappear into rounding error.

I found this out by testing the underflow flag for lots of random 
inputs.  I ended up writing my own code to compute the real part of 1/z 
without computing the imaginary part.  I have a working solution, but 
I'm sure it isn't optimal.  (I see there is code in the IEEE standard 
PDF document which I might copy, or I might copy ideas from the hypot 
function.)

Another problem I experienced - I sometimes need to compute log(0.5*y) 
when x is extremely small.  I was shocked to find that this produced 
spurious underflows when x was very close to DBL_MIN.  Of course with 
hindsight it is obvious why it happens, but these errors are hard to 
think of in advance, and hard to track down.

Right now I think I have working code.  But I also thought I had working 
code last week.  So who knows.

Anyway, thanks for reading my ramblings.  I do find this stuff 
interesting, so if anyone wants to reply and continue the conversation, 
please do so.

Stephen



More information about the freebsd-numerics mailing list