Date: Fri, 03 Aug 2012 15:35:41 -0500 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: freebsd-numerics@freebsd.org Subject: Re: Complex arg-trig functions Message-ID: <501C361D.4010807@missouri.edu> In-Reply-To: <5017111E.6060003@missouri.edu> References: <5017111E.6060003@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?501C361D.4010807>