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
Note also that this is on page 14 of the tutorial, and the tutorial does not have the unusual final digit... Cheers, -- Bob Bob McElrath (rsmcelrath@students.wisc.edu) Univ. of Wisconsin at Madison, Department of Physics
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.
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@Uni-Mainz.DE> <http://wwwthep.physik.uni-mainz.de/~kreckel/>
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
Hi, On Mon, 7 Jan 2002, Bob McElrath wrote:
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.
Oh, if you worry about that, then throwning away digits in the output is not going to help since it is merely cosmetic. If we could have two digits, then compared to 3.3, you are saying that 3.333 is better than 3.3334. Why bother? If you are doing longish computations and at the end your result is either NaN or zero or -10^5 while it should have been unity is something that falls under the "shit happens" category. It may happen with double- precision and proper IEEE rounding, don't forget that! The only way out seems to be interval arithmetic, but that is at least twice as costly and can only be defined rigorously on ring-operations. Another alternative would be to attach the precision to which any floating point object is correct to that object and correctly handle that precision. But that does not work either. Try this in Mathematica: In[1]:= x=1.11111111111111111111; In[2]:= Do[x=2*x-x,{50}]; In[3]:= If[x==x+1,"Mathematica fundamentally alters the mechanics of mathematics"] You'll see a funny thing happen. (BTW: the quote is from the book's cover.) And I bet you'll easily find similar things in Maple, too. William Kahan once said that it is a "folk theorem" that under any implementation of significance arithmetic there must exist chains of operations that lead to either an unwarranted loss of significance, or to a gain in significance, or both. [Quoted from Richard Fateman's famous Mma-Review.] My point is that numerical stability is a game for which there is no win-strategy. Hence, we shouldn't pretend providing one. Cheers -richy. -- Richard B. Kreckel <Richard.Kreckel@Uni-Mainz.DE> <http://wwwthep.physik.uni-mainz.de/~kreckel/>
On Mon, 7 Jan 2002, Richard B. Kreckel wrote:
Oh, if you worry about that, then throwning away digits in the output is not going to help since it is merely cosmetic. If we could have two digits, then compared to 3.3, you are saying that 3.333 is better than 3.3334. Why bother?
Sorry for sticking my oar in, but I think the reason the one inaccurate answer is better than the other is a human one, not a machine or a precision one. Either way you should describe the output behavior concisely, and with the current behavior, the description that should be used is "GiNaC gives you meaningless garbage at the end of your output under <x> circumstances." Not exactly something that sounds like a feature. ;) It's easier (I think) to spot that the number of digits isn't right, but nice to know that at least what you *do* see is correct. Just another county heard from. Cheers, Phil -- "Trying to do something with your life is like sitting down to eat a moose." --Douglas Wood
Hi, On Tue, 8 Jan 2002, Phil Mendelsohn wrote:
Oh, if you worry about that, then throwning away digits in the output is not going to help since it is merely cosmetic. If we could have two digits, then compared to 3.3, you are saying that 3.333 is better than 3.3334. Why bother?
Sorry for sticking my oar in, but I think the reason the one inaccurate answer is better than the other is a human one, not a machine or a precision one. Either way you should describe the output behavior concisely, and with the current behavior, the description that should be used is
"GiNaC gives you meaningless garbage at the end of your output under <x> circumstances."
No, no, no! This is not `meaningless garbage', although the human eye might think it is. This `garbage' is absolutely correct from the computer's point of view. See below...
Not exactly something that sounds like a feature. ;)
It's easier (I think) to spot that the number of digits isn't right, but nice to know that at least what you *do* see is correct.
Okay, we have two options, then: 1) Fix CLN, since this is where that output is really generated. 2) Apply a cosmetical patch to GiNaC that throws away the last numbers. Option 1) We first have to find an error here. Considering the 1/3 issue you may notice that the pure CLN output routine always seems to round up except for those cases where the precision is lower than 17 decimal digits. Why? The reason is that CLN uses whole "digit sequences" consisting of arrays of machine-size words for the mantissa, once a machine-size word cannot handle sign, mantissa and exponent for the accuracy that we have ordered. Such an array has necessarily a power of two binary digits. As such, 1/3 always looks like this in binary notation +.1010101010101010101010101010101010101...e-1 which has to be truncated by the machine to fit a power of two. In other words, the truncation has to take place after a 0 binary digit. But due to rounding, that last 0 must be converted to a 1. What we thus get is this, assuming words of 16-bits: +.10101010101010101010101010101011...e-1 Later, that internal representation is converted to decimal notation. Here are the first few cases, in which I give the mantissa as it is stored internally by CLN in binary notation in the first line, the "exact" output in decimal notation to two more digits and the output in decimal as it is generated by CLN. The last one is the one which seems to offend people when typing `evalf(1/3)' in ginsh. The number of decimal digits corresponds to cln::cl_float_format_t. 17 decimal Digits: 1010101010101010101010101010101010101010101010101010101010101011 0.3333333333333333333423 0.33333333333333333334 27 decimal Digits: 101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101011 0.33333333333333333333333333333543 0.333333333333333333333333333335 37 decimal Digits: 10101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101011 0.33333333333333333333333333333333333333382 0.333333333333333333333333333333333333334 46 decimal Digits: 1010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101011 0.333333333333333333333333333333333333333333333333447 0.3333333333333333333333333333333333333333333333334 56 decimal Digits: 101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101011 0.3333333333333333333333333333333333333333333333333333333333598 0.33333333333333333333333333333333333333333333333333333333336 66 decimal Digits: 10101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101011 0.3333333333333333333333333333333333333333333333333333333333333333333395 0.33333333333333333333333333333333333333333333333333333333333333333334 75 decimal Digits: 1010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101011 0.33333333333333333333333333333333333333333333333333333333333333333333333333333477 0.333333333333333333333333333333333333333333333333333333333333333333333333333335 The rounding is actually working wonderfully, in all these cases. Hence, there does not seem to be anything we must fix in CLN. Option 2) I am not entirely sure but to me this looks purely cosmetic and would only lull the user into a false sense of security. I still fail to see what's so wrong about 0.33333333333333333334. It gives you some valuable information, doesn't it? After all, when you throw away the last digit, we could interpret the number 0.3333333333333333333 as 0.33333333333333333327, but the original is much closer to reality. Best wishes -richy. -- Richard B. Kreckel <Richard.Kreckel@Uni-Mainz.DE> <http://wwwthep.physik.uni-mainz.de/~kreckel/>
On Tue, 8 Jan 2002, Richard B. Kreckel wrote:
No, no, no! This is not `meaningless garbage', although the human eye might think it is. This `garbage' is absolutely correct from the computer's point of view. See below...
You are, of course, correct about the mathematics. My point was more that the human eye is the intended audience of an _output_ routine.
Option 2) I am not entirely sure but to me this looks purely cosmetic and would only lull the user into a false sense of security.
Well, apparently it's another quanitzation error of a different source -- humans can have either a false sense of security or a false sense of panic. ;-)
I still fail to see what's so wrong about 0.33333333333333333334. It gives you some valuable information, doesn't it?
No, actually. I do recall the convention from Numerical Analysis where one more place than is accurate is written, so the fix is probably just to put that as an explicit footnote in the docs, that this is the convention followed. But the last digit is read as "ignore this digit, it really could be anything."
After all, when you throw away the last digit, we could interpret the number 0.3333333333333333333 as 0.33333333333333333327, but the original is much closer to reality.
Again, I think the feeling is that it isn't a matter of interpretation. If you specify 30 digits of accuracy, you want 30 accurate digits. In your example, if you knew you had an answer accurate to <x> digits, you'd be creating the error by interpreting an accurate output. (If I've counted your decimal places correctly.) It all comes down to how many places are *guaranteed*. Sure, it's cosmetic. But beauty is in the eye of the beholder. I promise not to say anything more on the subject. ;-) Cheers, Phil -- "Trying to do something with your life is like sitting down to eat a moose." --Douglas Wood
participants (3)
-
Bob McElrath
-
Phil Mendelsohn
-
Richard B. Kreckel