Optimal Power of N for FFTW

Specific Domains Signal and Image Processing
1

FFTW is fastest when it has small prime factors according to its manual. If there is the option to pad the input data, then with radix 2 FFT algorithms this is done with a “next power of 2” type function. However if the length is just over a power of 2, at large sizes this can lead to arrays of almost double the length.

I was wondering if anyone has experimented with a “next power of small prime” type of function? In this we would select a range of acceptable primes and then figure out the optimum based on some factor of memory and FFT time.

I have created an example of such a function where it selects the least amount of padding. Comments, suggestions and improvements are welcome.

using Primes

function nextn(n; maxprime = max(2, ceil(Int, n/10)))
    nextn = typemax(Int)            # initialize variable 
    for i in primes(maxprime)
        tmp = nextpow(i, n)
        if tmp < nextn
            nextn = tmp
        end
    end
    return nextn
end
3

Related: https://rosettacode.org/wiki/N-smooth_numbers#Julia

It might make sense to add a nextsmooth function to Primes.jl, analogous to nextprime, which returns the next n-smooth number for a given n.

2 Likes
4

Sounds like a good idea. Anyone want to make the PR? That said, it might be easier after I push more of the improvements to prime sieving that I’ve started.

1 Like
5

Is there any reason not to just use nextprod from base Julia?

julia> nextprod((2, 3, 5, 7), 257)
270
3 Likes
6

I’ve found 5-smooth numbers give a nice, small-increment set of N’s for FFTW.

1 Like
7

Right, I forgot that this function was still in Base (it was added in 2012 by @tim.holy when fft was still in Base. I think it was originally slated to be moved out of Base, but I guess it was overlooked in Julia 1.0?)

It might still be nice to have a higher-level nextsmooth(n, x) function in Primes.jl, which is where I would have expected to find this functionality, though maybe it can call nextprod under the hood (unless there is a better algorithm for large n?). (For small n it can just call nextprod with a lookup table of the corresponding primes.)

1 Like
8

Linking related thread.

1 Like
9

Could this instead be exported by FFTW.jl or AbstractFFTs.jl as a utility function with a sensible set of small primes as a default argument?

1 Like
10

An FFTW.nextfftlen(n) function would be reasonable, along with FFTW.nextrfftlen(n) for real-input FFTs (which favor even lengths). But I think this would be in addition to an N-smooth nextsmooth function in Primes.jl.

1 Like
11

DSP.jl has nextfastfft:

    nextfastfft(n)
Return the closest product of 2, 3, 5, and 7 greater than or equal
to `n`. FFTW contains optimized kernels for these sizes and
computes Fourier transforms of input that is a product of these
sizes faster than for input of other sizes.
2 Likes
12

Thanks all for the great information. As I have never made a pull request, I put in a request for nextsmooth function in the issues section of Primes.jl.