FFTUtil.icc

#ifndef ISCE_FFT_FFTUTIL_ICC
#error "FFTUtil.icc is an implementation detail of FFTUtil.h"
#endif
 
#include <cmath>
 
namespace isce3 { namespace fft {
 
template<typename T, typename std::enable_if<std::is_integral<T>::value>::type *>
inline
T nextPowerOfTwo(T n)
{
    if (n < 0) {
        throw isce3::except::DomainError(ISCE_SRCINFO(), "input must be non-negative");
    }
 
    if (n <= 1) {
        return 1;
    }
 
    T res = 2;
    --n;
    while (n >>= 1) {
        res <<= 1;
    }
 
    return res;
}
 
inline std::int32_t intpow(std::int32_t base, std::int32_t exponent)
{
    std::int32_t x = 1;
    for (std::int32_t i = 0; i < exponent; ++i) {
        x *= base;
    }
    return x;
}
 
// compute smallest m = 2^a * 3^b * 5^c >= n
inline std::int32_t nextFastPower(std::int32_t n)
{
    if (n < 0) {
        throw isce3::except::DomainError(ISCE_SRCINFO(), "input must be non-negative");
    }
    if (n <= 1) {
        return 1;
    }
 
    // Upper bound on each exponent happens when the others are zero, e.g.,
    // 3^max3 >= n and 5^max5 >= n
    // which is easy to compute with logarithms.
    double logn = std::log(n);
    std::int32_t max5 = std::ceil(logn / std::log(5));
    std::int32_t max3 = std::ceil(logn / std::log(3));
    std::int32_t mmin = INT32_MAX;
    for (std::int32_t x5 = 0; x5 <= max5; ++x5) {
        std::int32_t n5 = intpow(5, x5);
        for (std::int32_t x3 = 0; x3 <= max3; ++x3) {
 
            // Promote type since worst case m is roughly 15 * n^2, so we need
            // double the bits of n.
            std::int64_t m = static_cast<std::int64_t>(n5) * intpow(3, x3);
 
            // Go the rest of the way with factors of two.
            while (m < n) {
                m *= 2;
            }
            if (m < mmin) {
                mmin = static_cast<std::int32_t>(m);
            }
        }
    }
    return mmin;
}
 
}}