On Mon, 4 Feb 2013 13:35:57 +0100 Lionel Cons wrote:
> On 1 February 2013 16:09, David Korn <[email protected]> wrote:
> > Subject: Re: [ast-developers] Floating point oddities with pow()? Accuracy
> > problem?
> > --------
> >>
> >> Can someone explain why the the examples below yield a difference? IMO
> >> pow() should deliver the same result as a x*x*x*... sequence. Or not?
> >>
> >> ksh -c 'float x y ; ((x=y=.1 , y=pow(y,4.) , x=x*x*x*x)) ; printf
> >> "%30.30f\n" x-y'
> >> -0.000000000000000000000006617444
> >> ksh -c 'float x y ; ((x=y=.1 , y=pow(y,6.) , x=x*x*x*x*x*x)) ; printf
> >> "%30.30f\n" x-y'
> >> -0.000000000000000000000000103397
> >> ksh -c 'float x y ; ((x=y=.1 , y=pow(y,8.) , x=x*x*x*x*x*x*x*x)) ;
> >> printf "%30.30f\n" x-y'
> >> -0.000000000000000000000000001615
> >>
> >> Platform is Suse 12.1 on AMD64.
> >>
> >> Lionel
> >
> > The AMD64 has only 18 significant digits with long double and you are asking
> > for 30. This is caused by floating point round off.
> That was not the problem I was referring to. The problem I have is
> that ksh93 builtin pow() and a pow made from individual
> multiplications gives slightly different results, many digits after
> the comma.
the example ksh loop compares 2 completely different implementations of the
same operations
pow(v,N) is presumably intimate with the floating point implementation and
avoids the roundoff errors associated with the N multiplies used in the rolled
out form
similar divergence shows up in the C equivalent code below
run it with
./a.out # displays the floating point implementation LDBL_DIG
digits
./a.out 30 # displays 30 digits
aside from convenience, roundoff error is the reason why we have functions like
pow() and exp1m()
back in the day sliderules made roundoff error obvious
(because our eyes were more involved)
we now sometimes expect too much from our machines
when they are just doing what we instructed
--
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
int
main(int argc, char** argv)
{
int i, d;
long double x, y, z;
d = *++argv ? atoi(*argv) : LDBL_DIG;
for (i = -50, z = -5; i <= 50; i++, z += 0.1)
{
x = y = 0.1 * z;
x = x * x * x * x * x * x * x * x;
y = powl(y, 8);
printf("%*.*Lf %*.*Lf %+*.*Lf\n", d, d, x, d, d, y, d, d, x -
y);
}
return 0;
}
--
_______________________________________________
ast-developers mailing list
[email protected]
http://lists.research.att.com/mailman/listinfo/ast-developers
