Bug [1954846] bessel_i(1/2,0) -> divide by zero error
Subject: Bug [1954846] bessel_i(1/2,0) -> divide by zero error
From: Dieter Kaiser
Date: Mon, 7 Jul 2008 21:43:14 +0200
Below you can see a diff between the actual code and the code to eliminate the
usage of inifities and to test for a zero argument before expanding the half
integral Bessel functions. The results of this code I have presented earlier in
this thread.
Should I commit the changes and close the Bug report?
Index: bessel.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/bessel.lisp,v
retrieving revision 1.58
diff -u -r1.58 bessel.lisp
--- bessel.lisp 30 Apr 2008 13:50:11 -0000 1.58
+++ bessel.lisp 5 Jul 2008 10:30:47 -0000
@@ -892,12 +892,9 @@
((and (numberp order) (integerp order))
0)
((> ($realpart order) 0)
- 0)
- ((< ($realpart order) 0)
- ;; complex infinity
- '$infinity)
+ 0)
(t
- ;; the function is indeterminate
+ ;; in all other cases
(domain-error arg 'bessel_j))))
((bessel-numerical-eval-p order arg)
@@ -922,11 +919,16 @@
`((mtimes simp) -1 ((%bessel_j simp) ,(- order) ,arg))))
((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
- ;; When order is a fraction with a denominator of 2, we
- ;; can express the result in terms of elementary
- ;; functions.
- ;;
- (bessel-j-half-order rat-order arg))
+ (cond
+ ((and (numberp arg) (= arg 0))
+ ;; We dont't expand for a zero argument.
+ (if (> rat-order 0) 0 (domain-error arg 'bessel_j)))
+ (t
+ ;; When order is a fraction with a denominator of 2, we
+ ;; can express the result in terms of elementary
+ ;; functions.
+ ;;
+ (bessel-j-half-order rat-order arg))))
(t
(eqtest (list '(%bessel_j) order arg)
@@ -974,14 +976,7 @@
(rat-order nil))
(cond ((and (numberp arg) (= arg 0) (complex-number-p order))
- ;; We handle the different case for zero arg carefully.
- (cond ((and (numberp order) (zerop order))
- '$minf)
- ((not (= ($realpart order) 0))
- '$infinity)
- (t
- ;; the function is indeterminate
- (domain-error arg 'bessel_y))))
+ (domain-error arg 'bessel_y))
((bessel-numerical-eval-p order arg)
;; We have numeric order and arg and $numer is true, or
@@ -1003,14 +998,18 @@
`((mtimes simp) -1 ((%bessel_y simp) ,(- order) ,arg))))
((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
- ;; When order is a fraction with a denominator of 2, we
- ;; can express the result in terms of elementary
- ;; functions.
- ;;
- ;; Y[1/2](z) = -J[1/2](z) is a function of sin.
- ;; Y[-1/2](z) = -J[-1/2](z) is a function of cos.
-
- (bessel-y-half-order rat-order arg))
+ (cond
+ ((and (numberp arg) (= arg 0))
+ ;; We don't expand for a zero argument.
+ (domain-error arg 'bessel_y))
+ (t
+ ;; When order is a fraction with a denominator of 2, we
+ ;; can express the result in terms of elementary
+ ;; functions.
+ ;;
+ ;; Y[1/2](z) = -J[1/2](z) is a function of sin.
+ ;; Y[-1/2](z) = -J[-1/2](z) is a function of cos.
+ (bessel-y-half-order rat-order arg))))
(t
(eqtest (list '(%bessel_y) order arg)
@@ -1059,11 +1058,9 @@
(cond ((= order 0)
1)
((or (and (numberp order) (> order 0)) (integerp order))
- 0)
- ((and (numberp order) (< order 0))
- '$infinity)
+ 0)
(t
- ;; the function is indeterminate
+ ;; in all other cases domain-error
(domain-error arg 'bessel_i))))
((bessel-numerical-eval-p order arg)
@@ -1081,13 +1078,18 @@
(list '(%bessel_i) (- order) arg))
((and $besselexpand (setq rat-order (max-numeric-ratio-p order 2)))
- ;; When order is a fraction with a denominator of 2, we
- ;; can express the result in terms of elementary
- ;; functions.
- ;;
- ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
- ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
- (bessel-i-half-order rat-order arg))
+ (cond
+ ((and (numberp arg) (= arg 0))
+ ;; We don't expand for a zero argument.
+ (if (> rat-order 0) 0 (domain-error arg 'bessel_i)))
+ ( t
+ ;; When order is a fraction with a denominator of 2, we
+ ;; can express the result in terms of elementary
+ ;; functions.
+ ;;
+ ;; I[1/2](z) = sqrt(2/%pi/z)*sinh(z)
+ ;; I[-1/2](z) = sqrt(2/%pi/z)*cosh(z)
+ (bessel-i-half-order rat-order arg))))
(t
(eqtest (list '(%bessel_i) order arg)
@@ -1134,14 +1136,8 @@
(rat-order nil))
(cond ((and (numberp arg) (= arg 0) (complex-number-p order))
- ;; We handle the different case for zero arg carefully.
- (cond ((and (numberp order) (zerop order))
- '$inf)
- ((not (= ($realpart order) 0))
- '$infinity)
- (t
- ;; the function is indeterminate
- (domain-error arg 'bessel_k))))
+ ;; domain-error for all cases of zero arg
+ (domain-error arg 'bessel_k))
((bessel-numerical-eval-p order arg)
;; A&S 9.6.6
;; K[-v](x) = K[v](x)
@@ -1159,12 +1155,17 @@
((and $besselexpand
(setq rat-order (max-numeric-ratio-p order 2)))
- ;; When order is a fraction with a denominator of 2, we
- ;; can express the result in terms of elementary
- ;; functions.
- ;;
- ;; K[1/2](z) = sqrt(2/%pi/z)*exp(-z) = K[1/2](z)
- (bessel-k-half-order rat-order arg))
+ (cond
+ ((and (numberp arg) (= arg 0))
+ ;; We don't expand for a zero argument.
+ (domain-error arg 'bessel_k))
+ (t
+ ;; When order is a fraction with a denominator of 2, we
+ ;; can express the result in terms of elementary
+ ;; functions.
+ ;;
+ ;; K[1/2](z) = sqrt(2/%pi/z)*exp(-z) = K[1/2](z)
+ (bessel-k-half-order rat-order arg))))
(t
(eqtest (list '(%bessel_k) order arg)
exp)))))
CVS-Vorgang erfolgreich abgeschlossen
Dieter Kaiser