On 07/29/2012 07:08 AM, Bruce Evans wrote:
On Sat, 28 Jul 2012, Stephen Montgomery-Smith wrote:
On 07/28/2012 09:31 PM, Bruce Evans wrote:
On Sat, 28 Jul 2012, Stephen Montgomery-Smith wrote:
OK. This clog really seems to work.
x*x + y*y - 1 is computed with a ULP less than 0.8. The rest of the
errors seem to be due to the implementation of log1p. The ULP of the
final answer seems to be never bigger than a little over 2.
I really don't like this version...
I can understand your reticence. I'll let you work some more on clog,
and I will concentrate on the catrig stuff.
But did you see that I messed up my earlier error analysis of log1p(1+x)?
I didn't look at it closely, but just counted the number of operations and
multiplied by 0.5 and added a safety margin to get a similar value :-).
>> Also, I don't think the problem is due to the implementation of log1p.
>> If you do an error analysis of log(1+x) where x is about exp(-1)-1,
and
>> x is correct to within 0.8 ULP, I suspect that about 2.5 ULP is the
best
>> you can do for the final answer:
>>
>> relative_error(log(1+x)) = fabs(1/((1+x) log(1+x))) *
relative_error(x)
>> = 1.58 * relative_error(x)
>>
>> Given that log1p has itself a ULP of about 1, and relative error in
x is
>> 0.8, and considering x=exp(-1)-1, this gives a ULP at around
1.58*0.8+1
>> = 2.3. And that is what I observed.
log1p(x) does expand the error signficantly here, since since exp(-1)-1 is
too close to -1 for log1p to be a good function to use on it. The
expansion in units of ulps as x -> -1 seems to be fully exponential.
was surprised by this. exp(x) also expands errors significantly, but
only linearly in x (or is it x*log(x)?).
>> (Here "=" means approximately equal to.)
>
> And I should add that I just realized that ULP isn't quite the same as
> relative error, so an extra factor of up to 2 could make its way into
> the calculations.
In fact, I think I messed it up a bunch.
So let D(f(x)) denote the absolute error in f(x).
D(f(x)) = f'(x) Dx.
So
D(log(1+x)) = Dx/(1+x).
If x is a bit bigger than exp(-1)+1 = -0.63, which has ilogb = -1. If
ULP in calculating x is around 0.8, then
Dx = 0.8 * 2^(-d-1).
where d is the number of bits in the mantissa,
So D(log(1+x)) = Dx/0.37.
1/0.37 is best spelled as `e'.
Duh! Silly me.
_______________________________________________
[email protected] mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-bugs
To unsubscribe, send any mail to "[email protected]"
