Ticket #3426: kbessel_fixes.patch

File kbessel_fixes.patch, 2.1 KB (added by bober, 14 years ago)

  • sage/functions/special.py 

    # HG changeset patch
# User Jonathan Bober <bober@umich.edu>
# Date 1213480819 25200
# Node ID 15b19e6cb696e718b8311aca72e519eab9ca7b1e
# Parent  6a6766d05f3bc23e52c4ce7477132b9a7ce1607d
Some fixes to the bessel_K function. Previously this didn't work when
the result was a complex number. Now we ignore the imaginary part of
the result when we know that the answer should be real and we cast the
result to a complex number instead of a real number when it genuinely
is a complex number.

diff -r 6a6766d05f3b -r 15b19e6cb696 sage/functions/special.py
    a b def bessel_K(nu,z,algorithm="pari",prec= 
    572572        0.60
    573573        sage: bessel_K(1,1,"pari",100)
    574574        0.60190723019723457473754000154
     575        sage: bessel_K(10 * I, 10, prec = 100)
     576        9.8241574381992468026984622571e-8
     577        sage: bessel_K(10*I + 19, 10, prec=100)
     578        0.21750889400945826164540481362 + 4.3923792070483775870514070521*I
    575579    """
    576580    if algorithm=="scipy":
    577581        if prec != 53:
    … … def bessel_K(nu,z,algorithm="pari",prec= 
    584588        return sage_eval(ans)
    585589    elif algorithm == 'pari':
    586590        from sage.libs.pari.all import pari
    587         RR,a = _setup(prec)
    588         b = RR(pari(nu).besselk(z))
    589         pari.set_real_precision(a)
     591        # If vu is either real or purely imaginary, and z is real, then
     592        # the answer is known to be real, and the imaginary part of the
     593        # answer returned from pari is just numerical noise, so
     594        # we just discard it and return a real answer.
     595        #
     596        # Otherwise, we need to retain both the real and imaginary parts of the
     597        # answer and return a complex number.
     598        if (real(nu) == 0 or imag(nu) == 0) and imag(z) == 0:
     599            RR, a = _setup(prec)
     600            b = RR(real(pari(nu).besselk(z)))
     601            pari.set_real_precision(a)
     602        else:
     603            CC,a = _setup_CC(prec)
     604            b = CC(pari(nu).besselk(z))
     605            pari.set_real_precision(a)
    590606        return b
    591607    else:
    592608        raise ValueError, "unknown algorithm '%s'"%algorithm