Faster inverse gamma and factorial for the WP 34S

[Dieter]

(02-16-2015 05:46 PM)Bit Wrote:  The built-in digamma function that Pauli mentioned, and which isn't available by default, is more accurate. Perhaps up to 33 digits in double precision mode if I see it correctly.

That's what should be expected with the internal 39-digit precision of C-coded functions. ;-)

Back to the current solution: The Ψ code in the user code library was designed for 10-digit accuracy on a 41 series calculator. It uses the asymptotic method given in Abramowitz & Stegun (1972) 6.3.18 for x > 8 and four terms of the series. Smaller arguments are evaluated via Ψ(x+1) = Ψ(x)+1/x.

For the 34s I suggest using x > 16 (in SP) and six terms, which – when evaluated with sufficient precision – should provide an absolute error order of 10^–18. The results are even better with a slight tweak: instead of B₁₂ = –691/2730 I simply use –1/4. ;-)

Here is some experimental code that should work for positive arguments. It is not yet complete and cannot replace the current Ψ routine, but I think it points into the right direction:

Code:

LBL"PSI"
#016
DBL?
x²
STO 01
STO-01
x<> Y
1/x
STO+01
x<> L
INC X
x<? Y
BACK 005
STO 00
1/x
x²
FILL
#004
+/-
/
#011
1/x
+
x
#020
1/x
-
x
#021
1/x
+
x
#010
1/x
-
x
INC X
x
#012
/
RCL 00
STO+X
1/x
+
RCL 00
LN
RCL-01
x<> Y
-
RTN

Here are some results:

Code:

SP mode:
  1            XEQ"PSI"  -0,5772156649015329   exact
  2            XEQ"PSI"   0,4227843350984671   exact
 20            XEQ"PSI"   2,970523992242149    exact
1,46163211449  XEQ"PSI"  -2,949306481 E-8      7 digits correct (abs. err ~ 1E-15)

DP mode:
  1            XEQ"PSI"  -0,5772156649015328606065120900824037   32 digits correct
  2            XEQ"PSI"   0,4227843350984671393934879099175963   32 digits correct
 20            XEQ"PSI"   2,970523992242149050877256978825982    33 digits correct
1,46163211449  XEQ"PSI"  -2,9493065735632104820662566586 E-8     25 digits correct (abs. err ~ 3E-33)

(02-16-2015 05:46 PM)Bit Wrote:  If you don't have a build environment but would like to play with the built-in function, I can send you some binaries that have it enabled.

Thank you very much. But let me try first what can be done with user code. The mentioned results might give a first impression, now let's see what happens if the results move very close to zero, i.e. near the Gamma minimum at 1,4616...

I wonder what the internal digamma function might return for Ψ(1,4616321449768362) and ...363? Can you post the results with all 34 digits?

Dieter