Dear PDL people,

I work with long ndarrays (EEG recordings) and ran into this. I’m not quite 
sure whether it is a known issue.
PDL::FFT::fft() fails for some lengths that are products of 2, 3 and 5 only. 
Example: 354294 = 2 x 3^11.
$ perl -MPDL -MPDL::FFT -e '$L=354294; $x=sin(sequence($L)*0.37); 
($re,$im)=($x->copy,zeroes($L)); fft($re,$im); 
$q=sum($re*$re+$im*$im)/$L/sum($x*$x); print abs($q-1)<1e-6 ? "ok\n" : "FAILS 
$q\n"'
On macOS (PDL 2.104, and 2.106 tested from the build tree with -Mblib) it prints
  Error: FFTRADIX() - exceeded number of factors
  Error: FFTRADIX() - insufficient memory.
and then the process dies with “Abort trap: 6”.
On Linux (PDL 2.104) it prints the same two messages, but no exception is 
raised and the result is wrong (“FAILS 2.82251463473838e-06”). The same happens 
on another Ubuntu machine with perl 5.38.2 and PDL 2.085.
ifft() fails the same way for the same length (checked on macOS with 2.106).

What I see in lib/PDL/FFT/fftn.c (2.106; line numbers are the same in master 
3e1ec0e):
“#define NFACTOR 11” (line 193), so factor[] has 11 elements.
factorize() (lines 204-261) writes into factor[] without checking the count. 
Square factors (16, and squares of odd primes) are written a second time as the 
mirrored half, so 2 x 3^11 produces 12 entries and writes past the end of the 
array.
The “nFactor > NFACTOR” check (line 406) only runs after factorize() has 
returned.
In lib/PDL/FFT.pd the return value of fftn() (-1) is not checked, either in 
fft() (line 96) or in ifft() (line 139), so no error reaches Perl.
I haven’t confirmed it, but I assume the abort on macOS comes from detection of 
the out-of-bounds write on the stack.
Between 1000 and 2,000,000, 51 of the 493 lengths of the form 2^a 3^b 5^c 
exceed the limit.
The smallest are 157464 (2^3 x 3^9), 262440, 279936, 354294, 393660, 437400, 
466560 and 472392. 157464, 262440 and 279936 fail on PDL 2.106 (macOS), 2.104 
(Linux) and 2.085. 157500 works on 2.106 and 2.085.
This seems distinct from #436 (fixed in 2.085), since it still happens in 
2.085, 2.104 and 2.106.

PDL::FFTW3 transforms the same 354294-point signal correctly (0.203 on macOS, 
0.204 on Linux), so I’m using PDL::FFTW3 for now.
A possible fix would be to check the count inside factorize() before each 
write, and to barf in fft() and ifft() when fftn() returns an error.

Tested environment:
macOS Tahoe 26.6.2, arm64, perl 5.40.4, PDL 2.104 and 2.106
Ubuntu 24.04.3 (Parallels Desktop VM), arm64, perl 5.38.2, PDL 2.104

Best regards, 
Shugo

_______________________________________________
pdl-general mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pdl-general

Reply via email to