#define _USE_MATH_DEFINES
#include <complex>
#include <vector>
#include <cmath>


// Assumes the size of x is a power of 2
std::vector<std::complex<double> > fast_fourier_transform(std::vector<std::complex<double> > x, bool inverse = false) {
    
    std::vector<std::complex<double> > w(x.size(), 0.0);  // w[k] = e^{-i2\pi k/N}
    // Precalculate w[k] for faster FFT computation, do it in a 'binary-search' way to provide decent numeric accuracy
    w[0] = 1.0;
    for(int pow_2 = 1; pow_2 < (int)x.size(); pow_2 *= 2) {
        w[pow_2] = std::polar(1.0, 2*M_PI * pow_2/x.size() * (inverse ? 1 : -1) );
    }
    for(int i=3, last=2; i < (int)x.size(); i++) {
        // This way of computing w[k] guarantees that each w[k] is computed in at most log_2 N multiplications
        if(w[i] == 0.0) {
            w[i] = w[last] * w[i-last];
        } else {
            last = i;
        }
    }
    
    for(int block_size = x.size(); block_size > 1; block_size /= 2) {
        // Do the rearrangement for 'recursive' call block by block for each level'
        std::vector<std::complex<double> > new_x(x.size());
        
        for(int start = 0; start < (int)x.size(); start += block_size) {
            for(int i=0; i<block_size; i++) {
                new_x[start + block_size/2 * (i%2) + i/2] = x[start + i];
            }
        }
        x = new_x;
    }
    
    for(int block_size = 2; block_size <= (int)x.size(); block_size *= 2) {
        // Now compute the FFT 'recursively' level by level bottom-up
        std::vector<std::complex<double> > new_x(x.size());
        int w_base_i = x.size() / block_size;  // w[w_base_i] is the e^{-i2\pi / N} value for the N of this level
        
        for(int start = 0; start < (int)x.size(); start += block_size) {
            for(int i=0; i < block_size/2; i++) {
                new_x[start+i]              = x[start+i] + w[w_base_i*i] * x[start + block_size/2 + i];
                new_x[start+block_size/2+i] = x[start+i] - w[w_base_i*i] * x[start + block_size/2 + i];
            }
        }
        x = new_x;
    }
    return x;
}


int main() {
    std::vector<std::complex<double> > result = fast_fourier_transform({1, 4, 3, 6, 2, 3, 7, 4});
    return 0;
}