Gamma Function by Stieltjes Continued Fraction
09-02-2015, 10:53 AM (This post was last modified: 09-02-2015 08:58 PM by Dieter.)
Post: #2
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Gamma Function by Stieltjes Continued Fraction
(09-01-2015 10:31 PM)lcwright1964 Wrote:  In high precision environments where adequate guard digits are not an issue (e.g., quadruple precision on Free42 Decimal), the approximation gives 10 digits easily (a minimum of about 10.4 exact decimal digits [edd] to be precise), and usually 11 or more across the argument range 0 to 70.

The error decreases as x increases. So the largest errors occur with small x, like 0...7 – 12 where they may cause the 10th digit to round incorrectly. For x>13 the relative error drops below 1E–11 so that the 10-digit result shoud be within ±0,6 ULP (like the built-in transcendental functions).

(09-01-2015 10:31 PM)lcwright1964 Wrote:  I am keen on suggestions to save steps, because I would like to be able to fit the routine into the 49 steps of the HP33C/E or the HP25.

Here is a first idea. The two multiplications after step 40 obviously avoid overflow for large arguments. But they can be coded more straightforward:

Code:
 ...  40 Y^X  41 ENTER  42 RDN  43 *  44 R^  45 *  ...

Code:
 ...  40 Y^X  41 *  42 LastX  43 *  ...

Also on other calculators the final 0 X<>M may be replaced by a simple RCL M (resp. RCL 0), since the zero only clears the 41's M-register to avoid leaving strange characters in Alpha.

And finally here is another possible improvement in terms of accuracy. The crucial point of the Stieltjes method is the exponential of a number close to x (cf. step 33). If x exceeds 10, the argument of the e^x function is >10 so that 1 ULP equals 1 E–8, which in turn means that the result of e^x has a potential relative error of 1 E–8. So at this point we may lose up to two digits!

A possible solution is...

Code:
...  30 1/X  31 E^X  32 X<>Y  33 E^X  34 /  ...

...this method, i.e. instead of e^(a–x) evaluate e^a / e^x.

Try x=30, 50 or 70 with the original code and compare with the result of the proposed method: In both cases the latter yields one more valid digit, here it is even dead on in two out or three cases. You may even tweak the results slightly by replacing the multiplication by sqrt(2 Pi) with a division by sqrt(0,5/Pi) – this constant rounds better to 10 digits. ;-)

The only drawback is that for integer arguments below 10 the results are off in the last digit while the original code happens to round nicely to the exact integer solution. But then try x=10, 11, 12, 13, 14... Although I haven't done any comprehensive tests, it looks like the proposed changes may yield results with an error within 10 ULP or 1 unit in the 9th digit, even on a 10-digit calculator like the 41 or 67/97. But be sure to do your own tests.

Addendum: here are some results from GMST and a modified version.

Code:
  x         GMST           ULP    modified        ULP   exact ------------------------------------------------------------------------   2      1,000000000         0   1,000000000        0   1   3,5    3,323350973        +3   3,323350973       +3   3,323350970   4,1    6,812622860        -3   6,812622865       +2   6,812622863   5      24,00000000         0   24,00000001       +1   24   7      720,0000000         0   720,0000003       +3   720   8,3    9281,392518        -8   9281,392522       -4   9281,392526  10      362879,9998        -2   362880,0000        0   362880  13      479001598,3       -17   479001600,2       +2   479001600  14,567  2,756346267 E+10  +12   2,756346254 E+10  -1   2,756346255 E+10  20      1,216451005 E+17   +1   1,216451005 E+17  +1   1,216451004 E+17  33,3    7,487577618 E+35  +21   7,487577596 E+35  -1   7,487577597 E+35  40,5    1,286050246 E+47   -2   1,286050248 E+47   0   1,286050248 E+47  50      6,082818622 E+62  -24   6,082818646 E+62  +6   6,082818640 E+62  60      1,386831191 E+80   +6   1,386831185 E+80   0   1,386831185 E+80  70      1,711224527 E+98   +3   1,711224524 E+98   0   1,711224524 E+98  70,956  9,933230633 E+99  +14   9,933230624 E+99  +5   9,933230619 E+99  70,957  9,975586583 E+99  -19   9,975586601 E+99  -1   9,975586602 E+99

Dieter
 « Next Oldest | Next Newest »

 Messages In This Thread Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-01-2015, 10:31 PM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-02-2015 10:53 AM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-03-2015, 08:10 PM RE: Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-04-2015, 01:54 AM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-04-2015, 06:29 AM RE: Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-04-2015, 06:32 PM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-05-2015, 05:49 AM RE: Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-05-2015, 09:55 PM RE: Gamma Function by Stieltjes Continued Fraction - Gerald H - 09-06-2015, 06:40 AM RE: Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-06-2015, 03:07 PM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-09-2015, 09:21 PM RE: Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-10-2015, 10:51 PM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-11-2015, 08:14 PM RE: Gamma Function by Stieltjes Continued Fraction - lcwright1964 - 09-13-2015, 02:06 AM RE: Gamma Function by Stieltjes Continued Fraction - Dieter - 09-13-2015, 01:55 PM

User(s) browsing this thread: 1 Guest(s)