[with patch, needs work] bessel_K function is broken

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."

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

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)

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)

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

Cheers,

Michael

The try code does not in my opinion work. The issue here is correcting numerical instabilities, which attempting to coerce into RR will not do. Positive review still.

I am pretty confident that the try code works. Consider the following gp output. The current patch will try to coerce this to RR. If pari is right then this input after applying the patch will give a TypeError?.

? besselk(2,-1.121)
%1 = 1.234141459629829380224386595 - 0.5472316582663064541169798027*I

We probably eliminate the double evaluation of bessel function using pari which my current suggestion does.

I don't understand what double evaluation of the bessel function you are talking about. Furthermore the entire point is that Pari can give wrong answers in some cases because of numerical cases. Therefore we want to manually correct the numerical noise. Trying to coerce into RR instead of clearing it completely misses the point.

I refuse to believe such a large numerical instabilities. In fact maple gives the same answer as pari. I will investigate further and figure out. You cannot Willy nilly ignore the imaginary part.

And Mathematica does the same as pari.

I suggested the try: except: statement because it does not require deep understanding of Bessel's K function and probably in the end we might end up with this solution.

The try except statement does not significantly help the original complaint. Either the identity is true or we should give the patch a negative review and close it as invalid.

besselk(nu,x,{flag = 0})

K-Bessel function of index nu (which can be complex) and argument x. Only real and positive arguments x are allowed in the present version 2.3.3. If flag is equal to 1, uses another implementation of this function which is faster when x >> 1.

The library syntax is kbessel(nu,x,prec) and kbessel2(nu,x,prec) respectively.

Therefore pari gives incorrect answers for your negative x. Perhaps this is another bug and we should merge this ticket and open another one for adding some error checking to besselk?

The complaint is valid. In fact the current version is lot more broken that the patched version. The current version does not take any complex argument. The patched version only fails on negative real axis.

Ok lets patch this then and open a new ticket for the negative real axis case.

I think

if (real(nu) == 0 or imag(nu) == 0) and (imag(z) == 0)

should become

if (real(nu) == 0 or imag(nu) == 0) and (imag(z) == 0 and real(z) > 0)

Following the philosophy that wrong answers are worse than errors, this patch should not go in as is. Probably the code in Rishi's most recent comment is all that is needed for a fix, as long as we're not missing something else. I can't fix this at exactly this moment, but the fix should be trivial. Anyway, even though I posted the original patch, I have to give this a -1.

rishi's code does not prevent brokenness at all (in fact it is 100% equivalent to attempting to trying to return RR(answer, prec). The patch as is makes the answer "more correct," and then we can go back and write code (that makes use of this patch) to make it 100% correct. Alternatively, if someone wants to make a new patch that checks for real(z)>=0 in all cases and throws an error otherwise, I would give that a positive review. However, the modification of real(z)>0 is not sufficient to ensure correctness.

This ticket has been sitting around for a while without any movement. Change the title so that the reports pick up this ticket correctly.

Cheers,

Michael

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.

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.

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.

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.

apply instead of the previous patch

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

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.