Accuracy and the power function

02162014, 02:43 PM
(This post was last modified: 02162014 05:31 PM by Dieter.)
Post: #1




Accuracy and the power function
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 12digitvalue agrees with the true value rounded to 13 digits, so that the relative error in \(\large\pi\) is less than 6,6 E14.
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 12digit 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 12digit HP and see what you get. Dieter 

02162014, 03:11 PM
Post: #2




RE: Accuracy and the power function
(02162014 02:43 PM)Dieter Wrote: Let's assume \(2\large\pi\) = 6,28318530718, i.e. the best possible 12digit value. Now evaluate some powers of this: 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 12digit HP and see what you get. On the HP 50g, the errors on these calculations are 1 ULP, 1 ULP, and 4 ULP, respectively. <0ɸ0> Joe 

02162014, 05:45 PM
Post: #3




RE: Accuracy and the power function
(02162014 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 

02162014, 10:11 PM
Post: #4




RE: Accuracy and the power function
Methinks it may be explained by the HP15C 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 10digit mantissa, with intermediate internal math using 13digit 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^{20}; 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^{201}, which is calculated as 7.968419664 × 10^{95}. The last digit 4 should be 6 instead ... [jump to page 183] To understand the error in 3^{201}, 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. <0ɸ0> Joe 

02182014, 08:28 PM
Post: #5




RE: Accuracy and the power function
(02162014 10:11 PM)Joe Horn Wrote: Methinks it may be explained by the HP15C 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 10digit 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 13digit internal result is off by 1 ULP, the absolute error is 1E10. Which after the following exponentiation results in a relative error of 1E10 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 1E11, 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 1E12 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 

02192014, 12:33 AM
Post: #6




RE: Accuracy and the power function
(02182014 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. <0ɸ0> Joe 

02192014, 09:54 AM
(This post was last modified: 02192014 10:07 AM by Werner.)
Post: #7




RE: Accuracy and the power function
(02182014 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 1E12 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 15digit routines indeed are not exact to the last digit. In your example, %%LN(6.28318530718) is off by one ULP. Secondly, the 15digit routines (%%*, %%+, %%, %%/) do not round, they truncate. They are, after all, meant as intermediates to produce a correctly rounded 12digit result (for the elementary operations). Code:
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 15digit %%LN, 15digit %%EXP and 15digit rounding arithmetic would produce. Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L 

02192014, 11:36 AM
Post: #8




RE: Accuracy and the power function
(02192014 09:54 AM)Werner Wrote: In your example, %%LN(6.28318530718) is off by one ULP. Anything special you want to tell us using %% ?? d#? 

02192014, 11:47 AM
Post: #9




RE: Accuracy and the power function
For the uninitiated, the %%functions are the 48's internal 15digit versions of the User functions: LN, EXP, *, + , ...
Sorry! Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L 

02192014, 12:35 PM
Post: #10




RE: Accuracy and the power function
(02192014 09:54 AM)Werner Wrote: There are two things at work here.Ah, thank you very much, that's exactly what I have been looking for. I tried to simulate the error on a 34s in 16digit 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 Dieter 

02192014, 01:07 PM
Post: #11




RE: Accuracy and the power function
Truncation (and keeping track of lost nonzero digits with the Sticky Bit) in the 15digit routines is not only common practice, it's a necessity if you want to be able to produce a correctly rounded 12digit result for simple calculations.
Suppose the 15digit 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 16digit number rounded to 12 digits is 5.02000000121e15 Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L 

02202014, 05:17 PM
Post: #12




RE: Accuracy and the power function
(02192014 01:07 PM)Werner Wrote: Truncation (and keeping track of lost nonzero digits with the Sticky Bit) in the 15digit routines is not only common practice, it's a necessity if you want to be able to produce a correctly rounded 12digit 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 12digit 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 

02212014, 07:11 AM
Post: #13




RE: Accuracy and the power function
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 addandshift 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 nonzero 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 nonzero. 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 12digit HPs do get the correctly rounded answer for this example (or any other), yet they all truncate their 15digit result. Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L 

02232014, 02:19 PM
Post: #14




RE: Accuracy and the power function
(02212014 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 15digit 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 10digit devices with 13digit internal precision, for instance the 41series. Dieter 

02232014, 08:09 PM
(This post was last modified: 02242014 06:54 AM by Thomas Klemm.)
Post: #15




RE: Accuracy and the power function
(02232014 02:19 PM)Dieter Wrote: I would also be interested very much in information on this matter regarding 10digit devices with 13digit internal precision, for instance the 41series. Not for the HP41C but for the HP15C 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.
But this is handled separately and I have no trace added for this scenario. Though I didn't check the sources of the HP48 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 Add both numbers in a and b now that they are alligned: Code: a=00000000000499 b=05020000121000 c=00000000000513 Round half away from zero: Code: a=05020000121499 b=00000000000013 c=05020000121499 5.020000121e13 + 5009: Shift right register c. The 9 at the end of 5009 is lost: Code: a=05009000000513 b=00000000000513 c=00000000005009 Add both numbers in a and b now that they are alligned: Code: a=00000000000500 b=05020000121000 c=00000000000513 Round half away from zero: Code: a=05020000121500 b=00000000000013 c=05020000121500 5.020000124e13 + 4999: Shift right register c. The 9 at the end of 4999 is lost: Code: a=04999000000513 b=00000000000513 c=00000000004999 Add both numbers in a and b now that they are alligned: Code: a=00000000000499 b=05020000124000 c=00000000000513 Round half away from zero: Code: a=05020000124499 b=00000000000013 c=05020000124499 5.020000124e13 + 5009: Shift right register c. The 9 at the end of 5009 is lost: Code: a=05009000000513 b=00000000000513 c=00000000005009 Add both numbers in a and b now that they are alligned: Code: a=00000000000500 b=05020000124000 c=00000000000513 Round half away from zero: Code: a=05020000124500 b=00000000000013 c=05020000124500 

02232014, 09:06 PM
Post: #16




RE: Accuracy and the power function
(02232014 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 'halfway' situation, where the lost digits are exactly 500. The 10/13 calculators (41,15) did not implement roundtoeven, 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 15digit functions I'm afraid. Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)