Faster inverse gamma and factorial for the WP 34S

[Dieter]

(12-04-2014 06:14 AM)Bit Wrote:  I noticed that the inverse gamma library routine included with the 34S (discussed here) was quite slow for small inputs. I've added a very simple second algorithm that uses the secant method and converges faster for small inputs: < ~2.3 in single precision mode and < ~1e17 in double precision mode. For example, 1.6 vs 5.3 seconds in single precision mode or 2.1 vs 13.5 seconds in double precision mode with 0.89 as the input. While it's a noticeable improvement, it's probably far from optimal since I'm not an expert on numerical algorithms, so I'd be interested in better versions.

The standard lib that comes with the 34s firmware contains a digamma routine - so why not use the obvious and setup a simple Newton solver for Γ^-1(a):

x := x + (a/Γ(x) – 1) / Ψ(x)

The idea of using Ψ was already mentioned in the 2013 thread, but the solution posted there replaced Ψ with an approximation. Using "the real thing" seems to work better, at least for small arguments.

Since the digamma routine does not preserve the stack (at least the current version doesn't) two registers are used. R₀ stores the input a, R₁ holds the current approximation. The code is simple and straightforward:

Code:

LBL"IGM"
STO 00
#pi
#086
/
+
CNST 87   // that's sqrt(2*pi)
/
LN
RCL X
#eE
/
Wp
/
#1/2
+
LBL 00
STO 01
XEQ'Ψ'
RCL 00
RCL 01
Γ
-
RCL/L
RCL/Y
RCL 01
ENTER
RCL+Z
CNVG? 03
RTN
GTO 00

Timings:
An input of a = 0,89 returns the result in 1,9 (SP) resp. 2,7 (DP) seconds.

The above code is a first quick and dirty version, so take care. But it seems to run almost as fast as your version while it's much more compact and needs no separate handling of small or large arguments.

Dieter