Adding an FFT to be able to multiply two polynomials in time O(n log n).

1 files changed, 32 insertions, 0 deletions

diff --git a/math/fft.cpp b/math/fft.cpp
new file mode 100644
index 0000000..89018c5
--- /dev/null
+++ b/math/fft.cpp
@@ -0,0 +1,32 @@
+// s.size() muss eine Zweierpotenz sein!
+vector<cplx> fft(const vector<cplx> &a, bool inverse = 0) {
+	int logn = 1, n = a.size();
+	vector<cplx> A(n);
+	while ((1 << logn) < n) logn++;
+	for (int i = 0; i < n; i++) {
+		int j = 0;
+		for (int k = 0; k < logn; k++) j = (j << 1) | ((i >> k) & 1);
+		A[j] = a[i];
+	}
+	for (int s = 2; s <= n; s <<= 1) {
+		double angle = 2 * PI / s * (inverse ? -1 : 1);
+		cplx ws(cos(angle), sin(angle));
+		for (int j = 0; j < n; j+= s) {
+			cplx w = 1;
+			for (int k = 0; k < s / 2; k++) {
+				cplx u = A[j + k], t = A[j + s / 2 + k];
+				A[j + k] = u + w * t;
+				A[j + s / 2 + k] = u - w * t;
+				if (inverse) A[j + k] /= 2,	A[j + s / 2 + k] /= 2;
+				w *= ws;
+			}
+		}
+	}
+	return A;
+}
+
+// Polynome: a_0, a_1, ... & b_0, b_1, ...
+vector<cplx> a = {0,0,0,0,1,2,3,4}, b = {0,0,0,0,2,3,0,1};
+a = fft(a); b = fft(b);
+for (int i = 0; i < (int)a.size(); i++) a[i] *= b[i];
+a = fft(a,1); // a = a * b