Ticket #3426: trac3426-fix-bessel-fns.patch

File trac3426-fix-bessel-fns.patch, 6.6 KB (added by AlexGhitza, 11 years ago)

apply instead of the previous patch

  • sage/functions/special.py

    # HG changeset patch
    # User Alexandru Ghitza <aghitza@alum.mit.edu>
    # Date 1222843216 -36000
    # Node ID ed46740a68d3602b822352abf6378c70281b2b50
    # Parent  4805d6f2365704ea92de52d00a97375299950add
    trac 3426: fix Bessel functions in special.py
    
    diff -r 4805d6f23657 -r ed46740a68d3 sage/functions/special.py
    a b  
    380380    index (or "order") nu and argument z.
    381381
    382382    INPUT:
    383         nu -- a real (or complex, for pari) number
    384         z  -- a real (positive)
    385         algorithm - "pari" or "maxima" or "scipy"
    386         prec - real precision (for Pari only)
     383        nu -- a real or complex number
     384        z  -- a real or complex number
     385        algorithm -- "pari" (default) or "maxima" or "scipy"; note
     386            that "maxima" and "scipy" can only deal with real nu and
     387            positive real z
     388        prec -- real precision (for Pari only)
    387389   
    388390    DEFINITION:
    389391    \begin{verbatim}
     
    428430        sage: bessel_I(1,1,"scipy")
    429431        0.565159103992...
    430432
     433    A corner case:
     434        sage: bessel_I(0,0)
     435        1.00000000000000
     436
     437    If both nu and z are real, then bessel_I is known to be real if z
     438    is positive or if nu is an integer:
     439        sage: bessel_I(-2.1, 5.3)
     440        22.5840990831073
     441        sage: bessel_I(-2.1, -5.3)
     442        21.4787545976446 - 6.97887041932783*I
     443        sage: bessel_I(-2, -5.3)
     444        23.5424857046048
     445        sage: bessel_I(10*I, 10)
     446        1.19390297337705e6 - 688460.429969423*I
    431447    """
    432448    if algorithm=="pari":
    433449        from sage.libs.pari.all import pari
     450        R = RealField(prec)
     451        # this is a workaround for a bug in Pari 2.3.3
     452        if nu.is_zero() and z.is_zero():
     453            return R(1)
     454        C = ComplexField(prec)
    434455        try:
    435             R = RealField(prec)
     456            nu = ZZ(nu)
     457            # pari doesn't like negative integers here,
     458            # but luckily I(-nu, z) = I(nu, z) if nu is integral
     459            nu = nu.abs()
     460            K = R
     461        except TypeError:
     462            K = C
     463        try:
    436464            nu = R(nu)
    437465            z = R(z)
     466            if (z > 0):
     467                K = R
    438468        except TypeError:
    439             C = ComplexField(prec)
    440469            nu = C(nu)
    441470            z = C(z)
    442471            K = C
    443         K = z.parent()
    444         return K(pari(nu).besseli(z, precision=prec))
     472        pari_bes = pari(nu).besseli(z, precision=prec)
     473        if K is R:
     474            return K(pari_bes.real())
     475        else:
     476            return K(pari_bes)
    445477    elif algorithm=="scipy":
    446478        if prec != 53:
    447479            raise ValueError, "for the scipy algorithm the precision must be 53"
     
    508540        0.0583793793051868
    509541        sage: bessel_J(3,10,"scipy")
    510542        0.0583793793052...
     543
     544    A corner case:
     545        sage: bessel_I(0,0)
     546        1.00000000000000
     547
     548    If both nu and z are real, then bessel_J is known to be real if z
     549    is positive or if nu is an integer:
     550        sage: bessel_J(-2.1, 5.3)
     551        -0.123120872446240
     552        sage: bessel_J(-2.1, -5.3)
     553        -0.117094908031941 + 0.0380464419481583*I
     554        sage: bessel_J(-2, -5.3)
     555        -0.0547481464525993
     556        sage: bessel_J(10*I, 10)
     557        -119728.409042333 - 693848.987632644*I
    511558    """
    512559    if algorithm=="pari":
    513560        from sage.libs.pari.all import pari
     561        R = RealField(prec)
     562        # this is a workaround for a bug in Pari 2.3.3
     563        if nu.is_zero() and z.is_zero():
     564            return R(1)
     565        C = ComplexField(prec)
     566        fudge = 1
    514567        try:
    515             R = RealField(prec)
     568            nu = ZZ(nu)
     569            # pari doesn't like negative integers here,
     570            # but luckily J(-nu, z) = (-1)**nu J(nu, z) if nu is integral
     571            if nu < 0:
     572                nu = -nu
     573                fudge = (-1)**nu
     574            K = R
     575        except TypeError:
     576            K = C
     577        try:
    516578            nu = R(nu)
    517579            z = R(z)
     580            if (z > 0):
     581                K = R
    518582        except TypeError:
    519             C = ComplexField(prec)
    520583            nu = C(nu)
    521584            z = C(z)
    522585            K = C
    523         K = z.parent()
    524         return K(pari(nu).besselj(z, precision=prec))
     586        pari_bes = pari(nu).besselj(z, precision=prec)
     587        if K is R:
     588            return fudge*K(pari_bes.real())
     589        else:
     590            return fudge*K(pari_bes)
    525591    elif algorithm=="scipy":
    526592        if prec != 53:
    527593            raise ValueError, "for the scipy algorithm the precision must be 53"
     
    553619    otherwise.
    554620   
    555621    Sometimes bessel_K(nu,z) is denoted K_nu(z) in the
    556     literature. In Pari, nu can be complex and
    557     z must be real and positive.
     622    literature. In Pari, both nu and z can be complex; in SciPy they
     623    must be real and z mush be positive.
    558624
    559625    EXAMPLES:
    560626        sage: bessel_K(3,2,"scipy")
     
    567633        0.60
    568634        sage: bessel_K(1,1,"pari",100)
    569635        0.60190723019723457473754000154
     636
     637        sage: bessel_K(0,0)
     638        +Infinity
     639
     640    If both nu and z are real, then bessel_K is known to be real if
     641    z is positive:
     642        sage: bessel_K(-2.1, 5.3)
     643        0.00388720082599840
     644        sage: bessel_K(-2.1, -5.3)
     645        0.00369694767571369 - 70.9488385563182*I
     646        sage: bessel_K(-2, -5.3)
     647        0.00375340283798982 - 73.9609001368291*I
     648        sage: bessel_K(10*I, 10)
     649        9.82415743819925e-8
     650        sage: bessel_K(2, -1.121)
     651        1.23414145962983 - 0.547231658266306*I
    570652    """
    571653    if algorithm=="scipy":
    572654        if prec != 53:
     
    579661        return sage_eval(ans)
    580662    elif algorithm == 'pari':
    581663        from sage.libs.pari.all import pari
     664        R = RealField(prec)
     665        C = ComplexField(prec)
     666        if z.is_zero() and nu.is_zero():
     667            from sage.rings.infinity import Infinity
     668            return Infinity
     669        #TODO: sort out cases z = 0, nu != 0
    582670        try:
    583             R = RealField(prec)
    584671            nu = R(nu)
     672            # K(-nu, z) = K(nu, z) for all nu
     673            nu = nu.abs()
    585674            z = R(z)
     675            if (z > 0):
     676                K = R
     677            else:
     678                K = C
    586679        except TypeError:
    587             C = ComplexField(prec)
    588680            nu = C(nu)
    589681            z = C(z)
    590682            K = C
    591         K = z.parent()
    592         return K(pari(nu).besselk(z, precision=prec))
     683        pari_bes = pari(nu).besselk(z, precision=prec)
     684        if K is R:
     685            return K(pari_bes.real())
     686        else:
     687            return K(pari_bes)
    593688    else:
    594689        raise ValueError, "unknown algorithm '%s'"%algorithm
    595690