Complex Gamma Function... revisited Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 09-14-2009, 02:54 PM Valid for the complex plane plane where Gamma is defined (i.e. includes Re(z)<0, and excludes negative integers and zero). Here's the program listing, a nice example of the HP-41Z function set (can you tell it's FOCAL?? :-) Uses the Lanczos approximation, split in three parts as per the function names clearly show. More soon... ```LBL "ZG" CF 00 X<0? SF 00 X<0? ZNEG ZREPL ZGSUM LASTZ ZGLNC Z* Z<>W ZGPROD Z/ FC?C 00 GTO 00 ZINV Z<>W ZGNZG Z* LBL 00 ZAVIEW END ``` hugh steers Senior Member Posts: 536 Threads: 56 Joined: Jul 2005 09-14-2009, 05:03 PM I have had a lot of trouble with the Lanczos approximation. The coefficients of the series vary wildly and cause large swings in value. This causes a lot of cancellation which destroys accuracy. Currently, Im using increased precision to cope with this problem. I'd prefer a more stable formula than Lanczos. There's another method by Spouge, but you need a lot more series terms for convergence. it will be interesting to see how many digits you get on the 41. Marcus von Cube, Germany Posting Freak Posts: 3,283 Threads: 104 Joined: Jul 2005 09-15-2009, 01:40 PM Go and ask Viktor T. Toth at rskey.org! He's dedicated part of his life to Gamma. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 09-15-2009, 02:22 PM I think Toth calculates the Gamma function for real arguments only. This way he can apply the Gamma (or Ln Gamma) to a whole collection of programmable calculators. Egan Ford Posting Freak Posts: 1,619 Threads: 147 Joined: May 2006 09-15-2009, 04:06 PM Nice. How about a complex LogGamma and a complex Lambert W? hugh steers Senior Member Posts: 536 Threads: 56 Joined: Jul 2005 09-15-2009, 04:15 PM You kind-of get LnGamma anyway. Out of Lanczos you wind up with three numbers; S, A and B (say). Gamma(z) = S*exp(A*ln(B)-B) so using the same subroutine, you can have, LnGamma(z) = A*ln(B) - B + Ln(S) However, the latter might not be the best for larger numbers. Marcus von Cube, Germany Posting Freak Posts: 3,283 Threads: 104 Joined: Jul 2005 09-15-2009, 04:29 PM Namir, that depends on the calculator. On machines with decent complex support (like the TI-92) he calculates the complex Gamma function. His general article about Gamma is worth reading. Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 09-15-2009, 04:36 PM Here's where I'm going to show my utter lack of math finesse but what the heck: LBL "LNZG" ZG ZLN END Of course this is just a teaser, will look into a proper way during the week-end :) Cheers, 'AM Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 09-15-2009, 04:50 PM That exactly is the source I have used to programm it. http://www.rskey.org/gamma.htm To the careful reader it's clear that: ZGSUM computes the sum of the seven terms (qk * Z^k) ZGPROD is thhe product of (z+k) ZGLNC calculates the trascendent part of the formula, and ZGNZG corrects the expression for Re(z)<0 Cheers, Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 09-16-2009, 12:31 AM Quote: Nice. How about a complex LogGamma and a complex Lambert W? Wikipedia has an algorithm for Lambert's W function that seems okay and converges in the complex plane: here. This is the one I'm using in the 20b scientific firmware. Start with the approximation and then iterate to the solution - I've limited this to 20 times through. For complex (and real) LogGamma, I'm using Lanczos. Look for An Analysis of the Lanczos Gamma Approximation by Glendon Ralph Pugh pages 134 and 135 for the formula and coefficients of the approximation I'm using. - Pauli Les Wright Posting Freak Posts: 1,368 Threads: 212 Joined: Dec 2006 09-16-2009, 02:27 AM Hugh, I have done a lot of fiddling with Spouge, and the cancellation and digit loss issues persist there--only worse. You need a LOT of guard digits to preserve accuracy. I understand you double precision BCD-20 type has gone a long way to help here. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 09-16-2009, 10:31 AM Marcus, You are correct! I did revisit Toth's Gamma function math page (have not seen it for a while) and saw the Lancsoz approximation that is being talked about here. While series approximation are in general easy to work with, a few here and there are tough to tame! Using the old programmable calculators, series approximation were really time consuming. With today's machine's we have faster CPU and the approximation (again in general) are not too bad. Namir Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 09-22-2009, 08:51 AM 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. Steven Thomas Smith Junior Member Posts: 2 Threads: 0 Joined: Sep 2009 09-25-2009, 12:16 PM See GSL's complex Gamma implementation -- they use Lanczos for the right half-plane, Stirling for the left away from the negative real axis, and reflection to avoid the unstable parts -- actually, they use reflection to avoid Stirling altogether. It's a good implementation. Edited: 25 Sept 2009, 12:19 p.m. « Next Oldest | Next Newest »

 Possibly Related Threads... Thread Author Replies Views Last Post HP50g: Writing a function that returns a function Chris de Castro 2 949 12-10-2013, 06:49 PM Last Post: Han HP Prime: complex numbers in CAS. Alberto Candel 1 794 12-06-2013, 02:36 PM Last Post: parisse [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 1,584 12-05-2013, 07:40 AM Last Post: cyrille de Brébisson 17BII & 17BII+ Discounted Payback Period Revisited Tom Neudorfl 8 1,473 11-25-2013, 10:28 AM Last Post: Don Shepherd Touch periodic table on HP Prime - revisited Terje Vallestad 2 667 11-23-2013, 11:47 AM Last Post: Mic Complex Number Entry on Prime Jeff O. 19 2,652 11-16-2013, 12:34 PM Last Post: Jeff O. HP Prime complex results Javier Goizueta 0 501 10-06-2013, 12:59 PM Last Post: Javier Goizueta HP15c continued fraction for Ln(Gamma) Tom Grydeland 0 476 09-30-2013, 05:48 AM Last Post: Tom Grydeland HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 2,053 09-30-2013, 03:44 AM Last Post: From Hong Kong [HP-Prime xcas] operations with complex numbers + BUGs + Request CompSystems 9 1,605 09-08-2013, 10:40 PM Last Post: CompSystems

Forum Jump: