Richard B. Kreckel [kreckel@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@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@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@students.wisc.edu) Univ. of Wisconsin at Madison, Department of Physics