#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;
}


struct Polynomial {
    std::vector<double> a;
    Polynomial(std::vector<double> new_a) : a(new_a) {}
    
    Polynomial operator*(Polynomial r) {  // An FFT based polynomial multiplication algorithm
        int power_2 = 1;
        while(power_2 < (int)(a.size() + r.a.size() - 1)) {
            power_2 *= 2;  // Find the lowest power of 2 that can hold the size of the resulting polynomial
        }
        
        // Perform the FFT
        std::vector<std::complex<double> > x_l(power_2, 0.0);
        std::vector<std::complex<double> > x_r(power_2, 0.0);
        std::vector<std::complex<double> > product(power_2, 0.0);
        
        for(int i=0; i<(int)a.size(); i++) {
            x_l[i] = a[i];
        }
        for(int i=0; i<(int)r.a.size(); i++) {
            x_r[i] = r.a[i];
        }
        x_l = fast_fourier_transform(x_l);  // Do FFT, a.k.a find polynomial values at designated positions
        x_r = fast_fourier_transform(x_r);
        for(int i=0; i<power_2; i++) {  // Multiply the polynomial values at these positions
            product[i] = x_l[i] * x_r[i];
        }
        product = fast_fourier_transform(product, true);  // Do the inverse FFT, a.k.a polynomial interpolation
        
        std::vector<double> result_a(a.size() + r.a.size() - 1);
        for(int i=0; i<(int)result_a.size(); i++) {
            result_a[i] = product[i].real() / power_2;
        }
        return result_a;
    }
};


int main() {
    Polynomial x_1({1, 2, 5, 3});  // 3x^3 + 5x^2 + 2x + 1
    Polynomial x_2({7, 8, 5});  // 5x^2 + 8x + 7
    Polynomial result = x_1 * x_2;
    
    return 0;
}