FFT
FFTについて教わったことをまとめる。
として、
を求める為の方法。
物理においてFFTはこういう関数の微分を誤差が小さい方法で求めるときに使ったりする。
x,yを2進数展開したものを
のようにして書く記法をとる。
yの2進数展開をひっくり返した数列を として書く。
今後を と書くことにする。
の整数倍は指数関数の中で意味はなくなるのでxy/Nの小数部分だけを考える.
関数 を と書くことにすると、
は
と ごとに分けて書くことが出来る。
各 について見ると、 から への2点フーリエ変換とみることが出来るので
各桁の2点フーリエ変換を繰り返すことで は求まる。
について具体的に見ると
ここで を用いた。
についての和は
と のパターンはたくさんある。
これを番目まで繰り返すとが得られ、は
として求まる。
実装としては
をひっくり返すところをちゃんとするために
gを
とし,
と変換を繰り返す。
と をそれぞれ2進数とみなして
と書くことにすると、
この変換を , の全ての値について計算すれば良い。
atc001.contest.atcoder.jp
一応ACとれたのであってるはず。微妙だが。
#include <vector> #include <list> #include <map> #include <set> #include <deque> #include <stack> #include <bitset> #include <algorithm> #include <functional> #include <numeric> #include <utility> #include <sstream> #include <iostream> #include <iomanip> #include <cstdio> #include <cmath> #include <cstdlib> #include <cctype> #include <string> #include <cstring> #include <ctime> #include <queue> #include <complex> using namespace std; //conversion //------------------------------------------ inline int toInt(string s) { int v; istringstream sin(s); sin >> v; return v; } template<class T> inline string toString(T x) { ostringstream sout; sout << x; return sout.str(); } //math //------------------------------------------- template<class T> inline T sqr(T x) { return x * x; } //typedef //------------------------------------------ typedef vector<int> VI; typedef vector<VI> VVI; typedef vector<string> VS; typedef pair<int, int> PII; typedef long long LL; //container util //------------------------------------------ #define ALL(a) (a).begin(),(a).end() #define RALL(a) (a).rbegin(), (a).rend() #define PB push_back #define MP make_pair #define SZ(a) int((a).size()) #define EACH(i,c) for(typeof((c).begin()) i=(c).begin(); i!=(c).end(); ++i) #define EXIST(s,e) ((s).find(e)!=(s).end()) #define SORT(c) sort((c).begin(),(c).end()) //repetition //------------------------------------------ #define FOR(i,a,b) for(int i=(a);i<(b);++i) #define REP(i,n) FOR(i,0,n) //constant //-------------------------------------------- const double EPS = 1e-10; const double PI = acos(-1.0); //clear memory #define CLR(a) memset((a), 0 ,sizeof(a)) //debug #define dump(x) cerr << #x << '=' << (x) << endl; #define debug(x) cerr << #x << '=' << (x) << '('<<'L' << __LINE__ << ')' << ' ' << __FILE__ << endl; typedef complex<double> Comp; const Comp img = Comp(0, 1); void FFT(vector<Comp>& f, int n) { int N = 1 << n; vector<Comp> g(N); for (int m = 1; m <= n; m++) {//どこのbitに着目しているか m番目のbitに着目している int max_xt = (1 << (n - m)) - 1;//これはxtの値の取りうる範囲 2^{n-m} int max_yt = (1 << (m - 1)) - 1;//これはytの値の取りうる範囲 2^{m-1} int start_yt = (1 << (n - m + 1)); Comp a, b; //xt ytの取りうるすべての値を列挙 for (int xt = 0; xt <= max_xt; xt++)for (int yt = 0; yt <= max_yt; yt++) { a = f[yt*start_yt/*---yt 2^{n-m+1}---*/ + xt];//x_m==0の場合 b = f[yt*start_yt + (1 << (n - m)) + xt] * exp(PI*img*double(yt) / double(1 << (m - 1)));//x_m==1の場合 g[(yt << (n - m))/*---yt 2^{n-m}---*/ + xt] = a + b;//Y_m==0の場合 g[(1 << (n - 1))/*2^{n-1}*/ + (yt << (n - m))/*yt 2^{n-m}*/ + xt] = a - b;//Y_m==1の場合 } f = g;//次のループはこのgをfとして扱う } } //逆フーリエ変換をする //ζnをinverse(ζn)に変換(複素共益を取る) してFFTし最後にNで割るだけ void iFFT(vector<Comp>& f, int n) { int N = 1 << n; vector<Comp> g(N); for (int m = 1; m <= n; m++) { int max_xt = (1 << (n - m)) - 1; int max_yt = (1 << (m - 1)) - 1; int start_yt = (1 << (n - m + 1)); Comp a, b; for (int xt = 0; xt <= max_xt; xt++)for (int yt = 0; yt <= max_yt; yt++) { a = f[yt*start_yt + xt]; b = f[yt*start_yt + (1 << (n - m)) + xt] * exp(PI*(-img)*double(yt) / double(1 << (m - 1))); g[(yt << (n - m)) + xt] = a + b; g[(1 << (n - 1)) + (yt << (n - m)) + xt] = a - b; } f = g; } REP(i, N) { f[i] /= (double)N; } } #if 0 int main() { const int n = 3; int N = 1 << n; vector<Comp> f(N); const double k = 3.0; REP(x, N) { f[x] = exp(-2 * PI*img*k*double(x) / double(N)); } cout << "f=" << endl; REP(i, N) { cout << f[i] << " "; } cout << endl; auto a = FFT(f, n); cout << "FFT=" << endl; REP(i, N) { cout << a[i] << " "; } cout << endl; } #else int main() { int n; cin >> n; int a = 1; while (true) { if ((1 << a)>2 * n) { break; } a++; } int N = 1 << a; vector<Comp> A(N); vector<Comp> B(N); FOR(i, 0, n) { int x, y; cin >> x >> y; A[i] = Comp(x, 0), B[i] = Comp(y, 0); } FFT(A, a); FFT(B, a); REP(i, N) { A[i] *= B[i]; } iFFT(A, a); cout << 0 << endl; FOR(i, 0, 2 * n - 1) { cout << (int)round(A[i].real()) << endl; } } #endif