Threaded Mode | Linear Mode

Claudio L. · 10-07-2015, 03:18 PM

Reading the HHC 2015 paper by Wlodek about the Savage benchmark on the Prime, it occurred to me to run it in newRPL to see how it performs.
Despite a disappointing 15 seconds mark (this benchmarks hits newRPL very hard, due to the variable precision implementation of transcendental functions being significantly slower than fixed-precision counterparts), one thing that attracted my attention was that the results were not identical.
I expected newRPL to return the same exact value when set to 12 digits precision.
Old HP's consistently returned 2499.99948647, while newRPL returned 2499.99942402.

Intrigued by this, and assuming there had to be a bug in newRPL somewhere, I modified the benchmark to leave all 2500 partial results in the stack.
Turns out that results were identical all the way to iteration number 1470.
The value in the stack was exactly 1469.99992223, and the difference came from ATAN, which returns 1.57011605476 on the 50g, but 1.57011605475 in newRPL, and values diverged from there.
The actual result for ATAN(1469.99992223) is:
1.5701160547549999736...

So it's clearly a case of double-rounding. I think the internal 15-digit precision would round 49999 to 500, and therefore round the last digit to 6 instead of 5.

Now I forced the loop to start on the 50g from the corrected value, and this happened again for number 1647.99989715, also in the ATAN function.
In this case ATAN gives 1.57018953071499665..., which again causes the double rounding and values diverge again.

The surprise was... it diverged to the exact same value as before!
Result was 2499.99948647, even when starting from the "correct" sequence after value TAN(ATAN(1469.99992223))+1 = 1470.99991143 was used as starting point.

Why is this surprising? because the difference after doing ATAN and TAN on 1469.99992223 was 2160e-8, and for example the next error at 1647.99989715 is 2716e-8. I don't see how the sequence could "self-correct" to give the exact same result? Perhaps somebody has an explanation?

Dieter · (This post was last modified: 10-07-2015 07:38 PM by Dieter.)

(10-07-2015 03:18 PM)Claudio L. Wrote: The value in the stack was exactly 1469.99992223, and the difference came from ATAN, which returns 1.57011605476 on the 50g, but 1.57011605475 in newRPL, and values diverged from there.
The actual result for ATAN(1469.99992223) is:
1.5701160547549999736...

So it's clearly a case of double-rounding. I think the internal 15-digit precision would round 49999 to 500, and therefore round the last digit to 6 instead of 5.

Now I forced the loop to start on the 50g from the corrected value, and this happened again for number 1647.99989715, also in the ATAN function.
In this case ATAN gives 1.57018953071499665..., which again causes the double rounding and values diverge again.

I do not think this is caused by double rounding. The internal extended precision routines usually truncate the result after 13 resp. 15 digits. This would yield the correct result here: 1.57011605475499 => 1.57011605475.

But all this will only work as expected if (!) the internal result is exact in all 13 resp. 15 digits. Well... it isn't. If it was, there was not need for guard digits. ;-)

The internal routines for the transcendental functions actually are not that exact, that's why the final 12-digit result can be off by one ULP (or even more in some special cases). I cannot say much about the internal 15-digit routines of the 50G, but for the 41C the MCode console shows slight errors in the last (13th) digit even for simple cases like e¹ = 2,718281828458.

That's why in close cases (e.g. the one you mentioned) the last digit may be off by one. I don't think it's a rounding issue, it's simply a slight error in the 15-digit result.

Dieter

Claudio L. · 10-08-2015, 12:34 PM

(10-07-2015 07:37 PM)Dieter Wrote: That's why in close cases (e.g. the one you mentioned) the last digit may be off by one. I don't think it's a rounding issue, it's simply a slight error in the 15-digit result.

I think you are right, which makes it even more incredible that some random error in the last digit could "self-correct" back to the original value. I refuse to believe it's a coincidence, there has to be a mathematical explanation for that... but it's beyond my knowledge.

Dieter · (This post was last modified: 10-08-2015 07:13 PM by Dieter.)

(10-08-2015 12:34 PM)Claudio L. Wrote: I think you are right, which makes it even more incredible that some random error in the last digit could "self-correct" back to the original value. I refuse to believe it's a coincidence, there has to be a mathematical explanation for that... but it's beyond my knowledge.

Maybe we can solve this... if (!) I only understood what makes you wonder here.

From your initial post:

(10-07-2015 03:18 PM)Claudio L. Wrote: The surprise was... it diverged to the exact same value as before!
Result was 2499.99948647, even when starting from the "correct" sequence after value TAN(ATAN(1469.99992223))+1 = 1470.99991143 was used as starting point.

Why is this surprising? because the difference after doing ATAN and TAN on 1469.99992223 was 2160e-8, and for example the next error at 1647.99989715 is 2716e-8. I don't see how the sequence could "self-correct" to give the exact same result? Perhaps somebody has an explanation?

The difference after an tan(arctan x) for x=1,57... here is 1 ULP of x = 1E–11 divided by the first arctan derivative, i.e. 1/(x²+1). Which yields ~ 1E–11*(1470²+1) = 2,161E–5 resp. 1E–11*(1648²+1) = 2,716E–5. Which is what you get. This also shows that a very, very slight error in x (merely one single ULP due to an even lesser error in the 15-digit result) will cause large errors in the next step. And these errors are very similar for x=1469, 1470 or 1471. So if one result is 1 ULP high and one of the next is 1 ULP low this yields very similar errors that may compensate each other.

But maybe I did not understand your surprise correctly. In this case please provide a step-by-step analysis.

Dieter

Claudio L. · 10-08-2015, 08:31 PM

(10-08-2015 07:11 PM)Dieter Wrote: The difference after an tan(arctan x) for x=1,57... here is 1 ULP of x = 1E–11 divided by the first arctan derivative, i.e. 1/(x²+1). Which yields ~ 1E–11*(1470²+1) = 2,161E–5 resp. 1E–11*(1648²+1) = 2,716E–5. Which is what you get. This also shows that a very, very slight error in x (merely one single ULP due to an even lesser error in the 15-digit result) will cause large errors in the next step. And these errors are very similar for x=1469, 1470 or 1471. So if one result is 1 ULP high and one of the next is 1 ULP low this yields very similar errors that may compensate each other.

You are getting close, and I think you understood perfectly well what I mean.
So for 1469 there's +2161e-5, then no more errors until 1648, where you have other +2716e-5, then who knows where there are more errors like this, from there all the way up to 2499.
From 1 to 1469 there are no errors, then from there on there's only errors in a few numbers, their magnitudes are not the same, as you correctly estimated 1/(x²+1), so how in the world the value goes back to the same result?
[/quote]

(10-08-2015 07:11 PM)Dieter Wrote: But maybe I did not understand your surprise correctly. In this case please provide a step-by-step analysis.

I think you did, but here's the step-by-step:

The standard Savage Benchmark:
<< 1. 1. 2499. START SQ √ LN EXP ATAN TAN 1 + NEXT >>

Doing SQ √ in that order causes all results to be exact, so it has no effect for this analysis (it does for the speed benchmark of course).
The simplified loop is:
<< 1. 1. 2499. START LN EXP ATAN TAN 1 + DUP NEXT >>

in which I added a DUP at the end to keep all results in the stack.

The last value is 2499.99948647 (remember this one).

Now if you do:
1025 DROPN you'll see on your stack where the first problem happens.
1469.99992223, doing the loop by hand:
LN EXP gives you an exact result, back 1469.99992223
ATAN TAN 1 + gives the next value (incorrect): 1470.99993303

So, let's compute the correct value:
1469.99992223 ATAN 1E-11 - TAN 1 +
which gives 1470.99991143.

Let's now run the loop, but starting from the correct value and going up to 2499:

<< 1470.99991143 1471. 2499. START LN EXP ATAN TAN 1 + DUP NEXT >>

The only change is the initial value, and we start from iteration number 1471 going up to 2499.

Now look at the final value... 2499.99948647 (why the same?)

Here's another take: side by side computation of both sequences
<< { 1470.99991143 1470.99993303 } 1471. 2499. START LN EXP ATAN TAN 1 ADD DUP NEXT >>

The only change is that the initial value is now a list, and + was changed to ADD (on the 50g, those trying this on newRPL must keep the +).

I discovered that the lists get "in sync" again exactly at 1647, where the next error occurs.

The reduced loop:
<< { 1470.99991143 1470.99993303 } 1471. 1650. START LN EXP ATAN TAN 1 ADD DUP NEXT >>

shows exactly this: the first few values are identical, and then they are different.

Question still remains: why does it self-correct?

Dieter · 10-09-2015, 07:45 AM

(10-08-2015 08:31 PM)Claudio L. Wrote: Question still remains: why does it self-correct?

I do not think it "self-corrects", it just happends to round once up, once down. Depending on where you stop the loop the results sometimes match, sometimes they don't.

Take a look at tan(arctan/1470.99991143)). The exact arctan is 1,570 1165 1720 5390... Both the 35s and the 39gII, two 12-digit calculators, correctly round this to 1,57011651721. But the former returns the tangent of this (1470,9999 2140 4527...) as 1470,99992141 while the latter displays 1470,99992140. So there is always a chance that the last digit is off by one. I remember the accuracy discussion in the HP15C Advanced Functions Handbook. Here the trig functions are said to return results with an error between 0,5 and 1 ULP. Which means the last rounded digit can be one digit high or low.

If you google "savage benchmark" you will come across a PDF on technicalc.org that shows results for several calculators. It also includes a graph that shows the error after 1...5000 iterations. The essential point here is that the error graph does not continuously drop, as the accuracy of the result after n loops is expected to decrease, but – while generally decreasing – it bounces up and down. So at some points the accuracy actually increases (which is what you may call "error compensation") and at other points it decreases again. Which is what can be expected due to the more or less random errors in the last place that may occur any time. The final error after n loops depends on where you stop – it may be smaller or larger, and even smaller for a larger number of loops. ;-)

Dieter

Claudio L. · 10-09-2015, 12:53 PM

(10-09-2015 07:45 AM)Dieter Wrote:
(10-08-2015 08:31 PM)Claudio L. Wrote: Question still remains: why does it self-correct?

I do not think it "self-corrects", it just happends to round once up, once down. Depending on where you stop the loop the results sometimes match, sometimes they don't.

There's a lot more than that going on!
Take the first iteration, for example:
1470.99993303
LN EXP --> 1470.99993302 (there's -1 change introduced here)
ATAN TAN 1 + --> 1471.99994304

1470.99991143
LN EXP --> 1470.99991143 (but not on this path...)
ATAN TAN 1 + --> 1471.99992140

So deviations from the "exact" result happen somewhat randomly on both paths, sometimes it +/- 1 in the LN/EXP combo, sometimes it's +/-1 in the ATAN function, which ends up around +/-2e-4 in the final result, amplified by TAN.
By the time you get to iteration 1647, you've been following 2 different paths for 170 iterations, and you end up with 2 different numbers, both containing a number of random deviations in the last 4 digits. Magically, this iteration produces the exact same result, on both paths.

The last position where the 2 paths differ has:
1647.99989715 and 1647.99992427

Applying the exact same operators to these 2 different numbers produces the exact same result...
1647.99989715 LN EXP --> 1647.99989715
ATAN --> 1.57018953072

1647.99992427 LN EXP --> 1647.99992428
ATAN --> 1.57018953072

Turns out that:
tan(1.57018953071)=1647.99988358
tan(1.57018953072)=1647.99991074
tan(1.57018953073)=1647.99993790

This is a lot less magical than it appeared at the beginning (at least to me). All numbers between 1074 and 3790 in the last 4 digits will produce the same result. The difference 1e-11 in the tangent produces +/- 2716 in the last 4 digits. If we think the last 4 digits are completely random, then there is a chance 2716 in 10000, or roughly 1 in 4, that they will produce the same result. That's a very good chance that the paths will get in sync again.
I think the mystery is solved, the loop runs until these somewhat random errors produce values that are within the +/- 2716 in the last four digits, then the ATAN/TAN pair will coalesce these range of numbers into the same result, syncing back the paths.

Dieter · (This post was last modified: 10-09-2015 08:09 PM by Dieter.)

(10-09-2015 12:53 PM)Claudio L. Wrote: This is a lot less magical than it appeared at the beginning (at least to me). All numbers between 1074 and 3790 in the last 4 digits will produce the same result.
(...)
I think the mystery is solved...

So finally there is an answer. ;-)

But, maybe more important for you:

(10-07-2015 03:18 PM)Claudio L. Wrote: Old HP's consistently returned 2499.99948647, while newRPL returned 2499.99942402.

This is great news for newRPL users. Because it is the perfect result for a 12-digit machine. I did the Savage benchmark in a 34-digit environment (WP34s) with n-digit rounding after every single operation. For n=12 this is exactly the "perfect" result that should appear if every single intermediate result was correctly rounded to 12 digits.

So...

(10-09-2015 12:53 PM)Claudio L. Wrote: Intrigued by this, and assuming there had to be a bug in newRPL somewhere ...

...the bug is seems to be not in newRPL, but in the "classic" (Saturn based) implementation.

For comparison, here are the results for some other cases:

Code:

n=8     2508,9828

n=10    2499,970323

n=12    2499,99942402

n=15    2499,99999975044

n=16    2500,000000006490

This implies the 34s runs in rounding mode RM 0 (round to closest value, exactly 0,500000... rounds to next even digit). Classic HPs rounded 0.5 always up, which is RM 1 on the 34s. However, the results will only differ if there is an intermediate result where the n+1st digit is 5 and all following ones are zero. ;-)

Dieter

Rodger Rosenbaum · 11-03-2015, 05:19 PM

There was considerable discussion of all this here: http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=103356

Also discussion on a similar topic here: http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=103151

Claudio L. · 11-05-2015, 02:18 AM

(11-03-2015 05:19 PM)Rodger Rosenbaum Wrote: There was considerable discussion of all this here: http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=103356

Also discussion on a similar topic here: http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=103151

Thanks, my conclusion was the same as in that thread... only it took me 9 years more to figure it out!