A quick precision test

[everettr]

(06-03-2014 01:51 AM)Gerson W. Barbosa Wrote:  This is from an old Scientific American article, "How to handle numbers with thousands of digits, and why one might want to", by Fred Gruenberger (Computer Recreations, April 1984). ( "Come e perché trattare numeri con migliaia di cifre", in Le Scienze )

+1 for Fred Gruenberger, and the old Scientific American.

[Gerson W. Barbosa]

(06-03-2014 02:11 PM)everettr Wrote:  

(06-03-2014 01:51 AM)Gerson W. Barbosa Wrote:  This is from an old Scientific American article, "How to handle numbers with thousands of digits, and why one might want to", by Fred Gruenberger (Computer Recreations, April 1984). ( "Come e perché trattare numeri con migliaia di cifre", in Le Scienze )

+1 for Fred Gruenberger, and the old Scientific American.

Our university library had a subscription to Scientific American and lots of issues then. That's one of the few articles I still have a copy of. I hope they still have them. It might be cheaper to go there and scan the articles I am interested in than paying SciAm $7,99 each :-)

Gerson.

[Massimo Gnerucci]

(06-03-2014 09:39 PM)Gerson W. Barbosa Wrote:  

(06-03-2014 02:11 PM)everettr Wrote:  +1 for Fred Gruenberger, and the old Scientific American.

Our university library had a subscription to Scientific American and lots of issues then. That's one of the few articles I still have a copy of. I hope they still have them. It might be cheaper to go there and scan the articles I am interested in than paying SciAm $7,99 each :-)

Gerson.

I've got all "Le Scienze" issues from 1968 to 2008 on 2 DVDs, unfortunately not all sections are present, only main articles.
However I found complete text online.

This way is better, though...

[Paul Dale]

(06-03-2014 12:57 PM)Massimo Gnerucci Wrote:  This one?

That's the one. Shorter than I remember but priceless nonetheless.
I'm probably mixing things up again with the come from.

- Pauli

[Gerson W. Barbosa]

(06-03-2014 11:31 AM)Massimo Gnerucci Wrote:  

(06-03-2014 11:24 AM)walter b Wrote:  Quite natural, isn't it?

d:-)

To us, yes.
Cobol sentence mandatory in every program I wrote...

I must have written no more than half a dozen COBOL programs as a student, but I do remember that sentence (which I liked very much, BTW). Back in the day I would occasionally change the ROM of my MSX computer in order to exchange dots and commas in numeric outputs).

It was a pleasant surprise when I discovered the HP-15C, a replacement to a TI-59, had an equivalent of that built-in (ON/.). I had never bothered about decimal points in my previous calculators, but commas as thousands separators looked quite strange - there ought to have a way to change that. That was the first thing I looked for in the HP-15C manual, I think (no, that was the second, the first one was the whereabouts of a hidden equals key :-)

Cheers,

Gerson.

[Didier Lachieze]

(06-03-2014 10:20 PM)Gerson W. Barbosa Wrote:  It was a pleasant surprise when I discovered the HP-15C, a replacement to a TI-59, had an equivalent of that built-in (ON/.).

You have the exact same feature on the HP 41C with a Paname module

[Gerson W. Barbosa]

(06-03-2014 10:38 PM)Didier Lachieze Wrote:  

(06-03-2014 10:20 PM)Gerson W. Barbosa Wrote:  It was a pleasant surprise when I discovered the HP-15C, a replacement to a TI-59, had an equivalent of that built-in (ON/.).

You have the exact same feature on the HP 41C with a Paname module

I didn't know there was a Paname module for the HP 41C (but I remember the song, in a movie about Edith Piaf I watched recently :-)
CF 28 is fine too! (I didn't have a 41C back then)

Gerson.

[Dieter]

(06-03-2014 10:20 PM)Gerson W. Barbosa Wrote:  Back in the day I would occasionally change the ROM of my MSX computer in order to exchange dots and commas in numeric outputs).

I used a patched Turbo Pascal RTL file to make my preferred compiler do the same. ;-)

Quote:It was a pleasant surprise when I discovered the HP-15C, a replacement to a TI-59, had an equivalent of that built-in (ON/.). I had never bothered about decimal points in my previous calculators, but commas as thousands separators looked quite strange - there ought to have a way to change that.

That's why the 41C offers every option you like: with and without separators (flag 29), with a comma or dot as the decimal marker (flag 28). After a master clear my first keystrokes always are "CF 28" to restore the decimal comma.

BTW I would not care much about the "right" decimal marker if it was not hidden between two digits as it is nowadays. There were times when the dot/comma hat its own digit in the display, e.g. on the 67/97. I think this is much better than what we have on most calculators today. That's why I do not have a problem with numbers displayed in a fixed width font on a computer screen - the decimal marker is perfectly visible, both as a dot or as a comma.

On the other hand I can hardy distinguish the dot/comma setting on a 34s. The marker simply is too small. That's why I usually turn the thousands separators off – in a number with four dots it's hard to tell which of them is the decimal marker.

Dieter

[walter b]

(06-04-2014 01:08 PM)Dieter Wrote:  On the other hand I can hardy distinguish the dot/comma setting on a 34s. The marker simply is too small. That's why I usually turn the thousands separators off – in a number with four dots it's hard to tell which of them is the decimal marker.

That topic is going to be covered in a far friendlier way on the 43S.

d:-)

[Claudio L.]

(06-03-2014 12:36 AM)Paul Dale Wrote:  We always round to 16 or 34 digits when we return the result to the stack. These are the only precisions a user or a keystroke program can access. Unless you are coding numerics in C, this is all you get.

The 39 digits are used inside but only during the evaluation of a built in function. Once the result is known, it is rounded. The reason is simple: guard digits. I don't need to and generally cannot get all 39 digits correct but having the extra digits gives some headroom to get the rounded result correct or at least close.


- Pauli

Paul,

EDIT: Nevermind, I just read the post from G. Barbosa, showing that WP34S matches exactly the results of newRPL. I'm not deleting the information below just for the record, but my post is pretty much useless from here on... except to actually validate that WP34S and newRPL are on the same page.

I think he's onto something. Rounding should not kill so many digits.

I just took the newRPL demo (shameless advertising: download the demo from http://sf.net/p/newrpl ) and did:

Code:


34 SETPREC
1.0000001
<< 1 27 START DUP * NEXT >> EVAL

In human words: Set the precision to 34 digits to match the intermediate rounding you are using in WP34S. Then do the 27 squares by hand and the result is:

674530.4707410845593826891847277722 (if I didn't make a mistake copying by hand the numbers)

The "real" result at 34 digits is:
674530.4707410845593826891780297468

This is obtained from doing:

Code:


1.0000001 2 27 ^ ^

The difference is 6.698E-21, so there's 25 good digits out of the 34. Intermediate rounding only affected the last 9 digits.

So the original poster might have a good point. newRPL does not use decNumber but uses mpdecimal, also correctly rounded and we also use a few digits more but all intermediate calculations are rounded to the user precision so the results from newRPL and WP34S should match exactly, as both should produce 34 correctly rounded digits on each intermediate calculation.

I think his findings might be worth a second look.

Claudio

[walter b]

Claudio,

FYI, you're allowed to delete your post you rated as useless yourself.

[pito]

(06-04-2014 01:25 PM)Claudio L. Wrote:  [quote='Paul Dale' pid='13135' dateline='1401755816']
..
In human words: Set the precision to 34 digits to match the intermediate rounding you are using in WP34S. Then do the 27 squares by hand and the result is:

674530.4707410845593826891847277722 (if I didn't make a mistake copying by hand the numbers)

The "real" result at 34 digits is:
674530.4707410845593826891780297468

This is obtained from doing:

Code:


1.0000001 2 27 ^ ^

The difference is 6.698E-21, so there's 25 good digits out of the 34. Intermediate rounding only affected the last 9 digits.

Not sure whether doing 27 times x^2 or using Xpower2power27 is the same. Power functions are calculated in a different way usually (less precise I bet..).

I've calculated the stuff with my dusty vanilla decnum code (DECNUMDIGITS 34) and get:

Code:

// 27x x^2 TEST
cout<<"27xTEST"<<endl;
x = "1.0000001";
timer = millis;
for (i=1; i<=27; i++) {
  x = x * x;
}
timer = millis - timer;
cout<<"27xTest = "<<x<<endl;
cout<<"Elapsed x1 = "<<timer<<" ms"<<endl;

Code:

27xTEST
27xTest = 674530.4707410845593826891847277722
Elapsed x1 = 6 ms

[Claudio L.]

(06-04-2014 01:33 PM)walter b Wrote:  Claudio,

FYI, you're allowed to delete your post you rated as useless yourself.

I know, but that would remove the shameless advertising!

Claudio

[pito]

And the same with DECNUMDIGITS 60:

Code:

27xTEST
27xTest = 674530.470741084559382689178029746812844444143410342033363075
Elapsed x1 = 11 ms

[Claudio L.]

(06-04-2014 01:38 PM)pito Wrote:  Not sure whether doing 27 times x^2 or using Xpower2power27 is the same. Power functions are calculated in a different way usually (less precise I bet..).

I've calculated the stuff with my dusty vanilla decnum code (DECNUMDIGITS 34) and get:

Code:

// 27x x^2 TEST
cout<<"27xTEST"<<endl;
x = "1.0000001";
timer = millis;
for (i=1; i<=27; i++) {
  x = x * x;
}
timer = millis - timer;
cout<<"27xTest = "<<x<<endl;
cout<<"Elapsed x1 = "<<timer<<" ms"<<endl;

Code:

27xTEST
27xTest = 674530.4707410845593826891847277722
Elapsed x1 = 6 ms

Regarding your suspicion that it's not the same: you are right, it's not.
The power function is guaranteed to give you an accurate result with 34 good digits, for any 2 parameters (in this case, your guess that it's less precise is incorrect).
Doing 2 ^ 27, the result is less than 34 digits, so it's exact (no rounding error).
so 1.0000001 ^ (2^27) should give you an answer with all 34 digits correct, since both numbers are "exact" as given in the problem.

I don't know how decNumber does it in particular, but in general powers with small integer exponents are handled by a sequence of multiplications, while powers of real numbers or large numbers are usually done through logarithms:

Code:

x^a = exp( a * ln(x) )

Which is usually very accurate, as ln() and exp() are both guaranteed to give you 34+GUARD good digits on any number you throw at them.

If you assume that basic operations give you results that are guaranteed up to your user precision, your only source of error is the rounding of those results to your user precision. The less operations you perform the better, so you don't accumulate those rounding errors.

In this case, if you do:

Code:

 1.0000001 LN 2 27 ^ * EXP

You get the same result as doing

Code:

 1.0000001 2 27 ^ ^

but this is by pure chance (with this precision, this particular number behaves like this), since the intermediate rounded results of LN and * could've messed the last digit.

Claudio

[pito]

With 27x x^2 and with POW POW, DECNUMDIGITS 34:

Code:

// 27x POW POW based TEST
cout<<"27x POW BASED TEST"<<endl;
x = "1.0000001";
timer = millis;
x = pow(x, pow("2","27"));
timer = millis - timer;
cout<<"27xTest = "<<x<<endl;
cout<<"Elapsed x1 = "<<timer<<" ms"<<endl;

Code:

27xTEST
27xTest = 674530.4707410845593826891847277722
Elapsed x1 = 6 ms

27x POW BASED TEST
27xTest = 674530.4707410845593826891780297468
Elapsed x1 = 8 ms

And wolfram's alfa online (maybe the most precise) calculated with POWPOW:

Code:

674530.4707410845593826891780297468128444441434103420317423773278390177617568356​46924185036948314117161449392236580649487823578911213711850826688459566614043388​12537833649638778823165610131262467005046619414630533734997699545642384381930207​37082508929177604031150305792570919466251509008479500637190800014648201143130934​98728408041488454981733591279031941696561625036530446399539623519795473743965573​22542637960770829830176516928352599744557589814991829065879180290393434405570265​3738971...

Code:

27x x^2 34dig    674530.4707410845593826891847277722
POWPOW  34dig    674530.4707410845593826891780297468
27x x^2 60dig    674530.470741084559382689178029746812844444143410342033363075
wolfram POWPOW   674530.470741084559382689178029746812844444143410342031742377..

[Paul Dale]

(06-04-2014 01:25 PM)Claudio L. Wrote:  newRPL does not use decNumber but uses mpdecimal, also correctly rounded and we also use a few digits more but all intermediate calculations are rounded to the user precision so the results from newRPL and WP34S should match exactly, as both should produce 34 correctly rounded digits on each intermediate calculation.

By using extra guard digits you get a double rounding step which can result in incorrectly rounded results. A quick example...

If the true result of a computation is 123.456789 but our guard digit format only carries 4 digits, the guarded result will correctly be 123.5. If now, our returned result only has 3 digits, we'll get 124 which is incorrect. It should be 123.

Sure this is contrived but the same can happen regardless of the number of guard digits.

Double rounding is bad. Guard digits are good. Generally guard digits are more good than double rounding is bad, which is why they are used. However, getting a correctly rounded result is, in general, extremely difficult.


- Pauli

[Paul Dale]

(06-04-2014 02:17 PM)Claudio L. Wrote:  The power function is guaranteed to give you an accurate result with 34 good digits, for any 2 parameters (in this case, your guess that it's less precise is incorrect).

I very very much doubt this. Perhaps results with one ULP error but even then I doubt it is that good. The power function is a very difficult one to get right. The GNU C library implementation for 64 bit doubles, in bad cases, falls back to multiple precision calculations carrying over 1600 guard bits. I expect someone rather clever spent a lot of time and effort determining that number as an upper bound to guarantee correctness for the 53 mantissa bits in the result.

Quote:I don't know how decNumber does it in particular, but in general powers with small integer exponents are handled by a sequence of multiplications, while powers of real numbers or large numbers are usually done through logarithms:

Code:

x^a = exp( a * ln(x) )

Which is usually very accurate, as ln() and exp() are both guaranteed to give you 34+GUARD good digits on any number you throw at them.

Two points here:

1. Expressing power as x^a = exp( a * ln(x) ) is prone to digit loss. Don't expect an accurate result unless you have a lot of extra digits. See the 1600 above
2. The decNumber versions I've looked at do not guarantee ln() and exp() are correctly rounded, the source just claims close but usually within one ULP.
   

A good rule of thumb to follow is to assume correct rounding for the four arithmetic operations and square root -- these are mandated by the IEEE floating point standards these days. Everything else is close.

- Pauli

[Claudio L.]

(06-04-2014 10:09 PM)Paul Dale Wrote:  

(06-04-2014 02:17 PM)Claudio L. Wrote:  The power function is guaranteed to give you an accurate result with 34 good digits, for any 2 parameters (in this case, your guess that it's less precise is incorrect).

I very very much doubt this. Perhaps results with one ULP error but even then I doubt it is that good. The power function is a very difficult one to get right. The GNU C library implementation for 64 bit doubles, in bad cases, falls back to multiple precision calculations carrying over 1600 guard bits. I expect someone rather clever spent a lot of time and effort determining that number as an upper bound to guarantee correctness for the 53 mantissa bits in the result.

I don't know that much, but the documentation of mpdecimal says the error is 1 ULP, so we use 9 extra digits and then we round, so the 'bad rounding' case would have to be in the form 999999995 in order for the 1 ULP error to give a badly rounded result. Notice that I said 'accurate digits' not 'correctly rounded'. The 1 ULP case may still happen, with a chance of 1 in 10^9.

(06-04-2014 10:09 PM)Paul Dale Wrote:  

Quote:I don't know how decNumber does it in particular, but in general powers with small integer exponents are handled by a sequence of multiplications, while powers of real numbers or large numbers are usually done through logarithms:

Code:

x^a = exp( a * ln(x) )

Which is usually very accurate, as ln() and exp() are both guaranteed to give you 34+GUARD good digits on any number you throw at them.

Two points here:

1. Expressing power as x^a = exp( a * ln(x) ) is prone to digit loss. Don't expect an accurate result unless you have a lot of extra digits. See the 1600 above
2. The decNumber versions I've looked at do not guarantee ln() and exp() are correctly rounded, the source just claims close but usually within one ULP.
   

A good rule of thumb to follow is to assume correct rounding for the four arithmetic operations and square root -- these are mandated by the IEEE floating point standards these days. Everything else is close.

- Pauli

I agree, but "close" as in the last digit may be off by 1 due to rounding and only in some corner cases (like 1 in 10^9 chance), on most numbers you should get the right value.

I'd love to test some corner cases to see how the algorithms behave (like that case that needs 1600 digits, do you have any links to good material on this subject?). I did some tests of my own but I don't know what the corner cases are. I'd like to have some literature where they list cases that I can use to stress the algorithms and see the results.

Claudio

[Paul Dale]

(06-05-2014 01:18 AM)Claudio L. Wrote:  I'd love to test some corner cases to see how the algorithms behave (like that case that needs 1600 digits, do you have any links to good material on this subject?). I did some tests of my own but I don't know what the corner cases are. I'd like to have some literature where they list cases that I can use to stress the algorithms and see the results.

I doubt that the 1600 bits (not digits) is actually required in reality. I suspect rather, that a proof was devised that this was sufficient for all possible input values rather than finding the absolute worst case. There has been a fair bit of work done on actually locating the worst cases in the past few years.

A quick check on google leads to two links:

For binary: http://hal.inria.fr/inria-00072594/PS/RR-4044.ps
For decimal: http://perso.ens-lyon.fr/damien.stehle/d...malexp.pdf

There are plenty of others a search away.

The example at the bottom of the second page of the second paper is interesting. Forty seven correct decimal digits being required to correctly round towards zero -- 34 plus 9 guards won't get this right and this is only a 16 digit input, a 34 digit one would be worse.

As I've said many times, "numerical analysis is hard".

- Pauli