-----Urspr?ngliche Nachricht-----
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com]
> 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?
I have implemented no limit, but I think for big integers the results are not
very useful. Here are some examples for small values:
(%i17) gamma_incomplete(1/2,z);
(%o17) sqrt(%pi)*erfc(sqrt(z))
(%i18) gamma_incomplete(-1/2,z);
(%o18) 2*%e^-z/sqrt(z)-2*sqrt(%pi)*erfc(sqrt(z))
(%i20) gamma_incomplete(5/2,z);
(%o20) (3*z/4+3/2)*sqrt(z)*%e^-z+erfc(sqrt(z))
(%i21) gamma_incomplete(-5/2,z);
(%o21) -(-8*z^2/15+4*z/15-2/5)*%e^-z/z^(5/2)-8*sqrt(%pi)*erfc(sqrt(z))/15
(%i26) gamma_incomplete(a+1,z);
(%o26) z^a*%e^-z+a*gamma_incomplete(a,z)
(%i27) gamma_incomplete(a-1,z);
(%o27) -z^(a-1)*%e^-z/(a-1)-gamma_incomplete(a,z)/(1-a)
>> 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.
I would like to suggest to do the test in the simplifying routine simpgamma
before we call gamma-lanczos or the routines bffac or cbbfac. It seems to be
better not to call a routine for values which are not in the domain of the
function.
>> 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:
The overflow occurs within a call to exp. GCL does not generate an error in case
of an overflow, but returns a value which can be tested to be a float but is not
printable. But it is possible to test it against most-positive-double-float.
I have added a diff with the following changes:
1. Test for Float representation of zero and negative integers
2. Test for Bigfloat representation of zero and negative integers
3. Adding call to cbbfac for Complex Bigfloats
3. Test for overflow in gamma-lanczos
I have further modified slightly the calculation within gamma-lanczos to
avoid an overflow during at different places in the code.
The autoloading mechanism for cbbfac is not installed.
I have tested the code with GCl 2.6.8 and CLISP 2.46. The test_suite has almost
no problems. One test in expintegral_mac (Problem 40) has a slightly different
accuracy and fails.
Index: csimp2.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/csimp2.lisp,v
retrieving revision 1.21
diff -u -r1.21 csimp2.lisp
--- csimp2.lisp 14 Feb 2008 01:31:37 -0000 1.21
+++ csimp2.lisp 9 Sep 2008 18:21:29 -0000
@@ -179,13 +179,29 @@
(let* ((j (simpcheck (cadr x) z))
(jr ($realpart j))
(ji ($imagpart j)))
- (cond ((floatp j) (gammafloat j))
+ (cond ((and (floatp j)
+ (or (zerop j)
+ (and (< j 0)
+ (zerop (nth-value 1 (truncate j))))))
+ (merror "gamma(~:M) is undefined." j))
+ ((floatp j) (gammafloat j))
+ ((and ($bfloatp j)
+ (or (zerop1 j)
+ (and (eq ($sign j) '$neg)
+ (zerop1 (sub j ($truncate j))))))
+ (merror "gamma(~:M) is undefined." j))
(($bfloatp j) (mfuncall '$bffac (m+ j -1) $fpprec))
((and (numberp jr)
(numberp ji)
(or $numer (floatp jr) (floatp ji)))
(complexify (gamma-lanczos (complex (float jr)
(float ji)))))
+ ((and (mnump jr)
+ (mnump ji)
+ (or $numer ($bfloatp jr) ($bfloatp ji)))
+ (mfuncall '$cbffac
+ (add -1 ($bfloat jr) (mul '$%i ($bfloat ji)))
+ $fpprec))
((or (not (mnump j))
(ratgreaterp (simpabs (list '(%abs) j) 1 t) $gammalim))
(eqtest (list '(%gamma) j) x))
@@ -285,14 +301,16 @@
;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z)
;;
;; If z is a negative integer, Gamma(z) is infinity. Should
- ;; we test for this? Throw an error? What
+ ;; we test for this? Throw an error?
+ ;; The test will be done by the calling routine.
(/ (float pi)
(* (- z) (sin (* (float pi) z))
(gamma-lanczos (- z))))
(let* ((z (- z 1))
(zh (+ z 1/2))
(zgh (+ zh 607/128))
- (zp (expt zgh (/ zh 2)))
+ ;; Calculate log(zp) to avoid overflow at this point.
+ (lnzp (* (/ zh 2) (log zgh)))
(ss
(do ((sum 0.0)
(pp (1- (length c)) (1- pp)))
@@ -301,7 +319,19 @@
(incf sum (/ (aref c pp) (+ z pp))))))
(* (sqrt (float (* 2 pi)))
(+ ss (aref c 0))
- (* zp (exp (- zgh)) zp))))))
+ ;; Check for an overflow.
+ (let ((result (ignore-errors (exp (+ (- zgh) (* 2 lnzp))))))
+ (cond
+ (result
+ ;; ignore-errors don't work for GCL. GCL returns in case of
+ ;; an overflow a value which can be tested against the
+ ;; constant most-positive-double-float. We only test the
+ ;; realpart.
+ #+gcl
+ (if (> (realpart result) (/ most-positive-double-float 2))
+ (merror "Overflow in `gamma-lanczos'."))
+ result)
+ (t (merror "Overflow in `gamma-lanczos'.")))))))))
(defun gammafloat (a)
(realpart (gamma-lanczos (complex a 0.0))))
Success, CVS operation completed
Dieter Kaiser