Problems with gamma, lgamma, and alternatives

[Michael Creel]

In a message from last October, Tim Reluga wrote: 

>       Having recently had need for accurate evaluation of the gamma 
> function on the complex unit circle, I was disappointed to discover 
> Octave's current implementation doesn't handle complex arguements. 
>       I'm not sure how the current implementation works.  It sounds like 
> it may be legacy fortran code.  As such, there was no easy resolution to 
> my problem there. 
>       In the archives, I've found some old implementations based on 
> Stirling's asymptotic expansion.  Unfortunately, this is a divergent 
> series, and unable to provide arbitrary accuracy. 
>       Looking for alternatives, I've found Viktor Toth's web page 
> <http://www.rskey.org/gamma.htm> where he discusses using a Lanczos 
> Approximation in the complex domain.  Even better, Paul Godfrey has 
> explored the accuracy of the Lanczos approximation and already implemented 
> it in an octave-compatible fashion. See 
> <http://winnie.fit.edu/~gabdo/paulbio.html>.  I've emailed Mr. Godfrey, 
> and he has agreed to allow his implementation to be included with Octave, 
> in the maintainers are so inclined. If there is anything I can do to 
> facilitate it's incorporation into the next development release, let me 
> know. 
> 
>                                               Tim Reluga

I have been having some problems with the current implementation of gamma and 
lgamma in Octave - specifically crashes of the form

octave:1> gamma(1000)
 ***MESSAGE FROM ROUTINE DGAMMA IN LIBRARY SLATEC.
 ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  X SO BIG GAMMA OVERFLOWS
 *  ERROR NUMBER = 3
 *
 ***END OF MESSAGE

 ***JOB ABORT DUE TO FATAL ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     DGAMMA     X SO BIG GAMMA OVERF         3         2         1

These usually occur during intermediate iterations of a minimization routine, 
where it's difficult to control the value that will be passed to gamma. The 
problem is that everything crashes gracelessly here, with no hope for 
recovery. On the other hand, I have tried the routines by Paul Godfrey 
mentioned in the quote above. These do not cause crashes, but are much slower 
than the SLATEC versions, which I take it are Fortran or C. So, my question 
is, is there any way to get crash control and speed at the same time? I 
suppose writing the Godfrey routines in Fortran or C and linking them in 
would work. However, this is over my head at the moment - I don't do C or 
Fortran. Maybe the solution is to learn a little. Any advice would be 
appreciated. Thanks, Michael



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------