Problems with the Gamma function



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.

One open thing is to have a more carefully look at specific values. So I had a
look at the Gamma function too, because a correct implementation of the Gamma
function is important for the other Gamma functions.

The gamma function has for negative integers including zero the result infinity.
It is a pitty that we can not return infinities as a value because of the
missing arithmetic for infinities. So the negative integers are implemented as a
Maxima Error.

My problem is what we should expect and implement for a Float or Bigfloat number
representing an integer or a zero value. I have found the following results:

/* Integer 0 -> Maxima error */

(%i5) gamma(0);
gamma(0) is undefined
 -- an error.  To debug this try debugmode(true);

/* Float 0.0 -> Lisp error */

(%i6) gamma(0.0);
Maxima encountered a Lisp error:

 Error in MACSYMA-TOP-LEVEL [or a callee]: Zero divisor.

Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

/* Bigfloat 0.0b0 -> positive big result */

(%i7) gamma(0.0b0);
(%o7) 2.745565351356511b17

/* Here are the results for a negative integer and
   a float and Bigfloat representing an integer */

(%i8) gamma(-1);
gamma(-1) is undefined
 -- an error.  To debug this try debugmode(true);

(%i9) gamma(-1.0);
(%o9) -2.5653898036572368E+16

(%i10)gamma(-1.0b0);
(%o10) -1.372779803348707b17

(%i71) gamma(-500.0);
(%o71) 
Maxima encountered a Lisp error:

 Error in MSIZE-PREFIX [or a callee]: Bind stack overflow.

Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

But for some values we get also a zero as result:

(%i75) gamma(-110.0);
(%o75) 0.0

For big values we get further Lisp errors (this is reported as a bug
SF[2013650]):

(%i70) gamma(250.0);
(%o70) 
Maxima encountered a Lisp error:

 Error in PROGN [or a callee]: Can't print a non-number.

Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

So we get different results for negative integers and have problems when the
result overflows.

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.

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.

THIRD:
Perhaps it is interesting to implement code to allow further simplification for
the Gamma function. Depending on a flag e.g. $gammaexpand we could do an
expansion for expression like gamma(n*z) or gamma(n+z) where n is an integer.

Dieter Kaiser