FastFourierTransform

fft(a) computes hat f(k) = Σ_x a[x] exp(2π i · k x / N) for all k. N must be a power of 2. Useful for convolution: conv(a, b) = c, where c[x] = Σ a[i]b[x-i]. For convolution of complex numbers or more than two vectors: FFT, multiply pointwise, divide by n, reverse(start+1, end), FFT back. Rounding is safe if (Σ a_i² + Σ b_i²)log_2N < 9·10^14 (in practice 10^16; higher for random inputs). Otherwise, use NTT/FFTMod.

Time: O(N log N) with N = |A|+|B| (~ 1s for N=2²2) 35 lines somewhat tested

content/numerical/FastFourierTransform.h — Ludo Pulles, chilli, Simon Lindholm, source: http://neerc.ifmo.ru/trains/toulouse/2017/fft2.pdf (do read, it's excellent) Accuracy bound from http://www.daemonology.net/papers/fft.pdf

typedef complex<double> C;
typedef vector<double> vd;
void fft(vector<C>& a) {
	int n = sz(a), L = 31 - __builtin_clz(n);
	static vector<complex<long double>> R(2, 1);
	static vector<C> rt(2, 1);  // (^ 10% faster if double)
	for (static int k = 2; k < n; k *= 2) {
		R.resize(n); rt.resize(n);
		auto x = polar(1.0L, acos(-1.0L) / k);
		rep(i,k,2*k) rt[i] = R[i] = i&1 ? R[i/2] * x : R[i/2];
	}
	vi rev(n);
	rep(i,0,n) rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
	rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int k = 1; k < n; k *= 2)
		for (int i = 0; i < n; i += 2 * k) rep(j,0,k) {
			C z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled)
			a[i + j + k] = a[i + j] - z;
			a[i + j] += z;
		}
}
vd conv(const vd& a, const vd& b) {
	if (a.empty() || b.empty()) return {};
	vd res(sz(a) + sz(b) - 1);
	int L = 32 - __builtin_clz(sz(res)), n = 1 << L;
	vector<C> in(n), out(n);
	copy(all(a), begin(in));
	rep(i,0,sz(b)) in[i].imag(b[i]);
	fft(in);
	for (C& x : in) x *= x;
	rep(i,0,n) out[i] = in[-i & (n - 1)] - conj(in[i]);
	fft(out);
	rep(i,0,sz(res)) res[i] = imag(out[i]) / (4 * n);
	return res;
}