HP41: accuracy of 13-digit routines

[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