HP41: accuracy of 13-digit routines

[Á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