evalf()/Digits bugs
Bob McElrath
mcelrath at draal.physics.wisc.edu
Mon Jan 7 16:06:10 CET 2002
Richard B. Kreckel [kreckel at zino.physik.uni-mainz.de] wrote:
> Hi,
>
> On Sun, 6 Jan 2002, Bob McElrath wrote:
> > Some interesting (and wrong) calculations, probably rounding somewhere:
> >
> > (1)<mcelrath at hawk:/u/mcelrath> ginsh
> > ginsh - GiNaC Interactive Shell (GiNaC V1.0.3)
> > __, _______ Copyright (C) 1999-2001 Johannes Gutenberg University Mainz,
> > (__) * | Germany. This is free software with ABSOLUTELY NO WARRANTY.
> > ._) i N a C | You are welcome to redistribute it under certain conditions.
> > <-------------' For details type `warranty;'.
> >
> > Type ?? for a list of help topics.
> > > evalf(1/3);
> > 0.33333333333333333334
> > > Digits=30;
> > 30
> > > evalf(1/3);
> > 0.333333333333333333333333333333333333334
> > > Digits=40;
> > 40
> > > evalf(1/3);
> > 0.3333333333333333333333333333333333333333333333334
> > > Digits=50;
> > 50
> > > evalf(1/3);
> > 0.33333333333333333333333333333333333333333333333333333333336
> > >
>
> Objection, your Honor!
>
> You ordered 30 decimal digits and you got 39, of which the last one is
> incorrect. The 30 decimal digits you ordered are all correct. This
> pattern repeats itself.
>
> Paul Zimmermann is currently thoroughly investigating precision. He is
> comparing several packages (Pari, CLN, MPFR) and last I heard from
> him (about two weeks ago) he has not found any "errors" (according to
> above definition) yet.
>
> This is because CLN doesn't bother dealing with bits but only with
> machine-sized words -- for reasons that are more than obvious. So you
> always get some extra precision and there is just no code that throws
> away that extra precision.
>
> One of the cases closest to an actual error is for Digits=27, where you
> get 0.333333333333333333333333333335, which is actually 30 Digits. So
> it's still correct.
>
> Convinced?
> [ ] Yes, this sounds reasonable.
> [ ] No, I'll write a patch for GiNaC that rounds away extra digits in
> the output.
Well, it sounds reasonable, but is counter-intuitive. Consider Maple:
(0)<mcelrath at draal:/home/mcelrath> maple
|\^/| Maple V Release 5 (WMI Campus Wide License)
._|\| |/|_. Copyright (c) 1981-1997 by Waterloo Maple Inc. All rights
\ MAPLE / reserved. Maple and Maple V are registered trademarks of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
> evalf(1/3);
.3333333333
> Digits:=30;
Digits := 30
> evalf(1/3);
.333333333333333333333333333333
Here it not only returns the number of digits requested, but does not show odd
rounding errors. I'll see about a patch.
I worry a little bit that after some computations, your answer wouldn't be
accurate to Digits accuracy because of machine-word rounding effects.
Cheers,
-- Bob
Bob McElrath (rsmcelrath at students.wisc.edu)
Univ. of Wisconsin at Madison, Department of Physics
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 240 bytes
Desc: not available
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20020107/19d0b917/attachment.pgp
More information about the GiNaC-list
mailing list