FFT初步学习小结
来源：互联网 发布：看现货黄金数据的网址 编辑：程序博客网 时间：2024/05/02 02:35
FFT其实没什么需要特别了解的，了解下原理，（特别推荐算法导论上面的讲解），模板理解就行了。重在运用吧。
处理过程中要特别注意精度。
先上个练习的地址吧：
http://vjudge.net/vjudge/contest/view.action?cid=53596#overview
Problem A A * B Problem Plus
A*B的大数乘法，似乎大数模拟乘法不行的，得用FFT优化到nlogn，将一个数AnAn-1....A1A0,看做An*10^n+An-1*10^n-1+....A1*10+A0*10^0,这样就可以将两个数相乘当成多项式乘法了。
代码：
1 #include <cstdio> 2 #include <iostream> 3 #include <cstring> 4 #include <string> 5 #include <cstdlib> 6 #include <algorithm> 7 #include <vector> 8 #include <cmath> 9 #include <queue> 10 #include <map> 11 #include <set> 12 #include <complex> 13 #define pb push_back 14 #define in freopen("solve_in.txt", "r", stdin); 15 #define out freopen("solve_out.txt", "w", stdout); 16 #define pi (acos(-1.0)) 17 #define bug(x) printf("line %d :>>>>>>>>>>>>>>>\n", x); 18 #define pb push_back 19 using namespace std; 20 #define esp 1e-8 21 22 typedef complex<double> CD; 23 const int maxn = 50000+100; 24 char s1[maxn], s2[maxn]; 25 26 struct Complex{ 27 double x, y; 28 Complex(){} 29 Complex(double x, double y):x(x), y(y){} 30 }; 31 Complex operator + (const Complex &a, const Complex &b){ 32 Complex c; 33 c.x = a.x+b.x; 34 c.y = a.y+b.y; 35 return c; 36 } 37 Complex operator - (const Complex &a, const Complex &b){ 38 Complex c; 39 c.x = a.x-b.x; 40 c.y = a.y-b.y; 41 return c; 42 } 43 Complex operator * (const Complex &a, const Complex &b){ 44 Complex c; 45 c.x = a.x*b.x-a.y*b.y; 46 c.y = a.x*b.y+a.y*b.x; 47 return c; 48 } 49 50 inline void FFT(vector<Complex> &a, bool inverse) 51 { 52 int n = a.size(); 53 for(int i = 0, j = 0; i < n; i++) 54 { 55 if(j > i) 56 swap(a[i], a[j]); 57 int k = n; 58 while(j & (k>>=1)) j &= ~k; 59 j |= k; 60 } 61 double PI = inverse ? -pi : pi; 62 for(int step = 2; step <= n; step <<= 1) 63 { 64 double alpha = 2*PI/step; 65 Complex wn(cos(alpha), sin(alpha)); 66 for(int k = 0; k < n; k += step) 67 { 68 Complex w(1, 0); 69 for(int Ek = k; Ek < k+step/2; Ek++) 70 { 71 int Ok = Ek + step/2; 72 Complex u = a[Ek]; 73 Complex t = a[Ok]*w; 74 a[Ok] = u-t; 75 a[Ek] = u+t; 76 w = w*wn; 77 } 78 } 79 } 80 if(inverse) 81 for(int i = 0; i < n; i++) 82 a[i].x = (a[i].x/n); 83 } 84 vector<double> operator * (const vector<double> &v1, const vector<double> &v2) 85 { 86 int S1 = v1.size(), S2 = v2.size(); 87 int S = 2; 88 while(S < S1+S2) S <<= 1; 89 vector<Complex> a(S), b(S); 90 for(int i = 0; i < S; i++) 91 a[i].x = a[i].y = b[i].x = b[i].y = 0.0; 92 for(int i = 0; i < S1; i++) 93 a[i].x = v1[i]; 94 for(int i = 0; i < S2; i++) 95 b[i].x = v2[i]; 96 FFT(a, false); 97 FFT(b, false); 98 for(int i = 0; i < S; i++) 99 a[i] = a[i] * b[i];100 FFT(a, true);101 vector<double> res(S1+S2-1, 0.0);102 for(int i = 0; i < S1+S2-1; i++)103 res[i] = a[i].x;104 return res;105 }106 int ans[maxn*2+100];107 vector<double > v1, v2;108 int main()109 {110 111 while(scanf("%s%s", s1, s2) != -1)112 {113 114 v1.clear();115 v2.clear();116 int len1 = strlen(s1);117 int len2 = strlen(s2);118 for(int i = 0; s1[i]; i++)119 v1.pb((double)(s1[len1-1-i]-'0'));120 for(int i = 0; s2[i]; i++)121 v2.pb((double)(s2[len2-1-i]-'0'));122 v1 = v1*v2;123 memset(ans, 0, sizeof(ans));124 int carry = 0, top = 0;125 for(int i = 0; i < v1.size(); i++){126 carry += (int)(v1[i]+0.5);127 ans[top++] = carry%10;128 carry /= 10;129 }130 while(carry){131 ans[top++] = carry%10;132 carry /= 10;133 }134 while(top)135 {136 if(ans[top])137 break;138 top--;139 }140 for(int i = top; i >= 0; i--)141 printf("%d", ans[i]);142 cout<<endl;143 }144 return 0;145 }
Problem B 3-idiots
题意：给出一系列边长，问从中选出3条边构成三角形的概率。
分析：首先求出任意两边之和相同的取法有多少种，用每种长度边的数目作系数，FFT即可，然后就是去重了。
对所有边长排序，对于第i条边，sum[i]表示两边之和大于a[i]的边的对数，减去特殊情况，
1.两边均大于a[i],
2，两边一个大于a[i]，
3，两边中包含a[i]。
代码：
1 #include <cstdio> 2 #include <iostream> 3 #include <cstring> 4 #include <string> 5 #include <cstdlib> 6 #include <algorithm> 7 #include <vector> 8 #include <cmath> 9 #include <queue> 10 #include <map> 11 #include <set> 12 #include <complex> 13 #define pb push_back 14 #define in freopen("solve_in.txt", "r", stdin); 15 #define out freopen("solve_out.txt", "w", stdout); 16 #define pi (acos(-1.0)) 17 #define bug(x) printf("line %d :>>>>>>>>>>>>>>>\n", x); 18 #define pb push_back 19 using namespace std; 20 #define esp 1e-8 21 22 typedef complex<double> CD; 23 const int maxn = 100000+100; 24 struct Complex 25 { 26 double x, y; 27 Complex() {} 28 Complex(double x, double y):x(x), y(y) {} 29 }; 30 Complex operator + (const Complex &a, const Complex &b) 31 { 32 Complex c; 33 c.x = a.x+b.x; 34 c.y = a.y+b.y; 35 return c; 36 } 37 Complex operator - (const Complex &a, const Complex &b) 38 { 39 Complex c; 40 c.x = a.x-b.x; 41 c.y = a.y-b.y; 42 return c; 43 } 44 Complex operator * (const Complex &a, const Complex &b) 45 { 46 Complex c; 47 c.x = a.x*b.x-a.y*b.y; 48 c.y = a.x*b.y+a.y*b.x; 49 return c; 50 } 51 inline void FFT(vector<Complex> &a, bool inverse) 52 { 53 int n = a.size(); 54 for(int i = 0, j = 0; i < n; i++) 55 { 56 if(j > i) 57 swap(a[i], a[j]); 58 int k = n; 59 while(j & (k>>=1)) j &= ~k; 60 j |= k; 61 } 62 double PI = inverse ? -pi : pi; 63 for(int step = 2; step <= n; step <<= 1) 64 { 65 double alpha = 2*PI/step; 66 Complex wn = Complex(cos(alpha), sin(alpha)); 67 for(int k = 0; k < n; k += step) 68 { 69 Complex w(1, 0); 70 for(int Ek = k; Ek < k+step/2; Ek++) 71 { 72 int Ok = Ek + step/2; 73 Complex u = a[Ek]; 74 Complex t = w*a[Ok]; 75 a[Ek] = u+t; 76 a[Ok] = u-t; 77 w = wn*w; 78 } 79 80 } 81 } 82 if(inverse) 83 for(int i = 0; i < n; i++) 84 a[i].x = a[i].x/n; 85 } 86 87 vector<double> operator * (const vector<double> &v1, const vector<double> &v2) 88 { 89 int S1 = v1.size(), S2 = v2.size(); 90 int S = 2; 91 while(S < S1+S2) S <<= 1; 92 vector<Complex> a(S, Complex(0.0, 0.0)), b(S, Complex(0.0, 0.0)); 93 for(int i = 0; i < S1; i++) 94 a[i].x = v1[i]; 95 for(int i = 0; i < S2; i++) 96 b[i].x = v2[i]; 97 FFT(a, false); 98 FFT(b, false); 99 for(int i = 0; i < S; i++)100 a[i] = a[i] * b[i];101 FFT(a, true);102 vector<double> res(S1+S2-1, 0.0);103 for(int i = 0; i < S1+S2-1; i++)104 res[i] = a[i].x;105 return res;106 }107 int a[maxn];108 typedef long long LL;109 LL sum[maxn<<1];110 typedef long long LL;111 int n;112 int main()113 {114 in115 int T;116 for(int t = scanf("%d", &T); t <= T; t++)117 {118 LL ans =0;119 scanf("%d", &n);120 int mx = 0;121 for(int i = 0; i < (maxn<<1); i++)122 sum[i] = 0;123 for(int i = 0; i < n; i++)124 {125 scanf("%d", &a[i]);126 sum[a[i]]++;127 mx = max(a[i], mx);128 }129 int tmp = 2;130 while(tmp < 2*(mx+1)) tmp <<= 1;//上界 是mx + 1！131 vector<Complex> v2(tmp, Complex(0.0, 0.0));132 for(int i = 0; i <= mx; i++)133 v2[i] = (Complex((double)sum[i], 0.0));134 FFT(v2, 0);135 136 for(int i = 0; i < v2.size(); i++)137 {138 v2[i] = v2[i]*v2[i];139 }140 FFT(v2, 1);141 for(int i = 0; i <= mx*2; i++)142 {143 sum[i] = (LL)(v2[i].x+0.5);144 }145 for(int i = 0; i < n; i++)146 sum[a[i]*2]--;147 for(int i = 0; i <= 2*mx; i++)148 sum[i] /= 2;149 for(int i = 1; i <= 2*mx; i++)150 sum[i] += sum[i-1];151 sort(a, a+n);152 for(int i = 0; i < n; i++)153 {154 ans += (sum[mx*2]-sum[a[i]]);155 ans -= ((LL)(n-i-1)*(n-i-2)/2 + (n-1) + (LL)i*(n-1-i));//注意用long long保证精度156 }157 printf("%.7f\n", 1.0*ans/((LL)n*(n-1)*(n-2)/6));158 }159 return 0;160 }
Problem C K-neighbor substrings
给定A，B两个01串，求A中和B串长度相同的且哈密顿距离不超过K的不同子串个数。
分析：
求一个串和另一个串的哈密顿距离，也就是对应位置字符不同的个数，长度-同为1-同为0就是结果了。怎么样找同为1或者同为0的个数呢？由于当且仅当同为1，相乘结果为1，所以将B串反转，那么作FFT，对应系数为1表示相应位置同为1，然后统计系数的和就是两串同为1个数。同为0的个数计算也很简单，只需要将B串翻转一下，0变1,1变0，同样FFT就行了。
由于题目要求不同子串的个数，利用Hash就可以了。
代码：
1 #include <cstdio> 2 #include <iostream> 3 #include <cstring> 4 #include <string> 5 #include <cstdlib> 6 #include <algorithm> 7 #include <vector> 8 #include <cmath> 9 #include <queue> 10 #include <map> 11 #include <set> 12 #include <bitset> 13 #include <complex> 14 #define pb push_back 15 #define in freopen("solve_in.txt", "r", stdin); 16 #define out freopen("solve_out.txt", "w", stdout); 17 #define pi (acos(-1.0)) 18 #define bug(x) printf("line %d :>>>>>>>>>>>>>>>\n", x); 19 #define pb push_back 20 #define esp 1e-8 21 using namespace std; 22 23 typedef long long LL; 24 typedef unsigned long long ULL; 25 typedef map<ULL, int> MPLL; 26 const int maxn = 100000+100; 27 const int bb = 100009; 28 unsigned long long Hash[maxn], sum[maxn]; 29 MPLL mps; 30 bitset<maxn> b[2]; 31 int n, m, K; 32 struct Complex 33 { 34 double x, y; 35 Complex() {} 36 Complex(double x, double y):x(x), y(y) {} 37 }; 38 Complex operator + (const Complex &a, const Complex &b) 39 { 40 Complex c; 41 c.x = a.x+b.x; 42 c.y = a.y+b.y; 43 return c; 44 } 45 Complex operator - (const Complex &a, const Complex &b) 46 { 47 Complex c; 48 c.x = a.x-b.x; 49 c.y = a.y-b.y; 50 return c; 51 } 52 Complex operator * (const Complex &a, const Complex &b) 53 { 54 Complex c; 55 c.x = a.x*b.x-a.y*b.y; 56 c.y = a.x*b.y+a.y*b.x; 57 return c; 58 } 59 inline void FFT(vector<Complex> &a, bool inverse) 60 { 61 int n = a.size(); 62 for(int i = 0, j = 0; i < n; i++) 63 { 64 if(j > i) 65 swap(a[i], a[j]); 66 int k = n; 67 while(j & (k>>=1)) j &= ~k; 68 j |= k; 69 } 70 double PI = inverse ? -pi : pi; 71 for(int step = 2; step <= n; step <<= 1) 72 { 73 double alpha = 2*PI/step; 74 Complex wn = Complex(cos(alpha), sin(alpha)); 75 for(int k = 0; k < n; k += step) 76 { 77 Complex w(1, 0); 78 for(int Ek = k; Ek < k+step/2; Ek++) 79 { 80 int Ok = Ek + step/2; 81 Complex u = a[Ek]; 82 Complex t = w*a[Ok]; 83 a[Ek] = u+t; 84 a[Ok] = u-t; 85 w = wn*w; 86 } 87 88 } 89 } 90 if(inverse) 91 for(int i = 0; i < n; i++) 92 a[i].x = a[i].x/n+esp; 93 } 94 int x[maxn]; 95 96 void go(int S1, int S2) 97 { 98 int S = 2; 99 while(S < S1+S2) S <<= 1;100 vector<Complex> sa(S, Complex(0.0, 0.0)), sb(S, Complex(0.0, 0.0));101 for(int i = 0; i < S1; i++)102 sa[i].x = b[0][i];103 for(int i = 0; i < S2; i++)104 sb[i].x = b[1][i];105 FFT(sa, false);106 FFT(sb, false);107 for(int i = 0; i < S; i++)108 sa[i] = sa[i] * sb[i];109 FFT(sa, true);110 for(int i = 0; i <= n-m; i++)111 {112 x[i] += round(sa[i+m-1].x);113 }114 }115 void solve()116 {117 int ans = 0;118 go(n, m);119 b[0].flip();120 b[1].flip();121 go(n, m);122 for(int i = 0; i <= n-m; i++)123 {124 if(m-x[i] <= K)125 {126 int ok = 0;127 for(int k = 0; k < 1; k++)128 {129 ULL tmp = Hash[i]-Hash[i+m]*sum[m];130 if(mps.count(tmp))131 {132 ok++;133 }134 }135 if(ok == 0)136 {137 ans++;138 mps[Hash[i]-Hash[i+m]*sum[m]] = 1;139 // mps[1][Hash[1][0][i]-Hash[1][0][i+m]*sum[1][m]] = 1;140 }141 }142 }143 printf("%d\n", ans);144 }145 char A[maxn], B[maxn];146 void pre()147 {148 sum[0] = 1;149 for(int i = 1; i < maxn; i++)150 sum[i] = sum[i-1]*bb;151 }152 153 154 int main()155 {156 157 pre();158 int kase = 0;159 while(scanf("%d", &K), K >= 0)160 {161 printf("Case %d: ", ++kase);162 mps.clear();163 Hash[n] = 0;164 scanf("%s%s", A, B);165 n = strlen(A), m = strlen(B);166 for(int i = n-1; i >= 0; i--)167 {168 x[i] = 0;169 int t = A[i]-'a';170 b[0][i] = t;171 Hash[i] = Hash[i+1]*bb+t;172 }173 for(int i = m-1; i >= 0; i--)174 {175 int t = B[i]-'a';176 b[1][m-1-i] = t;177 }178 solve();179 }180 return 0;181 }
Problem D Linear recursive sequence
f(k) = 0(k <=0)
f(k) = af(k-p)+bf(k-q)
a,b,k<=10^9, p < q <= 10^4
分析：
具体用到了叉姐的论文《矩阵乘法递推的优化》，其实总结起来就是一个式子，W^i = b1*W^k-1+b2*W^k-2+......bk*E,也就是任何一个W^i都可以表示成E，W^1,W^2,
.....W^k-1的线性组合。求f(n)也就是求出f(n)对应的W^i,(i = n-k+1), 实际就是多次求一个多项式乘法，每次复杂度O(k^2)， 这样可以将矩阵乘法优化到k^2log(n)。
但是本题范围较大即使k^2log(n)也不够优化，而且递推系数只有2个，其他都为0，因此可以用FFT加速多项式乘法，O(qlogqlogn)已经很优化了。
代码：
1 #pragma comment(linker, "/STACK:16777216") 2 #include <cstdio> 3 #include <iostream> 4 #include <cstring> 5 #include <string> 6 #include <cstdlib> 7 #include <algorithm> 8 #include <vector> 9 #include <cmath> 10 #include <queue> 11 #include <map> 12 #include <set> 13 #include <bitset> 14 #include <complex> 15 #define inf 0x0f0f0f0f 16 #define pb push_back 17 #define in freopen("solve_in.txt", "r", stdin); 18 #define out freopen("solve_out.txt", "w", stdout); 19 #define pi (acos(-1.0)) 20 #define bug(x) printf("line %d :>>>>>>>>>>>>>>>\n", x); 21 #define pb push_back 22 #define esp 1e-8 23 typedef long long LL; 24 using namespace std; 25 26 const int M = 119; 27 int n, p, q, a, b; 28 29 struct Complex { 30 double x, y; 31 Complex() {} 32 Complex(double x, double y):x(x), y(y) {} 33 Complex operator + (const Complex &o)const { 34 return Complex(x+o.x, y+o.y); 35 } 36 Complex operator - (const Complex &o)const { 37 return Complex(x-o.x, y-o.y); 38 } 39 Complex operator * (const Complex &o)const { 40 return Complex(x*o.x-y*o.y, y*o.x+x*o.y); 41 } 42 }; 43 void add(double &x, double y) { 44 long long a = (LL)(x+.5), b = (LL)(y+.5); 45 a += b; 46 a%=M; 47 x = (double)a; 48 } 49 void FFT(vector<Complex> &a, bool inverse) { 50 int nn = a.size(); 51 for(int i =0, j = 0; i < nn; i++) { 52 if(j > i) 53 swap(a[i], a[j]); 54 int k = nn; 55 while(j &(k >>= 1)) j &=~k; 56 j |= k; 57 } 58 double PI = inverse ? -pi : pi; 59 for(int step = 2; step <= nn; step <<= 1) { 60 Complex wn(cos(PI*2/step), sin(PI*2/step)); 61 for(int j = 0; j < nn; j += step) { 62 Complex w(1, 0); 63 for(int Ek = j; Ek < j+step/2; Ek++) { 64 int Ok = Ek + step/2; 65 Complex u = a[Ek]; 66 Complex t = w*a[Ok]; 67 a[Ek] = u+t; 68 a[Ok] = u-t; 69 w = w*wn; 70 } 71 } 72 } 73 if(inverse) 74 for(int i = 0; i < nn; i++) 75 a[i].x = a[i].x/nn; 76 } 77 vector<double> operator *(const vector<double> &v1, const vector<double> &v2) { 78 int S1 = v1.size(); 79 int S2 = v2.size(); 80 int S = 2; 81 while(S < S1+S2) S <<= 1; 82 vector<Complex> aa(S, Complex(0.0, 0.0)), ab(S, Complex(0.0, 0.0)); 83 for(int i = 0; i < S1; i++) 84 aa[i].x = v1[i]; 85 for(int i = 0; i < S2; i++) 86 ab[i].x = v2[i]; 87 FFT(aa, false); 88 FFT(ab, false); 89 for(int i = 0; i < S; i++) 90 aa[i] = aa[i]*ab[i]; 91 FFT(aa, true); 92 vector<double> res(S1+S2-1, 0.0); 93 for(int i = 0; i < S1+S2-1; i++) { 94 res[i] = aa[i].x; 95 // cout<<aa[i].y<<endl; 96 // cout<<res[i]<<endl; 97 } 98 for(int i = S1+S2-2; i >= q; i--) { 99 add(res[i-p], a*res[i]);100 add(res[i-q], b*res[i]);101 }102 res.resize(q);103 return res;104 }105 const int maxn = (int)3e4+100;106 107 int h[maxn];108 void calPre() {109 memset(h, 0, sizeof h);110 h[0] = 1;111 for(int i = 1; i < 2*q-1; i++) {112 if(i<=p) h[i] = (h[i]+a)%M;113 else h[i] = (h[i]+a*h[i-p])%M;114 if(i-q <= 0) h[i] = (h[i]+b)%M;115 else h[i] = (h[i]+b*h[i-q])%M;116 }117 }118 int main() {119 120 while(scanf("%d%d%d%d%d", &n, &a, &b, &p, &q) == 5) {121 a %= M;122 b %= M;123 calPre();124 if(n < q) {125 cout<<h[n]<<endl;126 continue;127 }128 n=n-q+1;129 vector<double> Ma(q, 0.0), base(q, 0.0);130 Ma[0] = 1.0;131 base[1] = 1.0;132 while(n) {133 if(n&1) {134 Ma = Ma*base;135 }136 base = base*base;137 n>>=1;138 }139 //// Ma = Ma*base;140 //// cout<<base[0]<<base[1]<<endl;141 //// cout<<Ma[0]<<Ma[1]<<endl;142 double res = 0;143 for(int i = 0; i < q; i++) {144 add(res, (Ma[i]*h[q-1+i]));145 }146 cout<<(int)(res+.5)<<endl;147 }148 return 0;149 }
HDU G++超时，C++2000ms，真是无力了。还有连叉姐标程C++都会超时，G++才能过。不过好像极端数据，标程和我的代码都不能再规定时间跑完，大概15s左右？所以说数据还是很水的。
Problem E Cipher Message 3
题意：给定n个8二进制串，和m个8位二进制串构成2个序列，要求在n个串的序列中找到连续的一段，使得和m个二进制串匹配，匹配的意思是两个序列中对应的二进制串前7位相同，后面一位可以不同，但是会有额外的花费，求最后花费最小的一个匹配序列，而且要求该序列在n中尽量靠左。
分析:
训练赛时遇到的一道题，当时不会FFT，甚至不知道这玩意在竞赛里还能干啥，首先利用KMP找出所有能够匹配的位置，当时知道怎么找花费最小的了。利用FFT求出两个串在各个位置开始的哈密顿距离就可以了。转化和上面C题一样。
代码：
1 #include <cstdio> 2 #include <iostream> 3 #include <cstring> 4 #include <string> 5 #include <cstdlib> 6 #include <algorithm> 7 #include <vector> 8 #include <queue> 9 #include <ctime> 10 #include <map> 11 #include <set> 12 #include <cmath> 13 #include <bitset> 14 #define pb push_back 15 #define in freopen("solve_in.txt", "r", stdin); 16 #define out freopen("solve_out.txt", "w", stdout); 17 #define pi (acos(-1.0)) 18 #define bug(x) printf("line %d :>>>>>>>>>>>>>>>\n", x); 19 #define pb push_back 20 using namespace std; 21 #define esp 1e-8 22 23 const int maxn = 250000+100; 24 int sa[maxn], sb[maxn]; 25 int n, m; 26 int x[maxn]; 27 bitset<maxn> da, db; 28 struct Complex 29 { 30 double x, y; 31 Complex() {} 32 Complex(double x, double y):x(x), y(y) {} 33 inline Complex operator + (const Complex &b)const 34 { 35 return Complex(x+b.x, y+b.y); 36 } 37 inline Complex operator - ( const Complex &b)const 38 { 39 return Complex(x-b.x, y-b.y); 40 } 41 inline Complex operator * (const Complex &b)const 42 { 43 return Complex(x*b.x-y*b.y, x*b.y+y*b.x); 44 } 45 }; 46 47 48 inline void FFT(vector<Complex> &a, bool inverse) 49 { 50 int n = a.size(); 51 for(int i = 1, j = n/2; i < n-1; i++) 52 { 53 if(j > i) 54 swap(a[i], a[j]); 55 int k = n>>1; 56 while(j >= k) 57 { 58 j-=k; 59 k >>= 1; 60 } 61 j += k; 62 } 63 double PI = inverse ? -pi : pi; 64 for(int step = 2; step <= n; step <<= 1) 65 { 66 double alpha = 2*PI/step; 67 Complex wn(cos(alpha), sin(alpha)); 68 69 for(int k = 0; k < n; k+=step) 70 { 71 Complex w(1, 0); 72 for(int j = k; j < k+step/2; j++) 73 { 74 75 Complex u = a[j]; 76 Complex t = w*a[j+step/2]; 77 a[j] = u+t; 78 a[j+step/2] = u-t; 79 w = w*wn; 80 } 81 } 82 } 83 if(inverse) 84 for(int i = 0; i < n; i++) 85 a[i].x = a[i].x/n; 86 } 87 void go() 88 { 89 int S1 = n, S2 = m; 90 int tmp = max(S1, S2); 91 int S = 2; 92 while(S < S1+S2) S <<= 1; 93 vector<Complex> a(S, Complex(0.0, 0.0)), b(S, Complex(0.0, 0.0)); 94 for(int i = 0; i < S1; i++) 95 a[i].x = da[i]; 96 for(int i = 0; i < S2; i++) 97 b[i].x = db[i]; 98 FFT(a, false); 99 FFT(b, false);100 for(int i = 0; i < S; i++)101 {102 Complex c;103 c.x = (a[i].x*b[i].x)-(a[i].y*b[i].y);104 c.y = (a[i].x*b[i].y)+(a[i].y*b[i].x);105 a[i] = c;106 // a[i] = a[i]*b[i];107 }108 FFT(a, true);109 for(int i = 0; i <= n-m; i++)110 x[i] += round(esp+a[i+m-1].x);111 }112 int f[maxn];113 int match[maxn];114 int cnt = 0;115 void getFail(int a[], int n)116 {117 f[0] = f[1] = 0;118 int j;119 for(int i = 1; i < n; i++)120 {121 j = f[i];122 while(j && a[i] != a[j]) j = f[j];123 f[i+1] = (a[i] == a[j] ? j+1 : 0);124 }125 }126 void KMP(int T[], int P[], int n, int m)127 {128 getFail(P, m);129 int j = 0;130 for(int i = 0; i < n; i++)131 {132 while(j && T[i] != P[j]) j = f[j];133 if(T[i] == P[j]) j++;134 if(j == m)135 {136 match[cnt++] = i-m+1;137 }138 }139 }140 141 char bit[10];142 void input()143 {144 for(int i = 0; i < n; i++)145 {146 scanf("%s", bit);147 da[i] = bit[7]-'0';148 for(int ii = 0; ii < 7; ii++)149 sa[i] = sa[i]*2+bit[ii]-'0';150 }151 for(int i = 0; i < m; i++)152 {153 scanf("%s", bit);154 db[m-1-i] = bit[7]-'0';155 for(int ii = 0; ii < 7; ii++)156 sb[i] = sb[i]*2+bit[ii]-'0';157 }158 sa[n] = sb[m] = -1;159 }160 161 void solve()162 {163 srand(time(0));164 KMP(sa, sb, n, m);165 if(cnt == 0)166 {167 puts("No");168 return;169 }170 puts("Yes");171 go();172 da.flip();173 db.flip();174 go();175 int ans = m;176 int pos = 0;177 for(int i = 0; i < cnt; i++)178 {179 if(ans == 0)180 break;181 int tmp = match[i];182 if(m-x[tmp] < ans)183 ans = m-x[tmp], pos = tmp;184 }185 186 printf("%d %d\n", ans, pos+1);187 }188 int main()189 {190 191 scanf("%d%d", &n,&m);192 input();193 solve();194 return 0;195 }
其实还遇到一道题，里面也用到了FFT,HDU 4954 Permanent，不过先挖个坑。
总结：FFT看上去很繁琐，了解了发现其实也就是一个很好利用的工具。很多东西要尽量主动去学习，如果之前了解过FFT，说不定遇到关键问题就不会像上面卡壳了。对于遇到的新的东西，也要沉下心来，不要只想着切自己学过的，感兴趣的知识，那绝不是进步的方法。
0 0
- FFT初步学习小结
- Struts2初步学习小结
- FFT小结
- FFT小结
- Flash AS学习初步小结
- IDC初步学习小小结
- FFT、NTT小结
- fft 学习
- FFT学习
- 数学专题小结：FFT算法
- [阶段小结]数据透视表的初步学习
- FCN与图像语义分割小结（学习初步指引）
- 指针初步说明小结（C++学习笔记）
- FFT学习笔记 FFT详解 一小时学会FFT
- FFT,NTT学习笔记
- FFT算法学习
- FFT学习小记
- FFT频谱学习
- HDOJ多校联合第五场
- bitset学习小记
- NEERC 2010, Eastern subregional contest
- 通过抓网络包发现服务器有返回值，但是sokcet中的BufferReader获取不到的问题
- MemSQL Start[c]UP 2.0 - Round 2
- FFT初步学习小结
- html基础
- HDOJ 题目3908Triple（数学）
- 新建xib适配iphone4尺寸的注意
- 实现一个键值对存储：目录
- TP调试
- 九度OJ--题目1091：棋盘游戏
- ffplay.c解码分析
- jsplumb 学习(1)
