Problems with the Gamma function



Dieter Kaiser wrote:
> I have finished the numerical routines for the Incomplete Gamma function and
> added some code for argument simplification. Especially, I have implemented some
> code to expand expressions like gamma_incomplete(n,z), gamma_incomplete(a+n,z)
> or gamma_incomplete(n+1/2,z) for n an integer value.

Cool!  Is there a limit on the value of n?  Do I get a huge expression
if, say, n = 100?

> (%i71) gamma(-500.0);
> (%o71) 
> Maxima encountered a Lisp error:
> 
>  Error in MSIZE-PREFIX [or a callee]: Bind stack overflow.

This must be a gcl problem.  Cmucl signals a floating-point overflow.

> 
> But for some values we get also a zero as result:
> 
> (%i75) gamma(-110.0);
> (%o75) 0.0

Cmucl returns - 2.52133190461148e-164.

> FIRST:
> I would like to suggest to implement code to test for Float and Bigfloat numbers
> representing a negative integer. In this case we should give a Maxima Error as
> we do for negative integers.

Yes.  Shouldn't be too hard in gamma-lanczos to check to see if we have
a negative integer and signal an error.

> 
> SECOND:
> As mentioned in the bug report SF[2013650] the problems with the overflow of
> values for float numbers arise from the routine gamma-lanczos. Perhaps we could
> implement code to detect the overflow and return an appropriate error.

Yes.  Could just wrap gamma-lanczos inside ignore-errors and if NIL is
returned, signal an error instead.  But it seems that gcl doesn't work
as desired:

(%i1) :lisp (ignore-errors (* 1d308 1d308))

Maxima encountered a Lisp error:



Ray