How to use a mpmath function in a compiled (numba, cython) fonction ?

[Eric Ducasse]

Hi all,

I need to intensively use the following function :

import numpy as np
def besselk_ratio( order:np.int16, z:np.complex128, \

                   r:np.float64=1.0, dorder:np.int16=0, \

                   ctx=mpmath.mp) :
    ratio = ctx.besselk( order+dorder, r*z ) \

              / ctx.besselk( order, z )
    return np.complex128(  ratio )

It is not possible to directly use numpy.special.kv (or kve) function because the order can be more than 250 and  abs(besselk( order, z )) is over the max float64 value, since the ratio can be coded in complex128.

Thus, how can I create a ufunc (or a compiled function) from the mpmath.besselk function (or its source)?

Kind regards,

Eric Ducasse

Associate Professor

Mechanical and Engineering Institute of Bordeaux (I2M)

Ecole Nat. Sup. d'Arts et Metiers

eric.d...@u-bordeaux.fr

[Vinzent Steinberg]

Hi Eric,

I'm not sure it's possible to create a ufunc from mpmath functions. Did you try using gmpy as a backend for mpmath [0]? This should improve performance a bit.

If this is not enough, I would probably look into implementing a ufunc using a C library implementing besselk with enough precision, e.g. using Arb [1] (which was also created by mpmath's author).

Best regards,

Vinzent

[0]: https://mpmath.org/doc/current/setup.html#using-gmpy-optional

[1]: https://arblib.org/

[Eric Ducasse]

Dear Vinzent,

Thank you very much for your advice.

The use of Flint (``Arb was merged into FLINT in 2023.'', https://flintlib.org/doc/)  gives computation time divided by 10 with respect to mpmath computation time.

But it remains unfortunately around 100 times the numpy computation time, when the computation is possible.

Kind regards,

Eric