HP41: accuracy of 13-digit routines

[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