HP41: accuracy of 13-digit routines

[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