bessel_K function is broken

[sagetrac-bober]

Currently we have

sage: bessel_K(10 * I, 10)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/bober/sage-3.0.2/devel/sage-bober/sage/functions/<ipython console> in <module>()

/home/bober/sage/local/lib/python2.5/site-packages/sage/functions/special.py in bessel_K(nu, z, algorithm, prec)
    586         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)
    590         return b

/home/bober/sage-3.0.2/devel/sage-bober/sage/functions/real_mpfr.pyx in sage.rings.real_mpfr.RealField.__call__ (sage/rings/real_mpfr.c:3138)()

/home/bober/sage-3.0.2/devel/sage-bober/sage/functions/real_mpfr.pyx in sage.rings.real_mpfr.RealNumber._set (sage/rings/real_mpfr.c:5905)()

TypeError: Unable to convert x (='0.000000098241574381992468+0.E-161*I') to real number.
sage: bessel_K(10 * I, 10)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/bober/sage-3.0.2/devel/sage-bober/sage/functions/<ipython console> in <module>()

/home/bober/sage/local/lib/python2.5/site-packages/sage/functions/special.py in bessel_K(nu, z, algorithm, prec)
    586         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)
    590         return b

/home/bober/sage-3.0.2/devel/sage-bober/sage/functions/real_mpfr.pyx in sage.rings.real_mpfr.RealField.__call__ (sage/rings/real_mpfr.c:3138)()

/home/bober/sage-3.0.2/devel/sage-bober/sage/functions/real_mpfr.pyx in sage.rings.real_mpfr.RealNumber._set (sage/rings/real_mpfr.c:5905)()

TypeError: Unable to convert x (='0.000000098241574381992468+0.E-161*I') to real number.

In this case the result actually should be a real number, so we fix this by discarding the imaginary part of the result from pari. In other cases, however, the result is actually a complex number, and we shouldn't always be attempting to cast it to a real number (which the attached patch also fixes).

CC: @burcin @kcrisman @benjaminfjones

Component: calculus

Keywords: bessel, bessel_K

Reviewer: Karl-Dieter Crisman, Benjamin Jones

Issue created by migration from https://trac.sagemath.org/ticket/3426

[sagetrac-bober]

Attachment: kbessel_fixes.patch.gz

[craigcitro]

Changed keywords from bessel, bessel_K to bessel, bessel_K, editor_gfurnish

[sagetrac-bober]

comment:4

Regarding bessel_K being real for real argument and real or imaginary order, see, e.g. the appendix to H. Then, Maass cusp forms for large eigenvalues, Math. Comp. Volume 74, Number 249, pp. 363 - 381: "The K-Bessel function K_ir(x) is ... real for real arguments x and real or imaginary order ir."

[sagetrac-gregorybard]

comment:5

I can say that I agree that this is 100% correct.

[rishikesha]

comment:9

I think a solution of the following type would be better.

try:
        from sage.libs.pari.all import pari
        RR,a = _setup(prec)
        b = RR(pari(nu).besselk(z))
        pari.set_real_precision

except TypeError:
        CC,a = _setup(prec)
        b = CC(pari(nu).besselk(z))
        pari.set_real_precision(a)

[rishikesha]

comment:10

Probably the correct code would be

try:
        from sage.libs.pari.all import pari
        RR,a = _setup(prec)
        b = RR(pari(nu).besselk(z))
        pari.set_real_precision

except TypeError:
        CC,a = _setup_CC(prec)
        b = CC(pari(nu).besselk(z))
        pari.set_real_precision(a)

[sagetrac-mabshoff]

comment:11

Since Rishi commented on this it might be a good idea to discuss his comments.

Cheers,

Michael

[JohnCremona]

comment:26

I note that Alex added himself to the CC for this. Whatever is done for this issue absolutely must take into account the work done for #4096, so at the least I suggest that the author of this patch looks at that one and reworks this. Anything relying on Sage/pari precision questions is likely to be useless otherwise.

[aghitza]

comment:27

I looked up the definition and properties of the Bessel functions in
several references (Section 7.2 of the "Bateman Manuscript Project",
for instance).

I uploaded a brand new patch that implements the behavior described
there, namely returning a real number if the result is theoretically
known to be real, and a complex number otherwise. I added doctests
that document this behavior, and checked all of them against
Mathematica. I did this for all three Bessel functions that are
implemented in special.py using Pari, namely J, K, and I. I also put in a workaround for a silly Pari buglet that
complains about negative integer values of nu.

In the process I uncovered a couple of unrelated issues with
special.py and Bessel functions, for which I'll open separate tickets.

The patch is made against 3.1.3.alpha2.

[dandrake]

comment:28

Uh oh:

sage: bessel_J(0,0)
---------------------------------------------------------------------------
PariError                                 Traceback (most recent call last)

/home/drake/.sage/temp/klee/32521/_home_drake__sage_init_sage_0.py in <module>()
----> 1 
      2 
      3 
      4 
      5 

/opt/sage-3.1.3.alpha2/local/lib/python2.5/site-packages/sage/functions/special.pyc in bessel_J(nu, z, algorithm, prec)
    570             z = C(z)
    571             K = C
--> 572         pari_bes = pari(nu).besselj(z, precision=prec)
    573         if K is R:
    574             return fudge*K(pari_bes.real())

/opt/sage-3.1.3.alpha2/local/lib/python2.5/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:34414)()
   7865 
   7866 
-> 7867 
   7868 
   7869 

PariError:  (8)

Doing bessel_J(0, 0) in 3.1.2 works fine. I get similar errors with this patch for other Bessel functions, too.

[aghitza]

comment:29

Thanks for catching this. It actually comes from a bug in Pari:

? besselj(0, 0)
%1 = 1.000000000000000000000000000
? besselj(0.E-19, 0)
  *** besselj: gpow: 0 to a non positive exponent.

I've reported it upstream, but I will post a patch with a workaround while we wait.

[aghitza]

apply instead of the previous patch

[aghitza]

comment:30

Attachment: trac3426-fix-bessel-fns.patch.gz

OK, I've replaced my patch with one that fixes the issue reported by Dan.

[dandrake]

comment:31

Now bessel_J(0,0) works but I'm seeing other problems. I'm concentrating here on the "K" functions since that's what this ticket is about. This is all with the current patch applied to 3.1.4.

• bessel_K(0, -1) doesn't work at all; it just soaks up all the CPU and doesn't return. The correct answer is about 0.421024438240708 - 3.97746326050642*I.

• bessel_K(-1*I - 1, 0) gives an uninformative Pari error; the function isn't defined there. Mathematica says "ComplexInfinity"; can we give a better error message? Even just "function not defined at 0" would be fine.

• bessel_K(-1, -1) says it can't convert -0.601907230197235-1.77549968921218*I to a real number, which I agree with. :) It looks like Pari has the right answer but we're trying to convert it to a real when we shouldn't be.

• bessel_K(0, -1 - I) says "incorrect type", but Mathematica and Maple evaluate it just fine.

[rlmill]

comment:33

See #4626, which at least fixes the bessel_J problem.

[aghitza]

comment:34

This ticket is a huge mess :)

I now think that we should just use mpmath to evaluate Bessel functions, see

http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html

For the examples that Dan gave:

sage: from mpmath import *
sage: mp.dps = 25; mp.pretty = True
sage: besselk(0, -1)
(0.4210244382407083333356274 - 3.97746326050642263725661j)
sage: besselk(-1*I - 1, 0)
+inf
sage: besselk(-1, -1)
(-0.60190723019723457473754 - 1.775499689212180946878577j)
sage: besselk(0, -1-I)
(-1.479697108749625193260947 + 2.588306443392007370808151j)

[aghitza]

Changed keywords from bessel, bessel_K, editor_gfurnish to bessel, bessel_K

[rishikesha]

comment:35

From my quick experiments with the issues Bober was dealing with a year ago, I see that do not arise if we use mpmath, even when I set the precision to 5000.

I agree with Alex Ghitza.

I should say that the bessels functions for non integer indices have always bothered me. I believe computing will involve a log, and how do you consistently choose a branch.

[dandrake]

comment:36

I also agree with Alex -- both in that this ticket is a mess, and that we should just use mpmath. I'll see if I can work up a patch (although anyone else who wants to should feel free to do it).

[rishikesha]

comment:37

With some experiments, I saw that the branch of the log taken is the negative real axis. We should mention this in the documentation when it is implemented. I believe ddrake is working up a patch.

[dandrake]

comment:39

Here are some more problems with bessel_K: http://groups.google.com/group/sage-support/browse_thread/thread/84e5cd4fb1381ad5

In resolving this ticket, we should make sure to have doctests related to those problems!

[kcrisman]

comment:41

This would most likely be fixed by #4102.

[benjaminfjones]

comment:42

Yes, it will be fixed by #4102. I'll make a note to add a related doctest to that effect.

[kcrisman]

Reviewer: Karl-Dieter Crisman, Benjamin Jones

[kcrisman]

comment:43

Everything here is now in a doctest in #4102, including the stuff in the thread from three (!) years ago.