[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Problems with gamma, lgamma, and alternatives
From: |
Michael Creel |
Subject: |
Problems with gamma, lgamma, and alternatives |
Date: |
Tue, 25 Mar 2003 13:19:20 -0800 |
User-agent: |
KMail/1.5.1 |
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
-------------------------------------------------------------
- Problems with gamma, lgamma, and alternatives,
Michael Creel <=