pythran slower than scipy for convolution of long doubles

[lesshaste]

I was excited to see that pythran supports long doubles. I wrote the
following pythran code to compute the convolution of two arrays to
test the support.

#pythran export convolve(float128[], float128[])
import numpy as np

def convolve(x: np.float128, y: np.float128) -> np.float128:
    """
    Perform the convolution of two 1D arrays of type np.float128.

    Parameters:
    x (np.float128[]): First input array.
    y (np.float128[]): Second input array.

    Returns:
    np.float128[]: Convolved output array.
    """
    n = len(x)
    m = len(y)
    result_size = n + m - 1
    result = np.zeros(result_size, dtype=np.float128)

    for i in range(result_size):
        start = max(0, i - m + 1)
        end = min(i + 1, n)
        for j in range(start, end):
            result[i] += x[j] * y[i - j]

    return result

I then tested with the array

a = np.array(['1.e+401', '1.e+000', '1.e+401', '1.e+000']*10000,
dtype=np.float128)

from convolution import convolve
%timeit convolve(a, a)
5.49 s ± 75.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

scipy is abandoning its support for long doubles but for the moment is
still gives the correct answer if you use method='direct'

from scipy.signal import convolve as conv_scipy
%timeit conv_scipy(a, a, method='direct')
3.28 s ± 33.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

I was wondering what scipy is doing to get the speedup over pythran?

[lesshaste]

I tested it with float64s and the gap is still large.

#pythran export convolve(float64[], float64[])
import numpy as np

def convolve(x: np.float64, y: np.float64) -> np.float64:
    """
    Perform the convolution of two 1D arrays of type np.float64.

    Parameters:
    x (np.float64[]): First input array.
    y (np.float64[]): Second input array.

    Returns:
    np.float64[]: Convolved output array.
    """
    n = len(x)
    m = len(y)
    result_size = n + m - 1
    result = np.zeros(result_size, dtype=np.float64)

    for i in range(result_size):
        start = max(0, i - m + 1)
        end = min(i + 1, n)
        for j in range(start, end):
            result[i] += x[j] * y[i - j]

    return result

and then

a = np.array(['1.e+301', '1.e+000', '1.e+301', '1.e+000']*10000, dtype=np.float64)
In [3]: a = np.array(['1.e+301', '1.e+000', '1.e+301', '1.e+000']*10000, dtype=np.float64)

In [4]: %timeit convolve(a, a)
1.24 s ± 9.87 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

whereas with numpy it is very fast.

In [6]: %timeit np.convolve(a, a)
143 ms ± 3.58 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)