Accuracy and the power function

[Dieter]

Recently I posted a HP35s program that determines the Bernoulli numbers \(B_n\). The formula includes powers of \(2\large\pi\) with exponents up to 372. Since \(\large\pi\) carries just 12 digits, the results with large exponents are somewhat inaccurate even though the 12-digit-value agrees with the true value rounded to 13 digits, so that the relative error in \(\large\pi\) is less than 6,6 E-14.

While this deviation can be compensated quite easily, there is another error that can reach magnitudes I did not expect. In the mentioned program it here and there may cause error peaks while most other results are within 2 ULP.

Let's assume \(2\large\pi\) = 6,28318530718, i.e. the best possible 12-digit value. Now evaluate some powers of this:

\(6,28318530718^{-76}  =  2,179 3651 6466 35... · 10^{-61}\)
\(6,28318530718^{ 200}  =  4,324 8761 1407 46... · 10^{ 159}\)
\(6,28318530718^{ 372}  =  8,373 5772 1132 00... · 10^{ 296}\)

Please try this with your 12-digit HP and see what you get.

Dieter

[Joe Horn]

(02-16-2014 02:43 PM)Dieter Wrote:  Let's assume \(2\large\pi\) = 6,28318530718, i.e. the best possible 12-digit value. Now evaluate some powers of this:

\(6,28318530718^{-76}  =  2,179 3651 6466 35... · 10^{-61}\)
\(6,28318530718^{ 200}  =  4,324 8761 1401 74... · 10^{ 159}\)
\(6,28318530718^{ 372}  =  8,373 5772 1132 00... · 10^{ 296}\)

The middle example contains an extraneous digit near the end. It should read:
\(6,28318530718^{ 200} = 4,324\ 8761\ 1407\ 46... · 10^{ 159}\)

Quote:Please try this with your 12-digit HP and see what you get.

On the HP 50g, the errors on these calculations are 1 ULP, -1 ULP, and -4 ULP, respectively.

[Dieter]

(02-16-2014 03:11 PM)Joe Horn Wrote:  The middle example contains an extraneous digit near the end.

Right. Has been corrected now.

Quote:On the HP 50g, the errors on these calculations are 1 ULP, -1 ULP, and -4 ULP, respectively.

Same here on a 35s and on the 39G emulator. Hm, 4 ULP in a simple standard function. I admit I would not have expected an error that large. Not in an HP.

Dieter

[Joe Horn]

Methinks it may be explained by the HP-15C Advanced Functions Handbook (warning: this searchable PDF is 75.8 MB), which devotes 40 pages to "Accuracy of Numerical Calculations". Keep in mind that the 15C had only a 10-digit mantissa, with intermediate internal math using 13-digit mantissas. Page 179 includes the following:

Quote:A function that grows to ∞ or decays to 0 exponentially fast as its argument approaches ±∞ may suffer an error larger than one unit in its 10th significant digit, but only if its magnitude is smaller than 10^-20 or larger than 10²⁰; and though the relative error gets worse as the result gets more extreme (small or large), the error stays below three units in the last (10th) significant digit. The reason for this error is explained later. Functions so affected are e^x, y^x, x! (for noninteger x), SINH, and COSH for real arguments. The worst case known is 3²⁰¹, which is calculated as 7.968419664 × 10⁹⁵. The last digit 4 should be 6 instead ... [jump to page 183] To understand the error in 3²⁰¹, note that this is calculated as e^{201 ln(3)} = e^220.821... To keep the final relative error below one unit in the 10th significant digit, 201 ln(3) would have to be calculated with an absolute error rather smaller than 10^-10, which would entail carrying at least 14 significant digits for that intermediate value. The calculator does carry 13 significant digits for certain intermediate calculations of its own, but a 14th digit would cost more than it's worth. [bold emphasis added by -jkh-]

Since 13 internal digits can result in errors up to 3 ULP, I hypothesize that 15 internal digits (as most modern HP's use) can yield errors up to 4 ULP in those functions, if their inputs have sufficiently "extreme" magnitudes. I think your 3rd example above validates this hypothesis.

[Dieter]

(02-16-2014 10:11 PM)Joe Horn Wrote:  Methinks it may be explained by the HP-15C Advanced Functions Handbook (warning: this searchable PDF is 75.8 MB), which devotes 40 pages to "Accuracy of Numerical Calculations".

Ah, the 15C Advanced Functions Handbook. I have a PDF copy here and of course I read the section on accuracy more than once. By the way, the current PDF at hp.com (2012 version) has merely 3,8 MB.

Right, the power function \(a^b\) is evaluated as \(e^{b · ln a}\). But both the 10-digit calculators as well as those with 12 digits use three additional guard digits, so the expected error should be the same. The more important point here seems to be the much wider working range of the 12 digit devices:

For a result less than 1E+100 the largest possible exponent is 230,258... If the 13-digit internal result is off by 1 ULP, the absolute error is 1E-10. Which after the following exponentiation results in a relative error of 1E-10 or at most 1 ULP in the returned result.

If results up to 1E+500 are possible, the largest exponent is 1151,29..., i.e. there is one more digit left of the decimal point and therefore one digit less to the right. An error of merely one ULP in the 15th digit of the exponent equals 1E-11, so the final result may have a relative error of up to 10 ULP (!). Even if the exponent is correctly rounded to 15 digits, the result may be off by 5 ULP.

However, in the example with 4 ULP error the internally calculated exponent is 372*ln(2*pi) = 683,69... One ULP more or less equals an absolute error of 1E-12 or an (acceptable) error of 1 ULP in the returned result. If larger errors occur, one may conclude that 15 digit internal precision does not mean that all 15 are exact.

Dieter

[Joe Horn]

(02-18-2014 08:28 PM)Dieter Wrote:  By the way, the current PDF at hp.com (2012 version) has merely 3,8 MB.

Thanks! HP's copy is much cleaner than my scan, too. HP's version is currently downloadable from HERE.

[Werner]

(02-18-2014 08:28 PM)Dieter Wrote:  However, in the example with 4 ULP error the internally calculated exponent is 372*ln(2*pi) = 683,69... One ULP more or less equals an absolute error of 1E-12 or an (acceptable) error of 1 ULP in the returned result. If larger errors occur, one may conclude that 15 digit internal precision does not mean that all 15 are exact.

There are two things at work here.
First, the 15-digit routines indeed are not exact to the last digit. In your example, %%LN(6.28318530718) is off by one ULP.
Secondly, the 15-digit routines (%%*, %%+, %%-, %%/) do not round, they truncate. They are, after all, meant as intermediates to produce a correctly rounded 12-digit result (for the elementary operations).

Code:


          %%LN                    %% 372 %%*              %%EXP
exact   1.8378770664094112978.. 683.69026870430100..    8.37357721132000..e296
HP48    1.83787706640940        683.690268704296        8.37357721127809e296
trunc15 1.83787706640941        683.690268704300        8.37357721131159e296
round15                         683.690268704301        8.37357721131998e296

Not only is the %%LN one ULP off, but the subsequent multiplication by 372 actually ends in ..968, but is truncated to ..96. The last line demonstrates what a 15-digit %%LN, 15-digit %%EXP and 15-digit rounding arithmetic would produce.

Werner

[walter b]

(02-19-2014 09:54 AM)Werner Wrote:  In your example, %%LN(6.28318530718) is off by one ULP.
Secondly, the 15-digit routines (%%*, %%+, %%-, %%/) do not round, they truncate. They are, after all, meant as intermediates to produce a correctly rounded 12-digit result (for the elementary operations).

Code:


          %%LN                    %% 372 %%*              %%EXP
exact   1.8378770664094112978.. 683.69026870430100..    8.37357721132000..e296
HP48    1.83787706640940        683.690268704296        8.37357721127809e296
trunc15 1.83787706640941        683.690268704300        8.37357721131159e296
round15                         683.690268704301        8.37357721131998e296

Not only is the %%LN one ULP off, but the subsequent multiplication by 372 actually ends in ..968, but is truncated to ..96. The last line demonstrates what a 15-digit %%LN, 15-digit %%EXP and 15-digit rounding arithmetic would produce.

Anything special you want to tell us using %% ??

d#-?

[Werner]

For the uninitiated, the %%-functions are the 48's internal 15-digit versions of the User functions: LN, EXP, *, + , ...
Sorry!

Werner

[Dieter]

(02-19-2014 09:54 AM)Werner Wrote:  There are two things at work here.
First, the 15-digit routines indeed are not exact to the last digit. In your example, %%LN(6.28318530718) is off by one ULP.
Secondly, the 15-digit routines (%%*, %%+, %%-, %%/) do not round, they truncate.

Ah, thank you very much, that's exactly what I have been looking for. I tried to simulate the error on a 34s in 16-digit mode with a similar calculation. So the key here is the logarithm that's off by one ULP.

Truncation of the last (guard) digit(s) seems to be a common practice. IIRC the TI58/59 did so as well, and today you still can see truncation at work, for instance on the 35s. Try the tangent function with very small angles:

Code:

[MODE] DEG
0,00001   [tan]  →  1,74532920000 E-7
0,0000001 [tan]  →  1,74532000000 E-9

[MODE] RAD
0,0005  [tan]  →  5,00000041660 E-4
0,0001  [tan]  →  1,00000000330 E-4
0,00001 [tan]  →  1,00000000000 E-5

In degrees mode, arguments between 1E-3 and 1E-9 have approx. 11...6 valid digits, where the last is truncated or even 1 unit low.

Dieter

[Werner]

Truncation (and keeping track of lost non-zero digits with the Sticky Bit) in the 15-digit routines is not only common practice, it's a necessity if you want to be able to produce a correctly rounded 12-digit result for simple calculations.
Suppose the 15-digit routines rounded as well, then:
5.02000000121e15 + 4995 (exact result: 5.020000001214995e15) would be rounded to 15 digits as 5.02000000121500e15
5.02000000121500e15 rounded to even to 12 digits is 5.02000000122e15
yet the 16-digit number rounded to 12 digits is 5.02000000121e15

Werner

[Dieter]

(02-19-2014 01:07 PM)Werner Wrote:  Truncation (and keeping track of lost non-zero digits with the Sticky Bit) in the 15-digit routines is not only common practice, it's a necessity if you want to be able to produce a correctly rounded 12-digit result for simple calculations.

Thank you for this example.
Here's another one:

5.02000000124e15 + 5009 (exact result: 5.020000001245009e15) would be rounded to 15 digits as 5.02000000124501e15.
After rounding this to 12 digits (towards odd or even) the result is 5.02000000125e15.
Which is the correct 12-digit value.

5.02000000124e15 + 5009 (exact result: 5.020000001245009e15) would be truncated after 15 digits to 5.02000000124500e15.
After rounding this to 12 digits (towards even) the result is 5.02000000124e15.
Which is 1 ULP low.

So, as usual, it depends. ;-)

Dieter

[Werner]

Hi Dieter,
no it doesn't. That's where the Sticky Bit comes in: the CPU works with 15 digits only, aligning the mantissa's (in case of addition), or repeatedly add-and-shift in case of multiplication. In both cases it will keep track of the digits lost to the right, just to know whether they're all zero or not. That's why it's called a 'Sticky Bit': once you shifted a non-zero digit to the right, it is set, and remains set.
So, in your example:
5.02000000124e15 + 5009 (exact result: 5.020000001245009e15) would be truncated to 5.02000000124500, BUT the Sticky Bit would be set, as the '9' we lost is non-zero.
And the routine that rounds 15 digits to 12 will round a trailing '500' up when it is set, and round to even when it is clear.
That's why all 12-digit HPs do get the correctly rounded answer for this example (or any other), yet they all truncate their 15-digit result.

Werner

[Dieter]

(02-21-2014 07:11 AM)Werner Wrote:  That's where the Sticky Bit comes in

I see. So this sticky bit avoids some errors caused by simple truncation.

Since you seem to know the internals quite well: do you know of some kind of reference or summary regarding the accuracy of the internal 15-digit functions? I recently learned that the natural log may be more than 1 ULP off, so I wonder how accurate (or not) the other functions may be.

I would also be interested very much in information on this matter regarding 10-digit devices with 13-digit internal precision, for instance the 41-series.

Dieter

[Thomas Klemm]

(02-23-2014 02:19 PM)Dieter Wrote:  I would also be interested very much in information on this matter regarding 10-digit devices with 13-digit internal precision, for instance the 41-series.

Not for the HP-41C but for the HP-15C I analyzed the Debug Trace in Nonpareil.
As far as I understand there is no "Sticky Bit". Excessive digits are lost.
The result (13 digits) is then rounded half away from zero to 10 digits.
The trick used is to double these last 3 digits and check for a carry.

• 499+499=998 => no carry
• 500+500=1000 => carry
  

In case of a carry 1 is added to the result. Which could set the carry again.
But this is handled separately and I have no trace added for this scenario.

Though I didn't check the sources of the HP-48 I assume the same algorithm is used.

Kind regards
Thomas

5.020000121e13 + 4999:
Shift right register c. The 9 at the end of 4999 is lost:

Code:

 a=04999000000513 b=00000000000513 c=00000000004999
 10364: 1716 c sr    w                             
cycle 150136  P=3 q=4 carry=0  stat=.............. 
 a=04999000000513 b=00000000000513 c=00000000000499

Add both numbers in a and b now that they are alligned:

Code:

 a=00000000000499 b=05020000121000 c=00000000000513
 10376: 0456 a=a+b   w                             
cycle 150147  P=2 q=4 carry=0  stat=.............. 
 a=05020000121499 b=05020000121000 c=00000000000513

Round half away from zero:

Code:

 a=05020000121499 b=00000000000013 c=05020000121499
 07406: 0746 c=c+c   x                             
cycle 150170  P=c q=4 carry=0  stat=.............. 
 a=05020000121499 b=00000000000013 c=05020000121998
 07407: 0153 ?nc goto 07424                        
cycle 150171  P=c q=4 carry=0  stat=.............. 
 a=05020000121499 b=00000000000013 c=05020000121998
 07424: 0306 c=b     x                             
cycle 150172  P=c q=4 carry=0  stat=.............. 
 a=05020000121499 b=00000000000013 c=05020000121013

5.020000121e13 + 5009:
Shift right register c. The 9 at the end of 5009 is lost:

Code:

 a=05009000000513 b=00000000000513 c=00000000005009
 10364: 1716 c sr    w                             
cycle 154224  P=3 q=4 carry=0  stat=.............. 
 a=05009000000513 b=00000000000513 c=00000000000500

Add both numbers in a and b now that they are alligned:

Code:

 a=00000000000500 b=05020000121000 c=00000000000513
 10376: 0456 a=a+b   w                             
cycle 154235  P=2 q=4 carry=0  stat=.............. 
 a=05020000121500 b=05020000121000 c=00000000000513

Round half away from zero:

Code:

 a=05020000121500 b=00000000000013 c=05020000121500
 07406: 0746 c=c+c   x                             
cycle 154258  P=c q=4 carry=1  stat=.............. 
 a=05020000121500 b=00000000000013 c=05020000121000
 07407: 0153 ?nc goto 07424                        
cycle 154259  P=c q=4 carry=0  stat=.............. 
 a=05020000121500 b=00000000000013 c=05020000121000
 07410: 1072 c=c+1   m                             
cycle 154260  P=c q=4 carry=0  stat=.............. 
 a=05020000121500 b=00000000000013 c=05020000122000
 07411: 0133 ?nc goto 07424                        
cycle 154261  P=c q=4 carry=0  stat=.............. 
 a=05020000121500 b=00000000000013 c=05020000122000
 07424: 0306 c=b     x                             
cycle 154262  P=c q=4 carry=0  stat=.............. 
 a=05020000121500 b=00000000000013 c=05020000122013

5.020000124e13 + 4999:
Shift right register c. The 9 at the end of 4999 is lost:

Code:

 a=04999000000513 b=00000000000513 c=00000000004999
 10364: 1716 c sr    w                             
cycle 169908  P=3 q=4 carry=0  stat=.............. 
 a=04999000000513 b=00000000000513 c=00000000000499

Add both numbers in a and b now that they are alligned:

Code:

 a=00000000000499 b=05020000124000 c=00000000000513
 10376: 0456 a=a+b   w                             
cycle 169919  P=2 q=4 carry=0  stat=.............. 
 a=05020000124499 b=05020000124000 c=00000000000513

Round half away from zero:

Code:

 a=05020000124499 b=00000000000013 c=05020000124499
 07406: 0746 c=c+c   x                             
cycle 169942  P=c q=4 carry=0  stat=.............. 
 a=05020000124499 b=00000000000013 c=05020000124998
 07407: 0153 ?nc goto 07424                        
cycle 169943  P=c q=4 carry=0  stat=.............. 
 a=05020000124499 b=00000000000013 c=05020000124998
 07424: 0306 c=b     x                             
cycle 169944  P=c q=4 carry=0  stat=.............. 
 a=05020000124499 b=00000000000013 c=05020000124013

5.020000124e13 + 5009:
Shift right register c. The 9 at the end of 5009 is lost:

Code:

 a=05009000000513 b=00000000000513 c=00000000005009
 10364: 1716 c sr    w
cycle 173995  P=3 q=4 carry=0  stat=..............
 a=05009000000513 b=00000000000513 c=00000000000500

Add both numbers in a and b now that they are alligned:

Code:

 a=00000000000500 b=05020000124000 c=00000000000513
 10376: 0456 a=a+b   w                             
cycle 174006  P=2 q=4 carry=0  stat=.............. 
 a=05020000124500 b=05020000124000 c=00000000000513

Round half away from zero:

Code:

 a=05020000124500 b=00000000000013 c=05020000124500
 07406: 0746 c=c+c   x                             
cycle 174029  P=c q=4 carry=1  stat=.............. 
 a=05020000124500 b=00000000000013 c=05020000124000
 07407: 0153 ?nc goto 07424                        
cycle 174030  P=c q=4 carry=0  stat=.............. 
 a=05020000124500 b=00000000000013 c=05020000124000
 07410: 1072 c=c+1   m                             
cycle 174031  P=c q=4 carry=0  stat=.............. 
 a=05020000124500 b=00000000000013 c=05020000125000
 07411: 0133 ?nc goto 07424                        
cycle 174032  P=c q=4 carry=0  stat=.............. 
 a=05020000124500 b=00000000000013 c=05020000125000
 07424: 0306 c=b     x                             
cycle 174033  P=c q=4 carry=0  stat=.............. 
 a=05020000124500 b=00000000000013 c=05020000125013

[Werner]

(02-23-2014 02:19 PM)Dieter Wrote:  I see. So this sticky bit avoids some errors caused by simple truncation.

No, it is necessary to implement round to even. It allows to know whether you're in the 'half-way' situation, where the lost digits are exactly 500.
The 10/13 calculators (41,15) did not implement round-to-even, but simple round up (well 'away from zero'). Whether the lost digits are 500 or 500+ makes no difference there, see Tomas' examples.
I know nothing of the accuracy of the 15-digit functions I'm afraid.
Werner