Enhancement: Add prev_good_size to calculate fast FFT length smaller than N

[tpl2go]

Reference Issue: see https://github.com/scipy/scipy/pull/15321

What does this implement/fix? pocketfft currently implements a good_size_real and good_size_cmplx to find the smallest fast FFT length greater than a N. These functions are useful for zero-padding signals to a length suitable for fast FFT. However, in some situations it is easier to cut signals instead of zero-padding to reduce them to a fast FFT length. This PR implements prev_good_size_real and prev_good_size_cmplx to find the largest fast FFT length smaller than N.

[tpl2go]

Just in case you would like a test function to verify correctness of the implementation, I have written one here:

#include <assert.h>
#include "pocketfft_hdronly.h"

using namespace pocketfft::detail;

void test_prev_good_size_cmplx() {

    size_t prev_good_size = 1;
    size_t next_good_size = 1;
    size_t n1, n2, next_good_size2;

    util MyUtil;

    for (size_t n = 1; n <= 1000000; n++) {
        n1 = MyUtil.good_size_cmplx(n);
        n2 = MyUtil.prev_good_size_cmplx(n);

        printf("%d %d %d %d %d\n", n, n1, n2, prev_good_size, next_good_size);

        assert(n1 >= n);
        assert(n2 <= n);
        assert(next_good_size == n1);
        assert(prev_good_size == n2);

        if (n == next_good_size) {
            next_good_size2 = MyUtil.good_size_cmplx(n + 1);
            if (next_good_size2 == n + 1) {
                prev_good_size = n + 1;
                next_good_size = n + 1;
            } else {
                prev_good_size = next_good_size;
                next_good_size = next_good_size2;
            }
        }
        if (n + 1 == next_good_size) {
            prev_good_size = next_good_size;
        }
    }
}

void test_prev_good_size_real() {

    size_t prev_good_size = 1;
    size_t next_good_size = 1;
    size_t n1, n2, next_good_size2;

    util MyUtil;

    for (size_t n = 1; n <= 1000000; n++) {
        n1 = MyUtil.good_size_real(n);
        n2 = MyUtil.prev_good_size_real(n);

        printf("%d %d %d %d %d\n", n, n1, n2, prev_good_size, next_good_size);

        assert(n1 >= n);
        assert(n2 <= n);
        assert(next_good_size == n1);
        assert(prev_good_size == n2);

        if (n == next_good_size) {
            next_good_size2 = MyUtil.good_size_real(n + 1);
            if (next_good_size2 == n + 1) {
                prev_good_size = n + 1;
                next_good_size = n + 1;
            } else {
                prev_good_size = next_good_size;
                next_good_size = next_good_size2;
            }
        }
        if (n + 1 == next_good_size) {
            prev_good_size = next_good_size;
        }
    }
}

static void main(){
    test_prev_good_size_cmplx();
    test_prev_good_size_real();
}

[mreineck]

Thanks a lot for the PR. I'm traveling at the moment, but will have a look as soon as possible!

[mreineck]

BTW, do you want a copyright/author line in the header? If so, feel free to add one!

[tpl2go]

Thanks that is a really nice offer! I will send a commit in tomorrow. Will acknowledge Peter as he has contributed as well

[mreineck]

I'll merge now, thanks again!

[tpl2go]

Thanks a lot !