comparing cython on GSL library and scipy

[xiaoweiz]

I just did a naive comparison between the Bessel functions in the C-library GSL and their counterparts in Scipy as follows. For a large array, evaluate the Bessel function at each element and sum them up. 

Here's the IPython snippet for the version of Cython+GSL

%%cython -lgsl -lgslcblas

cdef extern from "gsl/gsl_sf_bessel.h":

    double gsl_sf_bessel_J0(double)

def test(double[::1] x):

    cdef int i, N

    N = len(x)

    cdef double s

    s = 0

    for i in range(N):

        s += gsl_sf_bessel_J0(x[i])

    return s

x = numpy.random.randn(1e6) so 

Then %timeit test(x) gives me 

10 loops, best of 3: 61.1 ms per loop

whereas calling scipy.special.j0(x) gives me 

10 loops, best of 3: 20.5 ms per loop

However, if I revise the above comparison by changing the evaluation of Bessel function to sin, then the Cython version is slightly faster. 
%%cython 

def test_2(double[::1] x):
    cdef int i, N
    N = len(x)
    cdef double s
    s = 0
    for i in range(N):
        s += numpy.sin(x[i])
    return s

%timeit test_2(x) gives me 
100 loops, best of 3: 14.3 ms per loop

whereas %timeit numpy.sin(x).sum() gives

100 loops, best of 3: 16.2 ms per loop

Can anyone share some insights on this?

[Robert Bradshaw]

Likely scipy is not using the same algorithm as GSL to compute the
Bessel functions; it's probably using a different, faster (more
approximate?) algorithm.

- Robert

[xiaoweiz]

Yes. That'd be my guess too. I later tried the pure C code version and it has roughly the same performance as the Cython version... 

在 2013年11月27日星期三UTC+8下午3时01分45秒，Robert Bradshaw写道：