16 views

Skip to first unread message

Feb 18, 2010, 1:25:25 AM2/18/10

to sage-support

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:

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

Feb 18, 2010, 4:26:53 AM2/18/10

to sage-s...@googlegroups.com

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:

>

[...]> 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

-------

Feb 18, 2010, 7:02:16 AM2/18/10

to sage-support

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

Feb 18, 2010, 7:27:00 AM2/18/10

to sage-s...@googlegroups.com

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

Feb 18, 2010, 8:27:02 AM2/18/10

to sage-s...@googlegroups.com

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)?

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

Feb 18, 2010, 1:59:30 PM2/18/10

to sage-support

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

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu