HP41: accuracy of 13-digit routines

[Dieter]

It is well known that the HP41 internally uses extended precision routines with 13 significant digits, so that the finally returned 10-digit result usually can be considered accurate in all its ten digits.

But what about the 13-digit routines? I read Ángel Martin's PDF on the "HP41 13-digit 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 13-digit 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 13-digit precision ...but not always 13-digit accuracy. I also own a non-HP RPN calculator (an SR-54NC) with 12-digit 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

[Thomas Klemm]

(09-01-2015 09:23 PM)Dieter Wrote:  In other words: do they really provide correctly rounded 13-digit 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 13-digit result isn't rounded correctly but cut off.

As far as I remeber rounding is just adding 500 to the 13-digit mantissa.
Thus in the example from above:

Code:

10.49954995499
+          500
10.49954995999

From this only the first 10 digits are used.

HTH
Thomas

[Joe Horn]

Thomas is correct. An example that's easier to key in is: \(\sqrt {5} \), which is equal to 2.236067977 499 7896964..., which the HP-41 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 10-digit models, but is still used by the newer 12-digit HP calculators, which calculate internally to 15 *truncated* digits, then round that result to 12 digits.

[Dieter]

(09-02-2015 04:45 AM)Thomas Klemm Wrote:  Thus the internal 13-digit result isn't rounded correctly but cut off.
...
From this only the first 10 digits are used.

(09-02-2015 06:11 AM)Joe Horn Wrote:  Thomas is correct.
...
This "truncate, then round" procedure was not only used by the old 10-digit models, but is still used by the newer 12-digit HP calculators, which calculate internally to 15 *truncated* digits, then round that result to 12 digits.

Thank you very much. So the 13-digit 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 13-digit 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

[Ángel Martin]

(09-02-2015 10:18 AM)Dieter Wrote:  Thank you very much. So the 13-digit 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 13-digit routines?

Typically that seems to be the case: the 13-digit result is held in the {A,B} registers, with the 13-digit 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 10-digit 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

[Michael Fehlhammer]

(09-02-2015 10:18 AM)Dieter Wrote:  Thank you very much. So the 13-digit routines truncate after 13 digits. I assume this also applies to the transcendental functions (log, exp, trig etc.).

Dieter

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.

[Ángel Martin]

(09-02-2015 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 X-register 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

[Dieter]

(09-02-2015 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 (...).

Operations that can suffer a slightly larger error, but still significantly smaller than one unit in the 10th significant digit of the result, include ∆%, →H, →RAD, →DEG, Px,y and Cx,y, LN, LOG, 10^x, ASIN, ACOS, ATAN (...) and finally SIN, COS, and TAN for real arguments in Degrees and Grads modes (...).

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. (...) Functions so affected, are e^x, y^x, x! (for non-integer x) (...).

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 nearly-exact 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 13-digit 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.

(09-02-2015 12:03 PM)Ángel Martin Wrote:  Yep, and even if they were how would you output it using just the standard X-register approach?

That's a different story. The question is what accuracy level can be expected from the internal 13-digit routines. Since some functions may not return more than 11 valid digits, this seems to be the limit for MCode using these limited-accuracy functions.

What do you think?

Dieter

[Ángel Martin]

(09-02-2015 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.

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 nearly-exact 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.

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 [SQR-13] followed by a 13-digit 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.

(09-02-2015 06:07 PM)Dieter Wrote:  What does this mean for a Gamma implementation in 13-digit 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 non-integer arguments. BTW there's also exponentials in the Lanczos formula, as well as sine factors in the reflection formula.

Cheers,
'AM

[Dieter]

(09-03-2015 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 [SQR-13] followed by a 13-digit 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 10-digit 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 [SQR-13] 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.

(09-03-2015 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.

(09-03-2015 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 non-integer 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.

(09-03-2015 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

[Ángel Martin]

(09-03-2015 12:36 PM)Dieter Wrote:  

(09-03-2015 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.

(09-03-2015 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 non-integer 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. ;-)

Great!

(09-03-2015 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/the-library/11

ÁM

[Dieter]

(09-03-2015 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.

(09-03-2015 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.

(09-03-2015 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)^(x-0,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 shift-and-divide 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...13-digit 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 13-digit 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 13-digit 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 exp-term. 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. #-)

(09-03-2015 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:
http://www.rskey.org/CMS/index.php/the-library/11

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

[Ángel Martin]

[attachment=2492]

(09-03-2015 06:17 PM)Dieter Wrote:  ... the basic formula is:

Γ(x) = [a0 + a1/x + a2/(x+1) + a3/(x+2) + a4/(x+3)] * ((x+c)/e)^(x-0,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) ]

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...13-digit 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 13-digit 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.

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:


X        SandMath           Dieter                  Wolfram-Alpha

-2.5        -0.9453087210    -0.9453087206        -0.9453087205
-1.5        2.3632731802    2.3632718000        2.3632718012
-0.5        -3.5449022520    -3.5449077020        -3.5449077018
0.1        9.513507699    9.5135077160        9.5135076987
0.5        1.7724538510    1.7724538520        1.7724538509
1.5        8.862269255 E-01    8.862269254 E-01        8,8622692545 E-01
5.5        5.234277778 E01    5.234277778 E01        5,234277778 E01
10.5        1.133278389 E06    1.133278389 E06        1.133278389 E06
20.5        5.406242982 E17    5.406242982 E17        5.406242982 E17
50.5        4.290462912 E63    4.290462912 E63        4.290462912 E63
70.5        1.429158828 E99    1.429158828 E99        1.429158828 E99

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

[Dieter]

(09-04-2015 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. :-)

(09-04-2015 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.

(09-04-2015 12:07 PM)Ángel Martin Wrote:  (as usual I can't get the table format right...

That's easy. Text bracketed by code-tags 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.)
------------------------------------------------------------------------------
-2.5      -9.453087210 E-01  -9,453087206 E-01    -9.453087205 E-01 (...0483)
-1.5       2.363271802        2.363271800          2.363271801      (...0121)
-0.5      -3.544907702       -3.544907702         -3.544907702      (...0181)
 0.1       9.513507699        9.513507716          9.513507699      (...9867)
 0.5       1.772453851        1.772453852          1.772453851      (...5091)
 1.5       8.862269255 E-01   8.862269254 E-01     8,862269255 E-01 (...5453)
 5.5       5.234277778 E+01   5.234277778 E+01     5,234277778 E+01 (...7846)
10.5       1.133278389 E+06   1.133278389 E+06     1.133278389 E+06 (...8895)
20.5       5.406242982 E+17   5.406242982 E+17     5.406242982 E+17 (...8234)
50.5       4.290462912 E+63   4.290462912 E+63     4.290462912 E+63 (...1235)
70.5       1.429158828 E+99   1.429158828 E+99     1.429158828 E+99 (...2830)

(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).

(09-04-2015 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

[Ángel Martin]

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

[Dieter]

(09-04-2015 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
a0 =  3.2743 15106 67  E-02
a1 =  1.1882 42945 545
a2 = -1.0667 84960 2
a3 =  1.8797 42975 77  E-01
a4 = -2.7122 78856     E-03

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.

(09-04-2015 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

[Dieter]

(09-01-2015 09:23 PM)Dieter Wrote:  [ HP41 13-digit routines ]
In other words: do they really provide correctly rounded 13-digit results under all circumstances? Or may the results be truncated instead of rounded? Or maybe even the last digit may be slightly off?

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:

• The four arithmetic functions seem to return the exact result, truncated after 13 digits. You can actually see how one digit of the result after another is calculated. So the error is +0 and –1 unit of the last digit.
• The same seems to apply to the SQRT function.
• The natural logarithm LN also seems to be exact, or at least nearly so. In some cases the result is even better than expected: ln 4 = 1,38629436111989... Here the HP41 yields 1,386294361120, i.e. the true result correctly rounded up. This seems to happen because internally ln 2 = 0,6931471805600 is used instead of 0,693147180559945...
• The base-10 logarithm LOG is evaluated as ln x / ln 10 where the latter equals 2,302585092994 (this value rounds very nicely to 13 digits). So the error should not be larger than it is for the natural log. But there are cases where log10 has a smaller magnitude than ln. Here one ULP of ln means 10/ln10 = 4,3 ULP of the log10 result. Since the final division is truncated, the result may be off by one more ULP. As far as I can see this can cause errors up to –5,3 ULP of the 13-digit result for the base-10 logarithm. For instance, lg 8,001 = 0,903144270409538...is returned as 0,9031442704091. The critical domains for such cases are e¹...10¹, e¹⁰...10¹⁰, ... etc.
• The power function y^x is evaluated as e^x · ln y. So the accuracy limits of the log and exp functions apply, but there is another problem: results larger than e¹⁰⁰ (approx. 2,688 E–43) require an intermediate result beyond 100, i.e. the natural log of the final result, even if evaluated as assurately as possible, may be off by 1 E–10. Which in turn means a relative error of the same order in the final result. Since an additional error may occur during the multiplication with the exponent, this may lead to an error of ln 10 = 2,3... = two or maybe three units in the last place of the 10-digit result.This is also mentioned in the HP15C Advanced Functions Handbook.
  

I did not do more sophisticated tests, and I am only beginning to understand what's going on in a HP41 during such calculations, but maybe the MCode experts may provide some more information. And/or corrections. ;-)

Dieter

[Dieter]

(09-02-2015 04:45 AM)Thomas Klemm Wrote:  As far as I remeber rounding is just adding 500 to the 13-digit mantissa.
Thus in the example from above:

Code:

10.49954995499
+          500
10.49954995999

From this only the first 10 digits are used.

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
=> 10.499 54995 998
=> 10.499 54995

Dieter

[Thomas Klemm]

(09-25-2015 10:46 PM)Dieter Wrote:  it seems that rounding is done by doubling the last three (guard) digits:

Code:

   10.499 54995 499
=> 10.499 54995 998
=> 10.499 54995

Yeah, that's probably easier than adding 500. The result is of course the same. Thanks for the refresher.

Cheers
Thomas