Date: Sun, 29 Jul 2012 22:08:45 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-bugs@FreeBSD.org, Stephen Montgomery-Smith <stephen@FreeBSD.org>, Bruce Evans <brde@optusnet.com.au> Subject: Re: bin/170206: complex arcsinh, log, etc. Message-ID: <20120729195931.J2598@besplex.bde.org> In-Reply-To: <5014A284.2060204@missouri.edu> References: <201207281550.q6SFoBdC073014@freefall.freebsd.org> <20120729121443.M1008@besplex.bde.org> <5014A284.2060204@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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'. I think exp(N) is right for exp(-N)+1.
I didn't think about the scaling of Dx.
> Since log(1+x) is a little bit bigger than -1, and so ilogb(log(1+x)) = -1.
>
> ULP(log(1+x)) = Dx/0.37 * 2^{d+1} = 0.8/0.37 = 2.2
Certainly there is exponential scaling. It seems to be fully exponential,
unlike for exp(x).
> Now add 1 for ULP in calculating log1p, and this only gives a ULP of
> 3.2. So the observed bound is actually better than expected. If one
> could get the ULP of log1p to be as good as possible (0.5), the best ULP
> one could get is 2.7. We still do a bit better than that.
This is easier to think about in binary with log2(). Consider values
-1+2**-N. We can't get closer to -1 than with N = d (maybe off by 1).
log2(1+this) is -d. If we are off by a full ulp in the arg, and towards
0, then the wrong arg is -1+2**-d+2**-d = -1+2**-(d-1). log2(1+this) is
-(d-1). The error has expanded by 1 ulp to a lot. Not by 2*d though
(in ulps). Only the relative error expanded by about 2**d (actually
by 2**d/d). It remains to calculate the error between -d and -(d-1) in
ulps. For d = 53, log2(53) is between 5 and 6, so all except the top
5 or 6 bits are wrong, so the error is between 2**47 and 2**48 ulps.
The details are more complicated for N = 1, and I didn't do them.
Your formula gives the above when modified to:
ULP(log(1+x)) = 1.0 * (exp2(N)/N)
Somehow the relative error us the same as the error in ulps (and the
approximate expansion factor).
This shows that we shouldn't consider using log1p() near -1, and I
think -1+exp(-1) is too near. x above -1/2 or -1/4 should be OK.
log1p() on -1+exp(-1) seems to only happen in your version when
I can't see any version does this. The worse case in your version
seems to be when ax*ax+ay*ay == 0.7. Then subtracting 1 gives -0.3.
My version does the extra precision calculation of ax*ax+ay*ay in
many more cases, but its worst case for using log1p() seems to be
when subtracting 1 gives -0.5. This may be too near -1. I wouldn't
like even a factor of 2 expansion in the error. But errors of a
couple of ulps will occur for different reasons using log() instead
of log1p(). Another thing to do is to find exactly when log1p()
starts being worse than log().
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120729195931.J2598>