FFT

//FFT
//nasarouf@cs.ubc.ca
//ref: AC Just a few more triangles! https://open.kattis.com/problems/moretriangles

#include <ccomplex>
#include <complex>
using namespace std;

typedef double ld;
const ld PI = acos(-1.L);
typedef complex<ld> Complex;

//precomputing these saves 50% time
map<pair<int, int>, Complex> wns;
void FFTinit(int N){
    for (int i = 1; i <= N; i *= 2)
        for (int dir = -1; dir <= 1; dir += 2)
            wns[make_pair(i, dir)] = polar((ld)1, dir * 2 * PI / i);
}

// *** Fast Fourier Transform (Recursive), ubc codearchive, heavily modified... ***
//array pointers instead of vectors, calculations done in-place
void rfft(int n, Complex* a, Complex* y, int dir = 1, int stride = 1) {
    if (n == 1) { y[0] = a[0]; return; }
    Complex wn = wns[make_pair(n, dir)], w = 1;
//  Complex wn = polar((ld)1, dir * 2 * PI / n), w = 1; //precomputed instead
    rfft(n / 2, a, y, dir, stride * 2);
    rfft(n / 2, a + stride, y + n / 2, dir, stride * 2);
    for (int k = 0; k < n / 2; k++, w *= wn) {
        Complex y1 = y[k] + w*y[n / 2 + k];
        Complex y2 = y[k] - w*y[n / 2 + k];
        y[k + n / 2] = y2;
        y[k] = y1;
    }
}
//end code archive
void ifft(int n, Complex* a, Complex* y) {
    for (int i = 0; i < n; i++) a[i] /= ld(n);
    rfft(n, a, y, -1);
}

/** example:

    int N = 1 << int(ceil(log2(sz)));

    FFTinit(N);
    vector<Complex> v(N, 0), c(N, 0), tmp(N, 0);

    //convolution
    for (ll i = 0; i < n; i++) tmp[i] = v[i];
    rfft(N, &tmp[0], &c[0]);
    for (ll i = 0; i < N; i++) tmp[i] = c[i] * c[i];
    ifft(N, &tmp[0], &c[0]);

*/