Use of C99 extra long double math functions after r236148

Stephen Montgomery-Smith stephen at missouri.edu
Sun Aug 12 23:09:45 UTC 2012


On 07/19/2012 01:12 PM, Bruce Evans wrote:
> On Thu, 19 Jul 2012, Stephen Montgomery-Smith wrote:
>
>> I have the ULP down to about 1.2 now.  I don't see how I can do
>> better, because I have to invoke log functions twice, and probably
>> each one has a ULP of about 0.6.
>>
>> Also I decided to use 1/2 log(x*x+y*y) when x and y are not too large.
>
> That's close to Apple complex.c clog().  Once you don't use hypot(),
> it is clearly best to use log1p():
>
>      log(sqrt(x*x + y*y)) = log(|x|) + 1/2 log(1 + (y*y)/(x*x))
>                   = log(|x|) + 1/2 log1p((y*y)/(x*x))
>
> where |x| >= |y| so that log1p()'s arg is as small as possible.
>
>> I am really rather proud of how I got around the large ULP when
>> hypot(x,y) is close to 1.  I would be glad if any of you could look at
>> the code when you get a chance.
>
> WIll look more closely later.  I see that you already use log1p() and much
> more.  Apple clog() uses not so much more, mainly by depending on extra
> precision in hardware.  The above also avoids overflow and use of hypot()
> for all finite x and and y, but is probably too simple.

I think their solution merely avoids the overflow/underflow problem, and 
was not meant to address the problem I worked on.  However, their 
solution will fail if y=1e100 and y=1e-100.

This caused me to realize that my solution failed to account for 
underflow, so here is my next iteration.




More information about the freebsd-numerics mailing list