evalf()/Digits bugs

[Bob McElrath]

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