From 191ca289d3c26813e0244251a983a62f2057347a Mon Sep 17 00:00:00 2001
From: Ralf Stephan <ralf@ark.in-berlin.de>
Date: Wed, 28 Jan 2015 09:18:47 +0100
Subject: [PATCH] 17678: special values of Bessel functions

---
 src/sage/functions/bessel.py | 101 +++++++++++++++++++++++++++++------
 1 file changed, 84 insertions(+), 17 deletions(-)

diff --git a/src/sage/functions/bessel.py b/src/sage/functions/bessel.py
index 8df019aaa68..d6cd9cae026 100644
--- a/src/sage/functions/bessel.py
+++ b/src/sage/functions/bessel.py
@@ -97,7 +97,7 @@
         sage: bessel_J(0, x)
         bessel_J(0, x)
         sage: bessel_J(0, 0)
-        bessel_J(0, 0)
+        1
         sage: bessel_J(0, x).diff(x)
         -1/2*bessel_J(1, x) + 1/2*bessel_J(-1, x)
 
@@ -194,11 +194,12 @@
 
 # remove after deprecation period
 from sage.calculus.calculus import maxima
-from sage.functions.other import real, imag
+from sage.functions.trig import sin, cos
+from sage.functions.other import real, imag, sqrt
 from sage.misc.sage_eval import sage_eval
 from sage.rings.real_mpfr import RealField
 from sage.plot.plot import plot
-from sage.rings.all import ZZ
+from sage.rings.all import ZZ, QQ
 
 
 class Function_Bessel_J(BuiltinFunction):
@@ -303,6 +304,34 @@ def __init__(self):
                                                   maxima='bessel_j',
                                                   sympy='besselj'))
 
+    def _eval_(self, n, x):
+        """
+        EXAMPLES::
+
+            sage: bessel_J(0, 0)
+            1
+            sage: bessel_J(5/2, 0)
+            0
+            sage: bessel_J(-5/2, 0)
+            +Infinity
+            sage: bessel_J(1/2, x)
+            sqrt(2)*sqrt(1/(pi*x))*sin(x)
+            sage: bessel_J(-1/2, x)
+            sqrt(2)*sqrt(1/(pi*x))*cos(x)
+        """
+        from sage.rings.infinity import infinity
+        if x == 0:
+            if n == 0:
+                return ZZ(1)
+            elif n > 0 or n in ZZ:
+                return ZZ(0)
+            else:
+                return infinity
+        if n == QQ(1)/2:
+            return sqrt(2/pi/x) * sin(x)
+        elif n == QQ(-1)/2:
+            return sqrt(2/pi/x) * cos(x)
+
     def _evalf_(self, n, x, parent=None, algorithm=None):
         """
         EXAMPLES::
@@ -475,6 +504,25 @@ def __init__(self):
                                                   maxima='bessel_y',
                                                   sympy='bessely'))
 
+    def _eval_(self, n, x):
+        """
+        EXAMPLES::
+
+            sage: bessel_Y(1, 0)
+            -Infinity
+            sage: bessel_Y(1/2, x)
+            -sqrt(2)*sqrt(1/(pi*x))*cos(x)
+            sage: bessel_Y(-1/2, x)
+            sqrt(2)*sqrt(1/(pi*x))*sin(x)
+        """
+        from sage.rings.infinity import infinity
+        if x == 0 and n >= 0:
+            return -infinity
+        if n == QQ(1)/2:
+            return -sqrt(2/pi/x) * cos(x)
+        elif n == QQ(-1)/2:
+            return sqrt(2/pi/x) * sin(x)
+
     def _evalf_(self, n, x, parent=None, algorithm=None):
         """
         EXAMPLES::
@@ -645,25 +693,40 @@ def _eval_(self, n, x):
         EXAMPLES::
 
             sage: y=var('y')
-            sage: bessel_I(y,x)
+            sage: bessel_I(y, x)
             bessel_I(y, x)
-            sage: bessel_I(0.0, 1.0)
-            1.26606587775201
+            sage: bessel_I(0, 0)
+            1
+            sage: bessel_I(7/2, 0)
+            0
+            sage: bessel_I(-7/2, 0)
+            +Infinity
             sage: bessel_I(1/2, 1)
             sqrt(2)*sinh(1)/sqrt(pi)
             sage: bessel_I(-1/2, pi)
             sqrt(2)*cosh(pi)/pi
         """
+        from sage.rings.infinity import infinity
         # special identities
-        if n == Integer(1) / Integer(2):
+        if x == 0:
+            if n == 0:
+                return ZZ(1)
+            elif n > 0 or n in ZZ:
+                return ZZ(0)
+            else:
+                return infinity
+        if n == QQ(1)/2:
             return sqrt(2 / (pi * x)) * sinh(x)
-        elif n == -Integer(1) / Integer(2):
+        elif n == -QQ(1)/2:
             return sqrt(2 / (pi * x)) * cosh(x)
 
+
     def _evalf_(self, n, x, parent=None, algorithm=None):
         """
         EXAMPLES::
 
+            sage: bessel_I(0.0, 1.0)
+            1.26606587775201
             sage: bessel_I(1,3).n(digits=20)
             3.9533702174026093965
         """
@@ -739,9 +802,9 @@ class Function_Bessel_K(BuiltinFunction):
         -1/2*bessel_K(3, x) - 1/2*bessel_K(1, x)
 
         sage: bessel_K(1/2, x)
-        bessel_K(1/2, x)
+        sqrt(1/2)*sqrt(pi)*e^(-x)/sqrt(x)
         sage: bessel_K(1/2, -1)
-        bessel_K(1/2, -1)
+        -I*sqrt(1/2)*sqrt(pi)*e
         sage: bessel_K(1/2, 1)
         sqrt(1/2)*sqrt(pi)*e^(-1)
 
@@ -817,23 +880,27 @@ def _eval_(self, n, x):
         """
         EXAMPLES::
 
-            sage: bessel_K(1,0)
-            bessel_K(1, 0)
-            sage: bessel_K(1.0, 0.0)
-            +infinity
-            sage: bessel_K(-1, 1).n(128)
-            0.60190723019723457473754000153561733926
+            sage: bessel_K(1, 0)
+            +Infinity
+            sage: bessel_K(1/2, x)
+            sqrt(1/2)*sqrt(pi)*e^(-x)/sqrt(x)
         """
+        from sage.rings.infinity import infinity
         # special identity
-        if n == Integer(1) / Integer(2) and x > 0:
+        if x == 0:
+            return infinity
+        if n == QQ(1)/2 or n == -QQ(1)/2 and x > 0:
             return sqrt(pi / 2) * exp(-x) * x ** (-Integer(1) / Integer(2))
 
+
     def _evalf_(self, n, x, parent=None, algorithm=None):
         """
         EXAMPLES::
 
             sage: bessel_K(0.0, 1.0)
             0.421024438240708
+            sage: bessel_K(-1, 1).n(128)
+            0.60190723019723457473754000153561733926
             sage: bessel_K(0, RealField(128)(1))
             0.42102443824070833333562737921260903614
         """