HP41: accuracy of 13digit routines

09012015, 09:23 PM
Post: #1




HP41: accuracy of 13digit routines
It is well known that the HP41 internally uses extended precision routines with 13 significant digits, so that the finally returned 10digit result usually can be considered accurate in all its ten digits.
But what about the 13digit routines? I read Ángel Martin's PDF on the "HP41 13digit OS routines", but while I found lots of useful information, nothing was said about the actual accuracy of these routines. In other words: do they really provide correctly rounded 13digit results under all circumstances? Or may the results be truncated instead of rounded? Or maybe even the last digit may be slightly off? Here I remember the good old TI58/59 with their 13digit precision ...but not always 13digit accuracy. I also own a nonHP RPN calculator (an SR54NC) with 12digit precision (and display), but the result sometimes does not carry more than 11 or even 10 valid digits. So my question is: how does the HP41 perform in this regard? Dieter 

09022015, 04:45 AM
Post: #2




RE: HP41: accuracy of 13digit routines
(09012015 09:23 PM)Dieter Wrote: In other words: do they really provide correctly rounded 13digit results under all circumstances? Or may the results be truncated instead of rounded? Or maybe even the last digit may be slightly off? This depends a little on the definition of "correctly rounded". You may try 104985/9999 which is: 10.499549954995... Correctly rounded to 13 digits would give us: 10.49954995500 If we'd round this to 10 digits we'd get: 10.49954996 But that's not what we get. Instead we get: 10.49954995 Thus the internal 13digit result isn't rounded correctly but cut off. As far as I remeber rounding is just adding 500 to the 13digit mantissa. Thus in the example from above: Code: 10.49954995499 HTH Thomas 

09022015, 06:11 AM
Post: #3




RE: HP41: accuracy of 13digit routines
Thomas is correct. An example that's easier to key in is: \(\sqrt {5} \), which is equal to 2.236067977 499 7896964..., which the HP41 says is 2.236067977, which shows that it did NOT round to the 13th digit, but truncated at the 13th digit, then rounded THAT to 10 digits.
This "truncate, then round" procedure was not only used by the old 10digit models, but is still used by the newer 12digit HP calculators, which calculate internally to 15 *truncated* digits, then round that result to 12 digits. <0ɸ0> Joe 

09022015, 10:18 AM
Post: #4




RE: HP41: accuracy of 13digit routines
(09022015 04:45 AM)Thomas Klemm Wrote: Thus the internal 13digit result isn't rounded correctly but cut off. (09022015 06:11 AM)Joe Horn Wrote: Thomas is correct. Thank you very much. So the 13digit routines truncate after 13 digits. I assume this also applies to the transcendental functions (log, exp, trig etc.). Now if I set a WP34s to RM 4 (truncate) and apply a RSD 13 to every result, this should yield the same result as the 41's 13digit routines? The reason I'm asking this is testing a Lanczos Gamma approximation with an optimized coefficient set, cf. the recent thread in the HP41 software library. Dieter 

09022015, 11:05 AM
(This post was last modified: 09022015 11:07 AM by Ángel Martin.)
Post: #5




RE: HP41: accuracy of 13digit routines
(09022015 10:18 AM)Dieter Wrote: Thank you very much. So the 13digit routines truncate after 13 digits. I assume this also applies to the transcendental functions (log, exp, trig etc.). Now if I set a WP34s to RM 4 (truncate) and apply a RSD 13 to every result, this should yield the same result as the 41's 13digit routines? Typically that seems to be the case: the 13digit result is held in the {A,B} registers, with the 13digit mantissa in B<12:0> and the sign and exponent in A. Once the math routine is done it transfers the mantissa from B to the C register, most usually using the instruction C=B M  thus it is a 10digit version made by direct truncation. I however have done some rounding *prior* to this truncation in some of my own math functions  hard to remember which ones though... Cheers, ÁM "To live or die by your own sword one must first learn to wield it aptly." 

09022015, 11:19 AM
Post: #6




RE: HP41: accuracy of 13digit routines
(09022015 10:18 AM)Dieter Wrote: Thank you very much. So the 13digit routines truncate after 13 digits. I assume this also applies to the transcendental functions (log, exp, trig etc.).The only reason why the 41C uses 13 digits internally ( although registers can store 10 digits only ) is to provide 10 digit precision for the results of the more complex functions; due to the used algorithms and error propagation you need these extra 3 digits. Therefore you should not assume that results of log, sin etc. yield more than 10 correct digits. 

09022015, 12:03 PM
(This post was last modified: 09022015 12:04 PM by Ángel Martin.)
Post: #7




RE: HP41: accuracy of 13digit routines
(09022015 11:19 AM)Michael Fehlhammer Wrote: [The only reason why the 41C uses 13 digits internally ( although registers can store 10 digits only ) is to provide 10 digit precision for the results of the more complex functions; due to the used algorithms and error propagation you need these extra 3 digits. Therefore you should not assume that results of log, sin etc. yield more than 10 correct digits. Yep, and even if they were how would you output it using just the standard Xregister approach? Of course HP could have added "double precision" routines like the 16C does for the division and multiplication  where the output would use two stack registers, logically X and Y but that's a different story... ÁM "To live or die by your own sword one must first learn to wield it aptly." 

09022015, 06:07 PM
(This post was last modified: 09032015 12:39 PM by Dieter.)
Post: #8




RE: HP41: accuracy of 13digit routines
(09022015 11:19 AM)Michael Fehlhammer Wrote: ... Therefore you should not assume that results of log, sin etc. yield more than 10 correct digits. Wait a minute. There is a chapter on accuracy in the HP15C Advanced Functions Handbook. The internal algorithms of the 41C and 15C arithmetic and transcendental functions should be very similar at least, if not identical, so let's assume the following also applies to the 41 series. The appendix (chapter "A Hierarchy of Errors") distinguishes several error levels, among them these three: Quote:Operations that deliver "correctly rounded" results whose error cannot exceed ½ unit in their last (10th) significant digit include the real algebraic operations +, –, *, ÷, x², √x, 1/x, and %, (...) and →H.MS. These results are the best that 10 significant digits can represent (...). So the accuracy level depends on the function, and sometimes even on its argument. Since the four basic arithmetic functions as well as x² and √x deliver correctly rounded results they must be virtually exact internally, i.e. in all 13 internal digits, maybe the last one truncated. This would yield a max. error of 0,501 ULP. The logs, 10^x and the trig functions (in DEG/GRAD) and their inverses may exhibit a slightly larger error, yet "significantly less" than 1 ULP. Maybe we can translate this to 0,6 ULP or 1 unit in the 11th place, or a relative error within 1E–11, maybe 0,7 ULP (rel. error < 2E–11). So there are 11 exact or nearlyexact digits. And finally the error may spill into the last of the 10 digits presented to the user, where the 10th digit may be off by 2 or even 3 units. What does this mean for a Gamma implementation in 13digit MCode? The functions used in the Lanczos algorithm are mostly +, –, x and / which are (nearly) exact. On the other hand the logs may only have 11 valid digits. So this is what can be expected: about 11 exact digits, and a final error similar to that of the internal log or trig functions. (09022015 12:03 PM)Ángel Martin Wrote: Yep, and even if they were how would you output it using just the standard Xregister approach? That's a different story. The question is what accuracy level can be expected from the internal 13digit routines. Since some functions may not return more than 11 valid digits, this seems to be the limit for MCode using these limitedaccuracy functions. What do you think? Dieter 

09032015, 06:22 AM
(This post was last modified: 09032015 06:23 AM by Ángel Martin.)
Post: #9




RE: HP41: accuracy of 13digit routines
(09022015 06:07 PM)Dieter Wrote: So the accuracy level depends on the function, and sometimes even on its argument. Since the four basic arithmetic functions as well as x² and √x deliver correctly rounded results they must be virtually exact internally, i.e. in all 13 internal digits, maybe the last one truncated. This would yield a max. error of 0,501 ULP. Sure, I meant that the accuracy will always be correct to the 10th decimal digit but I didn't imply that all algorithms would be equally effective in the digits beyond that. Definitely the transcendental functions are more complicated to implement, but not only those show weak points. I don't know if there were changes made to the 15C MCODE, but I can tell you that on the 41 the square root at least is far from being that accurate. If you do a chained calculation in MCODE with [SQR13] followed by a 13digit multiplication you not always get the same original argument to the 13th. digit.... and there's NOT such a thing as X^2 as individual function either, it uses the same multiplication code. (09022015 06:07 PM)Dieter Wrote: What does this mean for a Gamma implementation in 13digit MCode? The functions used in the Lanczos algorithm are mostly +, –, x and / which are (nearly) exact. On the other hand the logs may only have 11 valid digits. So this is what can be expected: about 11 exact digits, and a final error similar to that of the internal log or trig functions. the X! on the 15C is the gamma function, right? So there you also have it to test for the accuracy to the 10th decimal digit. As far as GAMMA in the SandMath I checked the accuracy for the natural arguments were ok (see the table in its manual) but didn't do any further checks for noninteger arguments. BTW there's also exponentials in the Lanczos formula, as well as sine factors in the reflection formula. Cheers, 'AM "To live or die by your own sword one must first learn to wield it aptly." 

09032015, 12:36 PM
(This post was last modified: 09032015 12:36 PM by Dieter.)
Post: #10




RE: HP41: accuracy of 13digit routines
(09032015 06:22 AM)Ángel Martin Wrote: I don't know if there were changes made to the 15C MCODE, but I can tell you that on the 41 the square root at least is far from being that accurate. If you do a chained calculation in MCODE with [SQR13] followed by a 13digit multiplication you not always get the same original argument to the 13th. digit.... Well, I hope so. ;) Squaring a root that is exact in all digits does not neccessarily yield the original number. For instance, there simply is no 10digit value for √2 that, if squared, would return 2 again. You get either 1,999999999 or 2,000000002. With 13 digits the same applies to √7. The (truncated) root is 2,645751311064, and if [SQR13] returns this result, everything is fine. Even though this value squared gives 6,999999999996 resp. ...997. Even if the root was correctly rounded to 13 digits instead of truncated, the square of this would not return 7 but 7,000000000002. (09032015 06:22 AM)Ángel Martin Wrote: the X! on the 15C is the gamma function, right? So there you also have it to test for the accuracy to the 10th decimal digit. Yes, on the 15C x! actually evaluates Gamma(x+1), just like the 34C did before. The 15C AFH says that the last digit may be off by 1 ULP, maybe 2 for large results. (09032015 06:22 AM)Ángel Martin Wrote: As far as GAMMA in the SandMath I checked the accuracy for the natural arguments were ok (see the table in its manual) but didn't do any further checks for noninteger arguments. Here and there the result may be off in the last digit, but the cases I found were just 1 ULP high or low. So it seems on par with the 15C. ;) BTW, have you read the recent thread on Spouge's Gamma approximation in the HP41 software library? Les and I have been talking about an improved Lanczos method, and I posted a set of coefficents resulting in a relative error less than 2 E–11 up to Gamma(71), and this with even two terms less than the current Sandmath implementation. (09032015 06:22 AM)Ángel Martin Wrote: BTW there's also exponentials in the Lanczos formula, as well as sine factors in the reflection formula. Yes. An that's what will limit the accuracy of the final result. Dieter 

09032015, 01:33 PM
(This post was last modified: 09032015 02:47 PM by Ángel Martin.)
Post: #11




RE: HP41: accuracy of 13digit routines
(09032015 12:36 PM)Dieter Wrote:(09032015 06:22 AM)Ángel Martin Wrote: the X! on the 15C is the gamma function, right? So there you also have it to test for the accuracy to the 10th decimal digit. Great! (09032015 12:36 PM)Dieter Wrote: BTW, have you read the recent thread on Spouge's Gamma approximation in the HP41 software library? Les and I have been talking about an improved Lanczos method, and I posted a set of coefficents resulting in a relative error less than 2 E–11 up to Gamma(71), and this with even two terms less than the current Sandmath implementation. I need to check how that one will pan out  how many coefficients are you using in the inproved version? I haven't really looked into the details yet... but if it's equal or better than the current implementation (up to 10 digits, remember...) then it'll be worth replacing it with the new set. PS. ok there are 4 coefficients instead of 7  so far so good, but what about the rest of the formula? c replaces the arbitrary "5" value, but does the rest remain unchanged?? Pls. let me know which of the formulas on Viktor's page would be the applicable one for the new coefficient set: http://www.rskey.org/CMS/index.php/thelibrary/11 ÁM "To live or die by your own sword one must first learn to wield it aptly." 

09032015, 06:17 PM
(This post was last modified: 09032015 06:58 PM by Dieter.)
Post: #12




RE: HP41: accuracy of 13digit routines
(09032015 01:33 PM)Ángel Martin Wrote: I need to check how that one will pan out  how many coefficients are you using in the inproved version? I have posted several coefficient sets for n=4 and n=5, i.e. 5 or 6 coefficients plus another one called c (3,8... resp. 5,081) which requires less precision. (09032015 01:33 PM)Ángel Martin Wrote: I haven't really looked into the details yet... but if it's equal or better than the current implementation (up to 10 digits, remember...) then it'll be worth replacing it with the new set. The n=4 set should be on par, maybe slightly better. With sufficient working precision the relative error for x = 1...71 is within ±1,6 E–11. The n=5 set should yield approx. two more digits. With sufficient working precision the relative error for x = 1...71 is within ±3,8 E–13. (09032015 01:33 PM)Ángel Martin Wrote: PS. ok there are 4 coefficients instead of 7  so far so good, but what about the rest of the formula? c replaces the arbitrary "5" value, but does the rest remain unchanged?? The formula can be found in post #23 and #24 ff. in the Spouge thread. Caution: Les' formula calculates Gamma, while mine is designed for x!, so the argument differs by one. Anyway, the basic formula is: Γ(x) = [a0 + a1/x + a2/(x+1) + a3/(x+2) + a4/(x+3)] * ((x+c)/e)^(x0,5) In an MCODE implementation, the power expression on the right may be rewritten: Γ(x) = [a0 + a1/x + a2/(x+1) + a3/(x+2) + a4/(x+3)] * exp[ (ln(x+c) – 1) * (x–0,5) ] This formula has two nice properties: It does not require any shiftanddivide operations for small x, and – at least in the cases discussed here – it does not overflow for large x. Here a0...a4 are the mentioned coefficients, and c is another constant that is approx. 3,8... for the n=4 case. A careful choice of c can substantially reduce the error. The optimized 12...13digit coefficients for n=4 and c = 3,8306414 are listed in post #33. These should yield a relative error not much more than 1,6 E–11. Since the 13digit MCODE routines seem to truncate their results (instead of rounding them) the error may be somewhat larger, resp. the "perfect" coefficients may slightly differ from those posted. There also is a coefficient set for n=5 in post #32, leading to a relative error less than 4 E–13 (if evaluated with sufficient precision). The values are given to 24 digits, maybe I can provide an optimized 13digit set. However, in real life the final accuracy is limited by a simple phenomenon, both in this approximation as well as in others. The crucial point is the power within the formula, resp. the expterm. Let's assume the argument of the exponential function has 13 exact digits. This value may be greater than 100 (about 230, to be precise), so even a correctly rounded resp. truncated value has a potential error of ±5 E–11 resp. ±1 E–10. After the exponentiation, this means quite exactly a relative error of the same magnitude. In other words: with three guard digits and results as large as 9 E+99 it does not get better than 10 digits. #) (09032015 01:33 PM)Ángel Martin Wrote: Pls. let me know which of the formulas on Viktor's page would be the applicable one for the new coefficient set: The original formula (with slightly different coefficients) was posted by Les, probably based on the Pugh paper mentioned earlier in that thread. Is looks like a rescaled, slighty modified Lanczos approximation. Dieter 

09042015, 12:07 PM
(This post was last modified: 09042015 12:28 PM by Ángel Martin.)
Post: #13




RE: HP41: accuracy of 13digit routines
[attachment=2492]
(09032015 06:17 PM)Dieter Wrote: ... the basic formula is: I finally got around to doing the changes in the GAMMA code, using the coefficients you came up with. See the attached MOD file with the modified verison of the function already compiled. Indeed the results are very much the same, with the exception of points around x=1, see the comparison table that also has the WolframAlpha colum as 'exact' answers: (as usual I can't get the table format right... Code:
The new code is 44 bytes shorter and on average it takes about 0.3 seconds less to execute than the original function based on Viktor's formula with 7 coefficients. Why it isn't more of a difference can be explained by the fact that even if the original had two more terms it was easier to calculate using Honer's method  which cannot be applied in the modified version. So all in all it's better  with a caveat in that region mentioned before. No idea why would that be the case, right? ÁM "To live or die by your own sword one must first learn to wield it aptly." 

09042015, 07:36 PM
(This post was last modified: 09042015 09:57 PM by Dieter.)
Post: #14




RE: HP41: accuracy of 13digit routines
(09042015 12:07 PM)Ángel Martin Wrote: I finally got around to doing the changes in the GAMMA code, using the coefficients you came up with. See the attached MOD file with the modified verison of the function already compiled. Ah, thank you very much. :) (09042015 12:07 PM)Ángel Martin Wrote: Indeed the results are very much the same, with the exception of points around x=1, see the comparison table that also has the WolframAlpha colum as 'exact' answers: Hm... for values below 1 the accuracy is not what I would have expected. But I think there is a simple reason for this. See below. (09042015 12:07 PM)Ángel Martin Wrote: (as usual I can't get the table format right... That's easy. Text bracketed by codetags is formatted with a fixed font. However, the message editor is always set to its standard proportional font... #) So you simply have to use another editor (notepad.exe will do), format your code section there and copy/paste it back to the message editor. At least that's the way I do it. This way you get... Code: x Sandmath Dieter Exact (10 dig.) (12 dig.) (note: your table listed a quite inaccurate Sandmath result for x=–0.5 while actually it is dead on. The new table lists the correct result). (09042015 12:07 PM)Ángel Martin Wrote: So all in all it's better  with a caveat in that region mentioned before. No idea why would that be the case, right? Well, actually I think I know what causes the problem... #) The original coefficient set was optimized for the factorial function with 0 ≤ x ≤ 70. For the Gamma function this translates to 1 ≤ x ≤ 71. So arguments less than 1 will cause substantially larger errors than those within the optimized domain. That's why it performs so poorly between x=0 and 1, for instance at x=0.1. On the other hand Gamma(1.1)/0.1 yields a correct result, as well as Gamma(1.5)/0.5 which now should return 1.772453851. I assume the reflection formula is used for x < 0.5 (where x = 1–x). So we actually need a new coefficient set that is optimized for 0.5 ≤ x ≤ 71. Sorry, I just overlooked this point. I'll see what I can come up with. This means some more hours with Excel and Free 42... #) Maybe you can do some more tests for x>1 or x<0 and compare the results with Sandmath's. The error maxima are near x=1.245, 2.395, 5.15, 12.25 and 34. Update: I now have set up a new approximation for 0.5 ≤ x ≤ 71 with a relative error of approx. 3 E–11. Dieter 

09042015, 09:53 PM
(This post was last modified: 09042015 09:53 PM by Ángel Martin.)
Post: #15




RE: HP41: accuracy of 13digit routines
it's good to have you there watching over  sorry about the bogus data point in x=0.5, it turns out i was using an older version of the SandMath on my CL  version control is not my forte... this weekend I'll update the CL flash library witht he latest revisions.
Great to hear you're on to an updated coefficient set  I look forward to it and getting the code squeaky clean with the new values  those .3 seconds and 44 bytes do make a difference  remember than GAMMA is used as subroutine in several other functions in the module, like Bessel and a few more. 'AM "To live or die by your own sword one must first learn to wield it aptly." 

09042015, 10:24 PM
(This post was last modified: 09052015 05:51 AM by Dieter.)
Post: #16




RE: HP41: accuracy of 13digit routines
(09042015 09:53 PM)Ángel Martin Wrote: Great to hear you're on to an updated coefficient set  I look forward to it and getting the code squeaky clean with the new values This is not the final result, but for the time being you may try these values: Code: c = 3.838 Evaluated exactly, the relative error should be less than 3,2 E–11 for x=0,5...71. At x=0.5 the result should even be more or less exact. Without this restriction the error may drop slightly below 3 E–11. (09042015 09:53 PM)Ángel Martin Wrote:  those .3 seconds and 44 bytes do make a difference  remember than GAMMA is used as subroutine in several other functions in the module, like Bessel and a few more. Hmmm... if GAMMA is used by other functions maybe a bit more accuracy would be helpful. Which means one more term and an estimated relative error near 1E–12. Dieter 

09252015, 01:53 PM
(This post was last modified: 09252015 01:57 PM by Dieter.)
Post: #17




RE: HP41: accuracy of 13digit routines
(09012015 09:23 PM)Dieter Wrote: [ HP41 13digit routines ] I now have played around a bit with V41's MCode console. Slowly, very slowly I am beginning to get a first idea of how it works. Here are some first insights:
Dieter 

09252015, 10:46 PM
Post: #18




RE: HP41: accuracy of 13digit routines
(09022015 04:45 AM)Thomas Klemm Wrote: As far as I remeber rounding is just adding 500 to the 13digit mantissa. Looking at the V41 MCode console output, it seems that rounding is done by doubling the last three (guard) digits: Code: 10.499 54995 499 Dieter 

09252015, 11:33 PM
Post: #19




RE: HP41: accuracy of 13digit routines  
« Next Oldest  Next Newest »

User(s) browsing this thread: