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