evalf()/Digits bugs

[Richard B. Kreckel]

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.

> Note also that this is on page 14 of the tutorial, and the tutorial does not
> have the unusual final digit...

Hmm, we ought to fix the tutorial then...   :)

Cheers
    -richy.
-- 
Richard B. Kreckel
<Richard.Kreckel at Uni-Mainz.DE>
<http://wwwthep.physik.uni-mainz.de/~kreckel/>