signal.correlate with method='fft' doesn't benefit from long doubles (but is expected to)

[figiel]

test_signaltools.py tests signal.correlate with long double types (both real and complex), which leads to failures on ARMv8. The failures happen only when the 'fft' method is in use and for types larger than 64-bit floats (np.longdouble). The failures are most likely caused by losing precision on conversions of input arrays to double in numpy's fftpack: https://github.com/numpy/numpy/blob/464f79eb1d05bf938d16b49da1c39a4e02506fa3/numpy/fft/fftpack_litemodule.c#L147.

The failure of test for np.clongdouble is related to the precision being lost on prior calculations on .real and .imag parts of complex number, both being long doubles, since correlate using 'auto' method on np.clongdouble doesn't select 'fft'.

The issue is concealed when the tests are running on intel, since the tests require much lesser precision there. This comes from the dynamic requirement regarding the precision implemented at

scipy/scipy/signal/tests/test_signaltools.py

Line 1509 in c2ec858

def equal_tolerance(self, res_dt):

and

scipy/scipy/signal/tests/test_signaltools.py

Line 1662 in c2ec858

def decimal(self, dt):

and also from the fact that the intel double floats offer only 80-bit precision which is reflected in numpy's datatype metadata:

Intel:

>>> np.finfo(np.longdouble)
finfo(resolution=1e-18, min=-1.189731495357231765e+4932, max=1.189731495357231765e+4932, dtype=float128)

ARMv8 (and likely all other architectures?):

>>> np.finfo(np.longdouble)
finfo(resolution=1e-33, min=-1.189731495357231765085759326628007e+4932, max=1.189731495357231765085759326628007e+4932, dtype=float128)

IMO either the tests' requirements towards precision for long doubles could be amended to resolve that, FFT implementation could be extended to support wider types or the tests for long doubles could be removed.

Reproducing code example:

Reproducible on embedded ARM v8 device running cross compiled scipy but also on Ubuntu aarch64 ran in qemu with native compiled scipy (I put some hints how to get it running on PC under #9515. To reproduce just run tests from test_signaltools.py.

Error message:

Test failures from amended test_signaltools.py to only run failing tests:

test_signaltools.py FF

======================================================================================= FAILURES =======================================================================================
_______________________________________________________________________ TestCorrelateReal.test_method[float128] ________________________________________________________________________

self = <test_signaltools.TestCorrelateReal object at 0x400fcacda0>, dt = <class 'numpy.float128'>

    def test_method(self, dt):
        if dt == Decimal:
            method = choose_conv_method([Decimal(4)], [Decimal(3)])
            assert_equal(method, 'direct')
        else:
            a, b, y_r = self._setup_rank3(dt)
            y_fft = correlate(a, b, method='fft')
            y_direct = correlate(a, b, method='direct')
    
>           assert_array_almost_equal(y_r, y_fft, decimal=self.equal_tolerance(y_fft.dtype))
E           AssertionError: 
E           Arrays are not almost equal to 16 decimals
E           
E           (mismatch 73.61111111111111%)
E            x: array([   0.,  184.,  504.,  912., 1360.,  888.,  472.,  160.,   46.,
E                   432., 1062., 1840., 2672., 1698.,  864.,  266.,  134.,  736.,
E                  1662., 2768., 3920., 2418., 1168.,  314.,  260.,  952., 1932.,...
E            y: array([2.2737367544323206e-13, 1.8400000000000017e+02,
E                  5.0399999999999994e+02, 9.1200000000000000e+02,
E                  1.3600000000000002e+03, 8.8800000000000023e+02,...

test_signaltools.py:67: AssertionError
_____________________________________________________________________ TestCorrelateComplex.test_rank3[complex256] ______________________________________________________________________

self = <test_signaltools.TestCorrelateComplex object at 0x401130a0f0>, dt = <class 'numpy.complex256'>

    def test_rank3(self, dt):
        a = np.random.randn(10, 8, 6).astype(dt)
        a += 1j * np.random.randn(10, 8, 6).astype(dt)
        b = np.random.randn(8, 6, 4).astype(dt)
        b += 1j * np.random.randn(8, 6, 4).astype(dt)
    
        y_r = (correlate(a.real, b.real)
               + correlate(a.imag, b.imag)).astype(dt)
        y_r += 1j * (-correlate(a.real, b.imag) + correlate(a.imag, b.real))
    
        y = correlate(a, b, 'full')
>       assert_array_almost_equal(y, y_r, decimal=self.decimal(dt) - 1)
E       AssertionError: 
E       Arrays are not almost equal to 21 decimals
E       
E       (mismatch 100.0%)
E        x: array([ 1.49073807792760147216 -0.396757518788941577036j,
E              -1.449283073249400990152-2.05675954843281385982j ,
E               1.845679171573719308005+2.420984271222354302969j, ...,...
E        y: array([ 1.490738077927598492778-0.396757518788940277688j,
E              -1.449283073249399977911-2.056759548432810613683j,
E               1.84567917157372068937 +2.420984271222352934849j, ...,...

test_signaltools.py:138: AssertionError
=============================================================================== 2 failed in 3.57 seconds ===============================================================================

Scipy/Numpy/Python version information:

1.3.0.dev0+b635563 1.15.4 sys.version_info(major=3, minor=5, micro=2, releaselevel='final', serial=0)

Also reproducible on v1.1.0.

[rgommers]

IMO either the tests' requirements towards precision for long doubles could be amended to resolve that, FFT implementation could be extended to support wider types or the tests for long doubles could be removed.

The test is wrong here. The intent is just to verify that longdouble works, not that it is more precise than double. Also none of the fft implementations in numpy or scipy support longdouble (it may work, but internally there will be downcasting to float64.