Behavior of a PARI function when called within SAGE

[eliot brenner]

I have a program that includes the command to evaluate the K-Bessel
function using Pari: it will sometimes be called with arguments like
the following:

bessel_K(5*i,320,"pari",100)

which evaluates the Bessel function at index (first variable) 5i,
value (second variable) 320, with 100 precision.

This exhibits some unexpected behavior when called from within SAGE,
as follows. The bessel_K fucntion, with fixed index, and growing
second variable exhibits exhibits exponential decay, as shown by these
(correct) calculations performed in Pari itself:

? besselk(5*I,10)

%21 = 0.0000052781217651491219933022050407280366666 + 0.E-63*I

? besselk(5*I,20)

%22 = 0.00000000031100590842180056295117848690804514567 +
1.0123223351797578055 E-56*I

? besselk(5*I,50)

%23 = 2.6618248851542254648076325142970767971 E-23 + 0.E-99*I

? besselk(5*I,100)

%24 = 4.1118977682833709016799697590620762478 E-45 + 0.E-102*I

? besselk(5*I,200)

%25 = 1.1515969255360280546036312473791321981 E-88 +
1.7485172637427202649 E-133*I

(the imaginary parts in the answer are really 0 and can easily be
removed).

however, the same calculations performed with the SAGE interface to
PARI yield:

sage: bessel_K(5*i,10,"pari",100)

5.2781217651491219933022050407e-6

sage: bessel_K(5*i,20,"pari",100)

3.1100590842180056295117848685e-10

sage: bessel_K(5*i,50,"pari",100)

2.6281390403939624075513589787e-23

sage: bessel_K(5*i,100,"pari",100)

-0.029088364960963220354623866640 - 0.040723710945348508499861545085*I

sage: bessel_K(5*i,200,"pari",100)

-4.3246147682420463563103557554e40 +
4.3246147682420463567825924037e40*I

So they agree until the second variable reaches at least 50, but then
SAGE starts reporting wrong answers.

Apparently, SAGE is interpreting very small floating point numbers as
very large ones, "wrapping around" the exponents. What I would like to
do is (obviously) have SAGE report the real answer, or failing that,
have it instead judge a very small floating point number to be 0 and
return 0 as the answer instead of a very large floating point.

Thanks,
Eliot

[Dan Drake]

On Wed, 17 Feb 2010 at 10:25PM -0800, eliot brenner wrote:
> I have a program that includes the command to evaluate the K-Bessel
> function using Pari: it will sometimes be called with arguments like
> the following:
>
> bessel_K(5*i,320,"pari",100)
>
> which evaluates the Bessel function at index (first variable) 5i,
> value (second variable) 320, with 100 precision.
>
> This exhibits some unexpected behavior when called from within SAGE,
> as follows. The bessel_K fucntion, with fixed index, and growing
> second variable exhibits exhibits exponential decay, as shown by these
> (correct) calculations performed in Pari itself:
>

[...]

> ? besselk(5*I,200)
>
> %25 = 1.1515969255360280546036312473791321981 E-88 +
> 1.7485172637427202649 E-133*I

[...]

> sage: bessel_K(5*i,200,"pari",100)
>
> -4.3246147682420463563103557554e40 +
> 4.3246147682420463567825924037e40*I
>
> So they agree until the second variable reaches at least 50, but then
> SAGE starts reporting wrong answers.

Yikes. That's awful.

There are some other ways in which that function is broken; see
http://trac.sagemath.org/sage_trac/ticket/3426.

I'll mention these problems on that ticket, so that when we fix it, we
put in some doctests to make sure this doesn't happen again.

Dan

--
--- Dan Drake
----- http://mathsci.kaist.ac.kr/~drake
-------

[eliot brenner]

Thanks!

In the meantime, any suggestions for a work-around (besides write the
entire program in PARI or C)?

Eliot

On Feb 18, 3:26 am, Dan Drake <dr...@kaist.edu> wrote:
> On Wed, 17 Feb 2010 at 10:25PM -0800, eliot brenner wrote:
> >  I have a program that includes the command to evaluate the K-Bessel
> > function using Pari: it will sometimes be called with arguments like
> > the following:
>
> > bessel_K(5*i,320,"pari",100)
>
> > which evaluates the Bessel function at index (first variable) 5i,
> > value (second variable) 320, with 100 precision.
>
> > This exhibits some unexpected behavior when called from within SAGE,
> > as follows. The bessel_K fucntion, with fixed index, and growing
> > second variable exhibits exhibits exponential decay, as shown by these
> > (correct) calculations performed in Pari itself:
>
> [...]
> > ? besselk(5*I,200)
>
> > %25 = 1.1515969255360280546036312473791321981 E-88 +
> > 1.7485172637427202649 E-133*I
> [...]
> > sage: bessel_K(5*i,200,"pari",100)
>
> > -4.3246147682420463563103557554e40 +
> > 4.3246147682420463567825924037e40*I
>
> > So they agree until the second variable reaches at least 50, but then
> > SAGE starts reporting wrong answers.
>
> Yikes. That's awful.
>

> There are some other ways in which that function is broken; seehttp://trac.sagemath.org/sage_trac/ticket/3426.

>
> I'll mention these problems on that ticket, so that when we fix it, we
> put in some doctests to make sure this doesn't happen again.
>
> Dan
>
> --
> ---  Dan Drake
> -----  http://mathsci.kaist.ac.kr/~drake
> -------
>

>  signature.asc
> < 1KViewDownload

[Burcin Erocal]

Hi Eliot,

On Thu, 18 Feb 2010 04:02:16 -0800 (PST)
eliot brenner <eliotp...@gmail.com> wrote:

> Thanks!
>
> In the meantime, any suggestions for a work-around (besides write the
> entire program in PARI or C)?

You can just use mpmath to compute besselk. The ticket Dan mentioned
(#3426) has these examples:

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)

BTW, does anyone think we need to open a ticket for the error in
converting back from pari. This (From Eliot's previous message) sounds
serious:

> Apparently, SAGE is interpreting very small floating point numbers as
> very large ones, "wrapping around" the exponents. What I would like to
> do is (obviously) have SAGE report the real answer, or failing that,
> have it instead judge a very small floating point number to be 0 and
> return 0 as the answer instead of a very large floating point.

Cheers,
Burcin

[Dan Drake]

On Thu, 18 Feb 2010 at 04:02AM -0800, eliot brenner wrote:
> In the meantime, any suggestions for a work-around (besides write the
> entire program in PARI or C)?

I would use mpmath, which (for the one example I tried) seems to be
agreeing with Pari:

sage: import mpmath
sage: mpmath.besselk(5*I, 200)
mpc(real='1.151596925536028e-88', imag='0.0')
sage: mpmath.mp.dps = 50
sage: mpmath.besselk(5*I, 200)
mpc(real='1.1515969255360280546036312473791321981420160724815173e-88', imag='0.0')

Try mpmath.besselk? for more info, and see the mpmath documentation on
the web -- it's quite nice.

[eliot brenner]

For now the following seems sufficient for what I need to do:

sage: import sage.libs.mpmath.all as mpmath
sage: w=mpmath.call(mpmath.besselk,5*I,200,prec=100)
sage: w
1.1515969255360280546036312474e-88
sage: type(w)
<type 'sage.rings.complex_number.ComplexNumber'>

thanks for your help.
Eliot

>  signature.asc
> < 1KViewDownload