Well, I'll leave Lambert for the real mathematicians out there but here's a pragmatic approach for the LnGamma case.-
The first program is the "brute force" approach, adding terms to the summ until their contribution isn't relevant to the (rounded) result. VERY slow, and not an appropriate method- just good for comparison purposes.
LnG(z) = -Gz - Lnz + SUM [(z/k)-Ln(1+z/k) |k=1,2.... ]
No doubt you'll recognized Stirling's approach embedded in the second method. The routine calculates LnG(z+7) for enhanced accuracy (about E-9), and adjusts it back appropriately.
MCODE functions from the 41Z module help bring down the execution time to manageable on the real 41, and of course irrelevant on emulators.
2LnG(z)=Ln(2pi)-Ln(z)+z{2Lnz + Ln[z*sinh(1/z)+(1/810*z^6)]-2}
LnG(z) = LnG(z+n) - Ln [ PROD(z+k)|k=1,2..(n-1)]
Method-2 is much faster, more accurate and takes less program
memory... can you ask for more? Yes, an MCODE version will be implemented on the 41Z during the next few weeks :)
Enjoy,
ÁM.
LBL "ZLNG" LBL "ZLNG2"
1 7
STO 02 +
RDN ZST0 (00)
ZSTO (00) 6
XEQ 05 CHS
LBL 00 Z^X
ZENTER^ 810
XEQ 05 ST/ Z
Z+ /
Z=WR? ZRCL (00)
GTO 02 ZINV
GTO 00 ZSINH
LBL 02 ZRCL (00)
ZRCL (00) Z*
ZLN Z+
Z- ZLN
ZRCL (00) ZRCL (00)
0,5772156649 ZLN
ST* Z ZDBL
* Z+
Z- 2
ZAVIEW -
RTN ZRCL (00)
LBL 05 Z*
ZRCL (00) ZRCL (00)
RCL 02 ZLN
ST/ Z Z-
/ PI
ZENTER^ ST+ X
1 LN
+ +
ZLN ZHALF
Z- ZRCL (00)
1 7
ST+ 02 -
RDN ZGPROD
END ZLN
Z-
ZAVIEW
END
Edited: 22 Sept 2009, 9:16 a.m.