Ticket #12455: trac_12455_newstyle_airy.patch

new file sage/functions/airy.py

@@ +1 @@
+#*****************************************************************************
+#       Copyright (C) 2010 Oscar Gerardo Lazo Arjona algebraicamente@gmail.com
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#    General Public License for more details.
+#
+#  The full text of the GPL is available at:
+#
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+from sage.symbolic.function import BuiltinFunction
+from sage.symbolic.expression import Expression
+from sage.symbolic.function import is_inexact
+from sage.structure.coerce import parent as sage_structure_coerce_parent
+
+from sage.functions.other import gamma
+from sage.rings.integer_ring import ZZ
+from sage.rings.real_double import RDF
+from sage.functions.special import meval
+from sage.calculus.var import var
+from sage.calculus.functional import derivative
+
+# from sage.rings.real_mpfr import RR
+from sage.rings.rational import Rational as R
+
+class FunctionAiryAiGeneral(BuiltinFunction):
+    def __init__(self):
+        r"""
+        The generalized derivative of the Airy Ai function
+
+        INPUT:
+        
+        - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to `z`. 
+          For `\alpha = n = 1,2,3,\ldots` this gives the derivative `\operatorname{Ai}^{(n)}(z)`, 
+          and for `\alpha = -n = -1,-2,-3,\ldots` this gives the `n`-fold iterated integral
+         
+        .. math ::
+         
+            f_0(z) = \operatorname{Ai}(z)
+         
+            f_n(z) = \int_0^z f_{n-1}(t) dt.
+
+        - ``x``.- The argument of the function.
+
+        - ``hold``.- Whether or not to stop from returning higher derivatives in terms of 
+          `\operatorname{Ai}(x)` and `\operatorname{Ai}'(x)`.
+        
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: airy_ai(-2,x)
+            airy_ai(-2, x)
+            sage: derivative(airy_ai(-2,x),x)
+            airy_ai(-1, x)
+            sage: airy_ai(n,x)
+            airy_ai(n, x)
+            sage: derivative(airy_ai(n,x),x)
+            airy_ai(n + 1, x)
+            sage: airy_ai(2,x,True)
+            airy_ai(2, x)
+            sage: derivative(airy_ai(2,x,hold_derivative=True),x)
+            airy_ai(3, x)
+        """
+
+        BuiltinFunction.__init__(self, "airy_ai", nargs=2,
+                latex_name=r"\operatorname{Ai}")
+
+    def _derivative_(self, alpha, *args, **kwds):
+        """
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: derivative(airy_ai(n,x),x) # indirect doctest
+            airy_ai(n + 1, x)
+        """
+
+        x=args[0]
+        return airy_ai_general(alpha+1,x)
+
+    def _eval_(self, alpha, *args):
+        """
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: airy_ai(-2,1.0) # indirect doctest
+            0.136645379421096
+            sage: airy_ai(n,1.0) # indirect doctest
+            airy_ai(n, 1.00000000000000)
+        """
+
+        x=args[0]
+        if not isinstance(x,Expression) and not isinstance(alpha,Expression):
+            if is_inexact(x):
+                return self._evalf_(alpha,x)
+        else:
+            return None
+
+    def _evalf_(self, alpha, x, **kwargs):
+        """
+        EXAMPLES::
+
+            sage: airy_ai(-2,1.0) # indirect doctest
+            0.136645379421096            
+        """
+        algorithm = kwargs.get('algorithm', None) or 'mpmath'
+        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
+        prec = parent.prec() if hasattr(parent, 'prec') else 53
+        if algorithm == 'mpmath':
+            import mpmath
+            from sage.libs.mpmath import utils as mpmath_utils
+            return mpmath_utils.call(mpmath.airyai, x, derivative=alpha, parent=parent)
+        elif algorithm == 'maxima':
+            raise NotImplementedError("general case not available in maxima")
+        else:
+            raise ValueError("unknown algorithm")
+
+class FunctionAiryAiSimple(BuiltinFunction):
+    def __init__(self):
+        """
+        The class for the Airy Ai function.
+
+        EXAMPLES::
+
+            sage: from sage.functions.airy import FunctionAiryAiSimple
+            sage: airy_ai_simple= FunctionAiryAiSimple()
+            sage: f=airy_ai_simple(x); f
+            airy_ai(x)
+        """
+
+        BuiltinFunction.__init__(self, "airy_ai", latex_name=r'\operatorname{Ai}',
+                conversions=dict(mathematica='AiryAi', maxima='airy_ai'))
+
+    def _derivative_(self, x, diff_param=None):
+        """
+        EXAMPLES::
+        
+            sage: derivative(airy_ai(x),x) # indirect doctest
+            airy_ai_prime(x)
+        """
+        return airy_ai_prime(x)
+
+    def _eval_(self, x):
+        """
+        EXAMPLES::
+        
+            sage: airy_ai(0) # indirect doctest
+            1/3*3^(1/3)/gamma(2/3)
+            sage: airy_ai(0.0) # indirect doctest
+            0.355028053887817
+            sage: airy_ai(I) # indirect doctest
+            airy_ai(I)
+            sage: airy_ai(1.0*I) # indirect doctest
+            0.331493305432141 - 0.317449858968444*I
+        """
+
+        if not isinstance(x,Expression):
+            if is_inexact(x):
+                return self._evalf_(x)
+            elif x==0:
+                r=ZZ(2)/3
+                return 1/(3**(r)*gamma(r))
+        else:
+            return None
+
+    def _evalf_(self, x, **kwargs):
+        """
+        EXAMPLES::
+        
+            sage: airy_ai(0.0) # indirect doctest
+            0.355028053887817
+            sage: airy_ai(1.0*I) # indirect doctest
+            0.331493305432141 - 0.317449858968444*I
+
+        We can use several methods for numerical evaluation::
+
+            sage: airy_ai(3).n(algorithm='maxima')
+            0.00659113935746
+            sage: airy_ai(3).n(algorithm='mpmath')
+            0.00659113935746072
+            sage: airy_ai(3).n(algorithm='mpmath',prec=100)
+            0.0065911393574607191442574484080
+            
+        TESTS::
+
+            sage: airy_ai(3).n(algorithm='maxima', prec=70)
+            Traceback (most recent call last):
+            ...
+            ValueError: for the maxima algorithm the precision must be 53
+        
+        """
+        algorithm = kwargs.get('algorithm', None) or 'mpmath'
+        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
+        prec = parent.prec() if hasattr(parent, 'prec') else 53
+        if algorithm == 'mpmath':
+            import mpmath
+            from sage.libs.mpmath import utils as mpmath_utils
+            return mpmath_utils.call(mpmath.airyai, x, parent=parent)
+        elif algorithm == 'maxima':
+            if prec != 53:
+                raise ValueError, "for the maxima algorithm the precision must be 53"
+            return RDF(meval("airy_ai(%s)" % RDF(x)))
+        else:
+            raise ValueError("unknown algorithm")
+
+class FunctionAiryAiPrime(BuiltinFunction):
+    def __init__(self):
+        """
+        The derivative of the Airy Ai function; see :func:`airy_ai`
+        for the full documentation.
+
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: airy_ai_prime(x)
+            airy_ai_prime(x)
+            sage: airy_ai_prime(0)
+            -1/3*3^(2/3)/gamma(1/3)
+        """
+
+        BuiltinFunction.__init__(self, "airy_ai_prime",
+                latex_name=r"\operatorname{Ai}'",
+                conversions=dict(mathematica='AiryAiPrime', maxima='airy_dai'))
+
+    def _derivative_(self, x, diff_param=None):
+        """
+        EXAMPLES::
+
+            sage: derivative(airy_ai_prime(x),x) # indirect doctest
+            x*airy_ai(x)
+        """
+        return x*airy_ai_simple(x)
+
+    def _eval_(self, x):
+        """
+        EXAMPLES::
+
+            sage: airy_ai(1,0) # indirect doctest
+            -1/3*3^(2/3)/gamma(1/3)
+            sage: airy_ai(1,0.0) # indirect doctest
+            -0.258819403792807
+        """
+        if not isinstance(x,Expression):
+            if is_inexact(x):
+                return self._evalf_(x)
+            elif x==0:
+                r=ZZ(1)/3
+                return -1/(3**(r)*gamma(r))
+        else:
+            return None
+
+    def _evalf_(self, x, **kwargs):
+        """
+        EXAMPLES::
+        
+            sage: airy_ai(1,0.0) # indirect doctest
+            -0.258819403792807
+
+        We can use several methods for numerical evaluation::
+
+            sage: airy_ai(1,4).n(algorithm='maxima')
+            -0.0019586409502
+            sage: airy_ai(1,4).n(algorithm='mpmath')
+            -0.00195864095020418
+            sage: airy_ai(1,4).n(algorithm='mpmath', prec=100)
+            -0.0019586409502041789001381409184
+
+        TESTS::
+
+            sage: airy_ai(1, 4).n(algorithm='maxima', prec=70)
+            Traceback (most recent call last):
+            ...
+            ValueError: for the maxima algorithm the precision must be 53
+
+
+        """
+        algorithm = kwargs.get('algorithm', None) or 'mpmath'
+        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
+        prec = parent.prec() if hasattr(parent, 'prec') else 53
+        if algorithm == 'mpmath':
+            import mpmath
+            from sage.libs.mpmath import utils as mpmath_utils
+            return mpmath_utils.call(mpmath.airyai, x, derivative=1, parent=parent)
+        elif algorithm == 'maxima':
+            if prec != 53:
+                raise ValueError, "for the maxima algorithm the precision must be 53"
+            return RDF(meval("airy_dai(%s)" % RDF(x)))
+        else:
+            raise ValueError("unknown algorithm")
+
+
+airy_ai_general=FunctionAiryAiGeneral()
+airy_ai_simple= FunctionAiryAiSimple()
+airy_ai_prime=  FunctionAiryAiPrime()
+
+
+def airy_ai(alpha,x=None, hold_derivative=False, *args, **kwds):
+    r""" 
+    The Airy Ai function `\operatorname{Ai}(x)` is one of the two
+    linearly independent solutions to the Airy differental equation
+    `f''(z) +f(z)x=0`, defined by the initial conditions:
+
+    .. math ::
+        \operatorname{Ai}(0)=\frac{1}{2^{2/3} \Gamma(\frac{2}{3})},
+
+        \operatorname{Ai}'(0)=-\frac{1}{2^{1/3} \Gamma(\frac{1}{3})}.
+
+    Another way to define the Airy Ai function is:
+
+    .. math::
+        \operatorname{Ai}(x)=\frac{1}{\pi}\int_0^\infty
+        \cos\left(\frac{1}{3}t^3+xt\right) dt.
+
+    INPUT:
+    
+    - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to `z`. 
+      For `\alpha = n = 1,2,3,\ldots` this gives the derivative `\operatorname{Ai}^{(n)}(z)`, 
+      and for `\alpha = -n = -1,-2,-3,\ldots` this gives the `n`-fold iterated integral
+     
+    .. math ::
+     
+        f_0(z) = \operatorname{Ai}(z)
+     
+        f_n(z) = \int_0^z f_{n-1}(t) dt.
+
+    - ``x``.- The argument of the function.
+
+    - ``hold_derivative``.- Whether or not to stop from returning higher derivatives 
+       in terms of `\operatorname{Ai}(x)` and `\operatorname{Ai}'(x)`.
+    
+
+    EXAMPLES::
+
+        sage: n,x=var('n x')
+        sage: airy_ai(x)
+        airy_ai(x)
+
+    It can return derivatives or integrals::
+    
+        sage: airy_ai(1,x)
+        airy_ai_prime(x)
+        sage: airy_ai(2,x)
+        x*airy_ai(x)
+        sage: airy_ai(2,x,True)
+        airy_ai(2, x)
+        sage: airy_ai(-2,x)
+        airy_ai(-2, x)
+        sage: airy_ai(n, x)
+        airy_ai(n, x)
+
+    It can be evaluated symbolically or numerically for real or complex values::
+
+        sage: airy_ai(0)
+        1/3*3^(1/3)/gamma(2/3)
+        sage: airy_ai(0.0)
+        0.355028053887817
+        sage: airy_ai(I)
+        airy_ai(I)
+        sage: airy_ai(1.0*I)
+        0.331493305432141 - 0.317449858968444*I
+
+    The functions can be evaluated numerically using either mpmath (the default)
+    or maxima, where mpmath can compute the values to arbitrary precision::
+
+        sage: airy_ai(2).n(prec=100)
+        0.034924130423274379135322080792
+        sage: airy_ai(2).n(algorithm='mpmath',prec=100)
+        0.034924130423274379135322080792
+        sage: airy_ai(2).n(algorithm='maxima')
+        0.0349241304233
+
+    And the derivatives can be evaluated::
+
+        sage: airy_ai(1,0)
+        -1/3*3^(2/3)/gamma(1/3)
+        sage: airy_ai(1,0.0)
+        -0.258819403792807
+
+    Plots::
+
+        sage: plot(airy_ai(x),(x,-10,5))+plot(airy_ai_prime(x),(x,-10,5),color='red')
+        
+    **References**
+     
+    - Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 10"
+
+    - http://en.wikipedia.org/wiki/Airy_function
+    """
+
+    #We catch the case with no alpha
+    if x==None:
+        x=alpha
+        return airy_ai_simple(x, **kwds)
+    #We raise an error if there are too many arguments
+    if len(args) > 0:
+        raise TypeError("Symbolic function airy_ai takes at most 3 arguments (%s given)" %
+                        (len(args)+3))
+
+    #We take care of all other cases.
+    if hold_derivative:
+        return airy_ai_general(alpha,x,**kwds)
+    elif alpha==0:
+        return airy_ai_simple(x, **kwds)
+    elif alpha==1:
+        return airy_ai_prime(x, **kwds)
+    elif alpha>1:
+        #we use a different variable here because if x is a
+        #particular value, we would be differentiating a constant
+        #which would return 0. What we want is the value of
+        #the derivative at the value and not the derivative of
+        #a particular value of the function.
+        v=var('v')
+        return derivative(airy_ai_simple(v,**kwds),v,alpha).subs(v=x)
+    else:
+        return airy_ai_general(alpha,x,**kwds)
+
+########################################################################
+########################################################################
+
+class FunctionAiryBiGeneral(BuiltinFunction):
+    def __init__(self):
+        r"""
+        The generalized derivative of the Airy Bi function
+
+        INPUT:
+        
+        - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to 
+          `z`. For `\alpha = n = 1,2,3,\ldots` this gives the derivative 
+          `\operatorname{Bi}^{(n)}(z)`, and for `\alpha = -n = -1,-2,-3,\ldots` this gives 
+          the `n`-fold iterated integral
+         
+        .. math ::
+         
+            f_0(z) = \operatorname{Bi}(z)
+         
+            f_n(z) = \int_0^z f_{n-1}(t) dt.
+
+        - ``x``.- The argument of the function.
+
+        - ``hold``.- Whether or not to stop from returning higher derivatives in terms of 
+          `\operatorname{Bi}(x)` and `\operatorname{Bi}'(x)`.
+        
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: airy_bi(-2,x)
+            airy_bi(-2, x)
+            sage: derivative(airy_bi(-2,x),x)
+            airy_bi(-1, x)
+            sage: airy_bi(n,x)
+            airy_bi(n, x)
+            sage: derivative(airy_bi(n,x),x)
+            airy_bi(n + 1, x)
+            sage: airy_bi(2,x,True)
+            airy_bi(2, x)
+            sage: derivative(airy_bi(2,x,hold_derivative=True),x)
+            airy_bi(3, x)
+        """
+
+        BuiltinFunction.__init__(self, "airy_bi", nargs=2,
+            latex_name=r"\operatorname{Bi}")
+
+    def _derivative_(self, alpha, *args, **kwds):
+        """
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: derivative(airy_bi(n,x),x) # indirect doctest
+            airy_bi(n + 1, x)
+        """
+
+        x=args[0]
+        return airy_bi_general(alpha+1,x)
+
+    def _eval_(self, alpha, *args):
+        """
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: airy_bi(-2,1.0) # indirect doctest
+            0.388621540699059
+            sage: airy_bi(n,1.0) # indirect doctest
+            airy_bi(n, 1.00000000000000)
+        """
+
+        
+        x=args[0]
+        if not isinstance(x,Expression) and not isinstance(alpha,Expression):
+            if is_inexact(x):
+                return self._evalf_(alpha,x)
+        else:
+            return None
+
+    def _evalf_(self, alpha, x, **kwargs):
+        """
+        EXAMPLES::
+
+            sage: airy_bi(-2,1.0) # indirect doctest
+            0.388621540699059
+
+        """
+        algorithm = kwargs.get('algorithm', None) or 'mpmath'
+        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
+        prec = parent.prec() if hasattr(parent, 'prec') else 53
+        if algorithm == 'mpmath':
+            import mpmath
+            from sage.libs.mpmath import utils as mpmath_utils
+            return mpmath_utils.call(mpmath.airybi, x, derivative=alpha, parent=parent)
+        elif algorithm == 'maxima':
+            raise NotImplementedError("general case not available in maxima")
+        else:
+            raise ValueError("unknown algorithm")
+
+class FunctionAiryBiSimple(BuiltinFunction):
+    def __init__(self):
+        """
+        The class for the Airy Bi function.
+
+        EXAMPLES::
+
+            sage: from sage.functions.airy import FunctionAiryBiSimple
+            sage: airy_bi_simple= FunctionAiryBiSimple()
+            sage: f=airy_bi_simple(x); f
+            airy_bi(x)
+        """
+
+        BuiltinFunction.__init__(self, "airy_bi", latex_name=r'\operatorname{Bi}',
+                conversions=dict(mathematica='AiryBi', maxima='airy_bi'))
+
+    def _derivative_(self, x, diff_param=None):
+        """
+        EXAMPLES::
+        
+            sage: derivative(airy_bi(x),x) # indirect doctest
+            airy_bi_prime(x)
+        """
+        return airy_bi_prime(x)
+
+    def _eval_(self, x):
+        """
+        EXAMPLES::
+        
+            sage: airy_bi(0) # indirect doctest
+            1/3*3^(5/6)/gamma(1/3)
+            sage: airy_bi(0.0) # indirect doctest
+            0.614926627446001
+            sage: airy_bi(I) # indirect doctest
+            airy_bi(I)
+            sage: airy_bi(1.0*I) # indirect doctest
+            0.648858208330395 + 0.344958634768048*I
+        """
+
+        if not isinstance(x,Expression):
+            if is_inexact(x):
+                return self._evalf_(x)
+            elif x==0:
+                one_sixth = ZZ(1)/6
+                return 1/(3**(one_sixth)*gamma(2*one_sixth))
+        else:
+            return None
+
+    def _evalf_(self, x, **kwargs):
+        """
+        EXAMPLES::
+        
+            sage: airy_bi(0.0) # indirect doctest
+            0.614926627446001
+            sage: airy_bi(1.0*I) # indirect doctest
+            0.648858208330395 + 0.344958634768048*I
+
+        We can use several methods for numerical evaluation::
+
+            sage: airy_bi(3).n(algorithm='maxima')
+            14.0373289637
+            sage: airy_bi(3).n(algorithm='mpmath')
+            14.0373289637302
+            sage: airy_bi(3).n(algorithm='mpmath', prec=100)
+            14.037328963730232031740267314
+
+        TESTS::
+
+            sage: airy_bi(3).n(algorithm='maxima', prec=100)
+            Traceback (most recent call last):
+            ...
+            ValueError: for the maxima algorithm the precision must be 53
+
+        """
+        algorithm = kwargs.get('algorithm', None) or 'mpmath'
+        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
+        prec = parent.prec() if hasattr(parent, 'prec') else 53
+        if algorithm == 'mpmath':
+            import mpmath
+            from sage.libs.mpmath import utils as mpmath_utils
+            return mpmath_utils.call(mpmath.airybi, x, parent=parent)
+        elif algorithm == 'maxima':
+            if prec != 53:
+                raise ValueError, "for the maxima algorithm the precision must be 53"
+            return RDF(meval("airy_bi(%s)" % RDF(x)))
+        else:
+            raise ValueError("unknown algorithm")
+
+class FunctionAiryBiPrime(BuiltinFunction):
+    def __init__(self):
+        """
+        The derivative of the Airy Bi function; see :func:`airy_bi`
+        for the full documentation.
+
+        EXAMPLES::
+
+            sage: x,n=var('x n')
+            sage: airy_bi_prime(x) # indirect doctest
+            airy_bi_prime(x)
+            sage: airy_bi_prime(0) # indirect doctest
+            3^(1/6)/gamma(1/3)
+        """
+
+        BuiltinFunction.__init__(self, "airy_bi_prime",
+                latex_name=r"\operatorname{Bi}'",
+                conversions=dict(mathematica='AiryBiPrime', maxima='airy_dbi'))
+
+    def _derivative_(self, x, diff_param=None):
+        """
+        EXAMPLES::
+
+            sage: derivative(airy_bi_prime(x),x) # indirect doctest
+            x*airy_bi(x)
+        """
+        return x*airy_bi_simple(x)
+
+    def _eval_(self, x):
+        """
+        EXAMPLES::
+
+            sage: airy_bi(1,0) # indirect doctest
+            3^(1/6)/gamma(1/3)
+            sage: airy_bi(1,0.0) # indirect doctest
+            0.448288357353826
+        """
+        if not isinstance(x,Expression):
+            if is_inexact(x):
+                return self._evalf_(x)
+            elif x==0:
+                one_sixth = ZZ(1)/6
+                return 3**(one_sixth)/gamma(2*one_sixth)
+        else:
+            return None
+
+    def _evalf_(self, x, **kwargs):
+        """
+        EXAMPLES::
+        
+            sage: airy_bi(1,0.0) # indirect doctest
+            0.448288357353826
+
+        We can use several methods for numerical evaluation::
+
+            sage: airy_bi(1,4).n(algorithm='maxima')
+            161.926683505
+            sage: airy_bi(1,4).n(algorithm='mpmath')
+            161.926683504613
+            sage: airy_bi(1,4).n(algorithm='mpmath', prec=100)
+            161.92668350461340184309492429
+
+        TESTS::
+
+            sage: airy_bi(1,4).n(algorithm='maxima', prec=70)
+            Traceback (most recent call last):
+            ...
+            ValueError: for the maxima algorithm the precision must be 53
+
+        """
+        algorithm = kwargs.get('algorithm', None) or 'mpmath'
+        parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x)
+        prec = parent.prec() if hasattr(parent, 'prec') else 53
+        if algorithm == 'mpmath':
+            import mpmath
+            from sage.libs.mpmath import utils as mpmath_utils
+            return mpmath_utils.call(mpmath.airybi, x, derivative=1, parent=parent)
+        elif algorithm == 'maxima':
+            if prec != 53:
+                raise ValueError, "for the maxima algorithm the precision must be 53"
+            return RDF(meval("airy_dbi(%s)" % RDF(x)))
+        else:
+            raise ValueError("unknown algorithm")
+
+airy_bi_general=FunctionAiryBiGeneral()
+airy_bi_simple= FunctionAiryBiSimple()
+airy_bi_prime=  FunctionAiryBiPrime()
+
+def airy_bi(alpha,x=None, hold_derivative=False, *args, **kwds):
+    r"""
+    The Airy Bi function `\operatorname{Bi}(x)` is one of the two
+    linearly independent solutions to the Airy differental equation
+    `f''(z) +f(z)x=0`, defined by the initial conditions:
+
+    .. math ::
+        \operatorname{Bi}(0)=\frac{1}{3^{1/6} \Gamma(\frac{2}{3})},
+
+        \operatorname{Bi}'(0)=\frac{3^{1/6}}{ \Gamma(\frac{1}{3})}.
+
+    Another way to define the Airy Bi function is:
+
+    .. math::
+        \operatorname{Bi}(x)=\frac{1}{\pi}\int_0^\infty
+        \left[ \exp\left( xt -\frac{t^3}{3} \right)
+        +\sin\left(xt + \frac{1}{3}t^3\right) \right ] dt.
+
+    INPUT:
+    
+    - ``alpha``.- Return the `\alpha`-th order fractional derivative with respect to `z`. 
+      For `\alpha = n = 1,2,3,\ldots` this gives the derivative `\operatorname{Bi}^{(n)}(z)`, 
+      and for `\alpha = -n = -1,-2,-3,\ldots` this gives the `n`-fold iterated integral
+     
+    .. math ::
+     
+        f_0(z) = \operatorname{Bi}(z)
+     
+        f_n(z) = \int_0^z f_{n-1}(t) dt.
+
+    - ``x``.- The argument of the function.
+
+    - ``hold_derivative``.- Whether or not to stop from returning higher derivatives in 
+      terms of `\operatorname{Bi}(x)` and `\operatorname{Bi}'(x)`.
+    
+
+    EXAMPLES::
+
+        sage: n,x=var('n x')
+        sage: airy_bi(x)
+        airy_bi(x)
+
+    It can return derivatives or integrals::
+    
+        sage: airy_bi(1,x)
+        airy_bi_prime(x)
+        sage: airy_bi(2,x)
+        x*airy_bi(x)
+        sage: airy_bi(2,x,True)
+        airy_bi(2, x)
+        sage: airy_bi(-2,x)
+        airy_bi(-2, x)
+        sage: airy_bi(n, x)
+        airy_bi(n, x)
+
+    It can be evaluated symbolically or numerically for real or complex values::
+
+        sage: airy_bi(0)
+        1/3*3^(5/6)/gamma(1/3)
+        sage: airy_bi(0.0)
+        0.614926627446001
+        sage: airy_bi(I)
+        airy_bi(I)
+        sage: airy_bi(1.0*I)
+        0.648858208330395 + 0.344958634768048*I
+
+    The functions can be evaluated numerically using either mpmath (the default)
+    or maxima, where mpmath can compute the values to arbitrary precision::
+
+        sage: airy_bi(2).n(prec=100)
+        3.2980949999782147102806044252
+        sage: airy_bi(2).n(algorithm='mpmath',prec=100)
+        3.2980949999782147102806044252
+        sage: airy_bi(2).n(algorithm='maxima')
+        3.29809499998
+
+    And the derivatives can be evaluated::
+
+        sage: airy_bi(1,0)
+        3^(1/6)/gamma(1/3)
+        sage: airy_bi(1,0.0)
+        0.448288357353826
+
+    Plots::
+
+        sage: plot(airy_bi(x),(x,-10,5))+plot(airy_bi_prime(x),(x,-10,5),color='red')
+        
+    **References**
+     
+    - Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 10"
+
+    - http://en.wikipedia.org/wiki/Airy_function
+    """
+
+    #We catch the case with no alpha
+    if x==None:
+        x=alpha
+        return airy_bi_simple(x, **kwds)
+    #We raise an error if there are too many arguments
+    if len(args) > 0:
+        raise TypeError("Symbolic function airy_ai takes at most 3 arguments (%s given)" %
+                        (len(args)+3))
+
+    #We take care of all other cases.
+    if hold_derivative:
+        return airy_bi_general(alpha,x,**kwds)
+    elif alpha==0:
+        return airy_bi_simple(x, **kwds)
+    elif alpha==1:
+        return airy_bi_prime(x, **kwds)
+    elif alpha>1:
+        #we use a different variable here because if x is a
+        #particular value, we would be differentiating a constant
+        #which would return 0. What we want is the value of
+        #the derivative at the value and not the derivative of
+        #a particular value of the function.
+        v=var('v')
+        return derivative(airy_bi_simple(v,**kwds),v,alpha).subs(v=x)
+    else:
+        return airy_bi_general(alpha,x,**kwds)
+

sage/functions/all.py

@@ -37 +37 @@
                      inverse_jacobi,
                      lngamma, exp_int, error_fcn, elliptic_e,
                      elliptic_f, elliptic_ec, elliptic_eu,
-                     elliptic_kc, elliptic_pi, elliptic_j,
-                     airy_ai, airy_bi)
+                     elliptic_kc, elliptic_pi, elliptic_j)
                         
 from orthogonal_polys import (chebyshev_T,
                               chebyshev_U,
@@ -64 +63 @@
                          kronecker_delta)
 
 from min_max import max_symbolic, min_symbolic
+
+from airy import airy_ai, airy_ai_prime, airy_bi, airy_bi_prime

sage/functions/special.py

@@ -72 +72 @@
    kind; they are two linearly independent solutions of Bessel's
    differential equation. They are named for Hermann Hankel.
 
--  Airy function The function `Ai(x)` and the related
-   function `Bi(x)`, which is also called an Airy function,
-   are solutions to the differential equation
-
-   
-   .. math::
-
-         y'' - xy = 0, 
-
-   known as the Airy equation. They belong to the class of 'Bessel functions of
-   fractional order'. The initial conditions
-   `Ai(0) = (\Gamma(2/3)3^{2/3})^{-1}`,
-   `Ai'(0) = -(\Gamma(1/3)3^{1/3})^{-1}` define
-   `Ai(x)`. The initial conditions
-   `Bi(0) = 3^{1/2}Ai(0)`, `Bi'(0) = -3^{1/2}Ai'(0)`
-   define `Bi(x)`.
-
-   They are named after the British astronomer George Biddell Airy.
-
 -  Spherical harmonics: Laplace's equation in spherical coordinates
    is:
    
@@ -335 +316 @@
 
 - http://en.wikipedia.org/wiki/Bessel_function
 
-- http://en.wikipedia.org/wiki/Airy_function
-
 - http://en.wikipedia.org/wiki/Spherical_harmonics
 
 - http://en.wikipedia.org/wiki/Helmholtz_equation
@@ -406 +385 @@
 
     Then after using one of these functions, it changes::
 
-        sage: from sage.functions.special import airy_ai
-        sage: airy_ai(1.0)
-        0.135292416313
+        sage: from sage.functions.special import spherical_harmonic
+        sage: spherical_harmonic(3,2,1,2)
+        15/4*sqrt(7/30)*e^(4*I)*sin(1)^2*cos(1)/sqrt(pi)
         sage: sage.functions.special._done
         True
     """
@@ -426 +405 @@
     
     TEST::
 
-        sage: from sage.functions.special import airy_ai
-        sage: airy_bi(1.0)
-        1.20742359495
+        sage: from sage.functions.special import spherical_bessel_J
+        sage: spherical_bessel_J(2.,3.)
+        0.298637497076
+
     """
     return maxima(x).sage()
 
@@ -579 +559 @@
 
     return NewMaximaFunction()
 
-
-def airy_ai(x):
-   r"""
-   The function `Ai(x)` and the related function `Bi(x)`,
-   which is also called an *Airy function*, are
-   solutions to the differential equation
-
-   .. math::
-
-      y'' - xy = 0,
-
-   known as the *Airy equation*. The initial conditions
-   `Ai(0) = (\Gamma(2/3)3^{2/3})^{-1}`,
-   `Ai'(0) = -(\Gamma(1/3)3^{1/3})^{-1}` define `Ai(x)`.
-   The initial conditions `Bi(0) = 3^{1/2}Ai(0)`,
-   `Bi'(0) = -3^{1/2}Ai'(0)` define `Bi(x)`.
-
-   They are named after the British astronomer George Biddell Airy.
-   They belong to the class of "Bessel functions of fractional order".
-
-   EXAMPLES::
-   
-       sage: airy_ai(1.0)        # last few digits are random
-       0.135292416312881400
-       sage: airy_bi(1.0)        # last few digits are random
-       1.20742359495287099
-
-   REFERENCE:
-
-   - Abramowitz and Stegun: Handbook of Mathematical Functions,
-     http://www.math.sfu.ca/~cbm/aands/
-
-   - http://en.wikipedia.org/wiki/Airy_function
-   """
-   _init()
-   return RDF(meval("airy_ai(%s)"%RDF(x)))
-
-def airy_bi(x):
-   r"""
-   The function `Ai(x)` and the related function `Bi(x)`,
-   which is also called an *Airy function*, are
-   solutions to the differential equation
-
-   .. math::
-
-      y'' - xy = 0,
-
-   known as the *Airy equation*. The initial conditions
-   `Ai(0) = (\Gamma(2/3)3^{2/3})^{-1}`,
-   `Ai'(0) = -(\Gamma(1/3)3^{1/3})^{-1}` define `Ai(x)`.
-   The initial conditions `Bi(0) = 3^{1/2}Ai(0)`,
-   `Bi'(0) = -3^{1/2}Ai'(0)` define `Bi(x)`.
-
-   They are named after the British astronomer George Biddell Airy.
-   They belong to the class of "Bessel functions of fractional order".
-
-   EXAMPLES::
-   
-       sage: airy_ai(1)        # last few digits are random
-       0.135292416312881400
-       sage: airy_bi(1)        # last few digits are random
-       1.20742359495287099
-
-   REFERENCE:
-
-   - Abramowitz and Stegun: Handbook of Mathematical Functions,
-     http://www.math.sfu.ca/~cbm/aands/
-
-   - http://en.wikipedia.org/wiki/Airy_function
-   """
-   _init()
-   return RDF(meval("airy_bi(%s)"%RDF(x)))
-
-
 def bessel_I(nu,z,algorithm = "pari",prec=53):
     r"""
     Implements the "I-Bessel function", or "modified Bessel function,