FFT初步学习小结
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 in 115 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，说不定遇到关键问题就不会像上面卡壳了。对于遇到的新的东西，也要沉下心来，不要只想着切自己学过的，感兴趣的知识，那绝不是进步的方法。
