add a python float version of special funtions?

[Ondrej Certik]

Hi,

is it a good or bad idea to also add a pure python float versions of
special functions to mpmath? It should be much faster, wouldn't it?

Would it be difficult (meaningful) to implement automatic switching
between python floats and mpmath mpf? Like you suggested in the issue
recently.

I'd like to use mpmath for anything, that scipy does if I am not
concerned about speed that much. Like the integration routines. Also I
have some runge-kutta ODE integrators,
so I wonder what is your strategy about float vs. mpf.

I guess a good strategy is to just implement it using python floats
but in a way, so that it is dead easy to use mpmath as a drop-in
replacement?

Ondrej

[Fredrik Johansson]

SciPy already provides a truckload of machine precision special
functions, with excellent (fast and robust) implementations. It'd be
hard to top that.

You can pass to SciPy functions to mpmath functions if you wrap them.

For example, let's try integrating a Bessel function in SciPy:

>>> from scipy.special import j0
>>> from scipy.integrate import quadrature
>>>
>>> quadrature(j0, 2, 4, tol=1e-18, maxiter=1000)
Took 15 points.
(-0.40103613373641989, 0.0)

and then with mpmath:

>>> from scipy.special import j0
>>> from mpmath import mpf, quadts
>>>
>>> quadts(lambda x: mpf(j0(float(x))), 2, 4)
mpf('-0.40103613373642022')

As reference, Mathematica gives:

In[1]:= N[Integrate[BesselJ[0,x], {x,2,4}], 25]
Out[1]= -0.4010361337364200872132472

Mpmath can't go higher than ~15 digits, since the SciPy version of j0
is only accurate to this precision, but using mpmath's integrator
seems to avoid rounding errors during the integration, and as a
consequence the result is two digits more accurate than SciPy's
integrator.

It seems you can't pass j0 to quadts directly, because of this:

>>> j0(mpf(2.5))
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: function not supported for these types, and can't coerce safely to su
pported types

What's the solution here? Define a special method in mpf?

Fredrik

[Ondrej Certik]

Why not to implement j0 in mpmath with any precision and correct
rounding? That will solve everything.

Ondrej

[Fredrik Johansson]

On Mon, Mar 10, 2008 at 6:23 PM, Ondrej Certik <ond...@certik.cz> wrote:

> Why not to implement j0 in mpmath with any precision and correct
> rounding? That will solve everything.

Of course. It's just a matter of actually doing it :-)

Fredrik

[Fredrik Johansson]

Ok, there, you can have your Bessel function :)

>>> from mpmath import *
>>> mp.dps = 50
>>> print quadts(j0, 2, 4)
-0.40103613373642008721324715264599746684510802121868

Mathematica output:

In[1]:= NIntegrate[BesselJ[0,x], {x, 2, 4}, WorkingPrecision->53]
Out[1]= -0.4010361337364200872132471526459974668451080212187

In fact, speed is pretty good, even without psyco:

In [1]: from mpmath import *
In [2]: mp.dps = 50
In [3]: %timeit quadts(j0, 2, 4)
10 loops, best of 3: 54.5 ms per loop

In[5]:= Timing[NIntegrate[BesselJ[0,x], {x, 2, 4}, WorkingPrecision->53];]
Out[5]= {0.012998 Second, Null}
In[6]:= Timing[NIntegrate[BesselJ[0,x], {x, 2, 4}, WorkingPrecision->53];]
Out[6]= {0.013998 Second, Null}

So mpmath is only 4x slower than Mathematica at computing this
integral. Actually probably more like 8x, since I'm running
Mathematica remotely on a machine which I think is about half as fast
as mine.

Fredrik

[Ondrej Certik]

Excellent, thanks for implementing it. :)

Ondrej