Spin-weighted transforms fail N=41

[eschnett]

This code fragment fails for me.

It works fine when I define N to a different value (anything different from 41). For N=41, the result is either wrong or contains nans. Sometimes, this problem only occurs the second time, hence I added a for loop.

I have FastTransforms v0.13.3 installed.

using FastTransforms
using LinearAlgebra

s = 0

N = 41
M = 2N-1

for iter in 1:2

F = ones(ComplexF64, N, M);

P = plan_spinsph2fourier(F, s)
PA = plan_spinsph_analysis(F, s)
C = copy(F);
lmul!(PA, C);
ldiv!(P, C);

P = plan_spinsph2fourier(C, s)
PS = plan_spinsph_synthesis(C, s)
F′ = copy(C);
lmul!(P, F′);
lmul!(PS, F′);

@assert !any(isnan, F′)
@assert maximum(abs.(F′ - F)) ≤ sqrt(eps())

end

[eschnett]

The almost trivial case N=1 fails as well. It always outputs C[1,1] = 0.

[MikaelSlevinsky]

Thanks for these bug reports. Let's start with the N=41 NaN bug. This test suggests the bug is restricted to executing spinsph_analysis and that it happens for more values of N:

using FastTransforms, LinearAlgebra

for N in 2:1000
    M = 2N-1
    F = ones(ComplexF64, N, M);
    PA = plan_spinsph_analysis(F, 0)
    C = copy(F);
    lmul!(PA, C)
    if norm(C[1] - sqrt(2π)) > sqrt(eps())
        println("This is the failure: ", N)
    end
    if N÷100 == N/100
        println("N = ", N)
    end
end

My output is

This is the failure: 25
This is the failure: 61
This is the failure: 85
N = 100
This is the failure: 113
This is the failure: 145
This is the failure: 181
This is the failure: 182
N = 200
This is the failure: 215
This is the failure: 254
This is the failure: 265
N = 300
This is the failure: 303
This is the failure: 313
N = 400
N = 500
N = 600
N = 700
N = 800
N = 900
N = 1000

Changing the FT_FFTW_FLAGS in ftinternal.h to FFTW_ESTIMATE | FFTW_DESTROY_INPUT results in fewer errors in the same loop:

This is the failure: 61
N = 100
This is the failure: 145
This is the failure: 181
N = 200
This is the failure: 265
N = 300
N = 400
This is the failure: 421
This is the failure: 481
N = 500
N = 600
N = 700
N = 800
This is the failure: 841
N = 900
This is the failure: 925
N = 1000

Two possible hypotheses: 1. there's an issue with the use of the variables X and Y, which get casted to and from double *, ft_complex * and fftw_complex *. 2. the bug is in FFTW.

[MikaelSlevinsky]

As far as casting between double * and fftw_complex * is concerned, this is in use with sph_analysis, which doesn't seem to give errors with the same type of test.

[MikaelSlevinsky]

The issues seems to be fixed when building from source with the latest commit MikaelSlevinsky/FastTransforms@29d9e87. (I guess the real issue is that I don't know how to use FFTW in C.)

I think I'll go ahead and tag a new release.

[eschnett]

Thank you!