首页 > 代码库 > HDU 4609 3-idiots ——（FFT）
HDU 4609 3-idiots ——（FFT）
这是我接触的第一个关于FFT的题目，留个模板。
这题的题解见：http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html。
FFT的模板如下：
1 #include<bits/stdc++.h> 2 using namespace std; 3 const double pi = atan(1.0)*4; 4 struct Complex { 5 double x,y; 6 Complex(double _x=0,double _y=0) 7 :x(_x),y(_y) {} 8 Complex operator + (Complex &tt) { return Complex(x+tt.x,y+tt.y); } 9 Complex operator - (Complex &tt) { return Complex(x-tt.x,y-tt.y); } 10 Complex operator * (Complex &tt) { return Complex(x*tt.x-y*tt.y,x*tt.y+y*tt.x); } 11 }; 12 Complex a[15],b[15]; 13 void fft(Complex *a, int n, int rev) { 14 // n是(大于等于相乘的两个数组长度)2的幂次 ； 比如长度是5 ,那么 n = 8 2^2 < 5 2^3 > 5 15 // 从0开始表示长度，对a进行操作 16 // rev==1进行DFT，==-1进行IDFT 17 for (int i = 1,j = 0; i < n; ++ i) { 18 for (int k = n>>1; k > (j^=k); k >>= 1); 19 if (i<j) std::swap(a[i],a[j]); 20 } 21 for (int m = 2; m <= n; m <<= 1) { 22 Complex wm(cos(2*pi*rev/m),sin(2*pi*rev/m)); 23 for (int i = 0; i < n; i += m) { 24 Complex w(1.0,0.0); 25 for (int j = i; j < i+m/2; ++ j) { 26 Complex t = w*a[j+m/2]; 27 a[j+m/2] = a[j] - t; 28 a[j] = a[j] + t; 29 w = w * wm; 30 } 31 } 32 } 33 if (rev==-1) { 34 for (int i = 0; i < n; ++ i) a[i].x /= n,a[i].y /= n; 35 } 36 } 37 int main(){ 38 a[0] = Complex(0,0); // a[0]: x的0次项。 39 a[1] = Complex(1,0); 40 a[2] = Complex(2,0); 41 a[3] = Complex(3,0); 42 43 b[0] = Complex(3,0); 44 b[1] = Complex(2,0); 45 b[2] = Complex(1,0); 46 b[3] = Complex(0,0); 47 fft(a,8,1); 48 fft(b,8,1); 49 for(int i = 0 ; i < 8 ; i ++){ 50 a[i] = a[i] * b[i]; 51 } 52 fft(a,8,-1); 53 for(int i = 0 ; i < 8 ; i ++){ 54 cout << i << " " << a[i].x << endl;; 55 } 56 /* 57 * fft:nlogn求两个多项式相乘,(原来要n^2) 58 * 59 * f1 = 0x^0 + 1x^1 + 2x^2 + 3x^3 60 * f2 = 3 + 2x^1 + x^2 + 0x^3; 61 * 62 * f1 = a + b + c + d; 63 * f2 = x + y + z + k; 64 * f3 = __ 65 * dp[1] dp[2] dp[3] dp[4] dp[5]; 66 * c[0] c[1] c[2] c[3]; 67 */ 68 /* 69 0 0 * x^0 70 1 3 * x^1 71 2 8 * x^2 72 3 14 73 4 8 -----> 0x^0 + 3x^1 + 8x^2 ...... 3x^5 74 5 3 75 6 0 * x^6 76 7 0 * x^7 77 * */ 78 return 0; 79 }
关于这个模板有几点需要注意的：
1.在系数转化成整数时，会有精度误差，需要加eps。
2.假设a和b之前的长度都是n，卷积以后的大小应该是2*n-1，再考虑到fft中第二个参数n必须是大于等于卷积长度的2的幂次，因此最后的数组长度必须是n的4倍，也就是说这里的a数组大小应该开4倍才行。
最后，本题的AC代码如下：
1 #include<bits/stdc++.h> 2 using namespace std; 3 const double pi = atan(1.0)*4; 4 const int N = 1e5 + 5; 5 typedef long long ll; 6 7 struct Complex { 8 double x,y; 9 Complex(double _x=0,double _y=0) 10 :x(_x),y(_y) {} 11 Complex operator + (Complex &tt) { return Complex(x+tt.x,y+tt.y); } 12 Complex operator - (Complex &tt) { return Complex(x-tt.x,y-tt.y); } 13 Complex operator * (Complex &tt) { return Complex(x*tt.x-y*tt.y,x*tt.y+y*tt.x); } 14 }; 15 Complex a[N*4],b[N]; 16 void fft(Complex *a, int n, int rev) { 17 // n是(大于等于相乘的两个数组长度)2的幂次 ； 比如长度是5 ,那么 n = 8 2^2 < 5 2^3 > 5 18 // 从0开始表示长度，对a进行操作 19 // rev==1进行DFT，==-1进行IDFT 20 for (int i = 1,j = 0; i < n; ++ i) { 21 for (int k = n>>1; k > (j^=k); k >>= 1); 22 if (i<j) std::swap(a[i],a[j]); 23 } 24 for (int m = 2; m <= n; m <<= 1) { 25 Complex wm(cos(2*pi*rev/m),sin(2*pi*rev/m)); 26 for (int i = 0; i < n; i += m) { 27 Complex w(1.0,0.0); 28 for (int j = i; j < i+m/2; ++ j) { 29 Complex t = w*a[j+m/2]; 30 a[j+m/2] = a[j] - t; 31 a[j] = a[j] + t; 32 w = w * wm; 33 } 34 } 35 } 36 if (rev==-1) { 37 for (int i = 0; i < n; ++ i) a[i].x /= n,a[i].y /= n; 38 } 39 } 40 41 int A[N]; 42 ll num[N*2], sum[N*2]; 43 44 int main(){ 45 /*a[0] = Complex(0,0); // a[0]: x的0次项。 46 a[1] = Complex(1,0); 47 a[2] = Complex(2,0); 48 a[3] = Complex(3,0); 49 50 b[0] = Complex(3,0); 51 b[1] = Complex(2,0); 52 b[2] = Complex(1,0); 53 b[3] = Complex(0,0); 54 fft(a,8,1); 55 fft(b,8,1); 56 for(int i = 0 ; i < 8 ; i ++){ 57 a[i] = a[i] * b[i]; 58 } 59 fft(a,8,-1); 60 for(int i = 0 ; i < 8 ; i ++){ 61 cout << i << " " << a[i].x << endl;; 62 }*/ 63 //a[0] = Complex(0, 0); b[0] = Complex(0, 0); 64 int T; scanf("%d",&T); 65 while(T--) 66 { 67 int n; scanf("%d",&n); 68 memset(num,0,sizeof num); 69 for(int i=1;i<=n;i++) 70 { 71 scanf("%d",A+i); 72 num[A[i]]++; 73 } 74 sort(A+1, A+1+n); 75 int len = A[n]; 76 int LIM = 1; 77 while(1) 78 { 79 if(LIM >= len*2+1) break; 80 else LIM <<= 1; 81 } 82 for(int i=0;i<=len;i++) 83 { 84 a[i] = Complex(num[i], 0); 85 } 86 for(int i=len+1;i<LIM;i++) 87 { 88 a[i] = Complex(0, 0); 89 } 90 fft(a,LIM,1); 91 for(int i=0;i<LIM;i++) a[i] = a[i] * a[i]; 92 fft(a,LIM,-1); 93 // finish fft 94 len = len * 2; 95 for(int i=0;i<=len;i++) num[i] = (ll)(a[i].x + 0.5); 96 for(int i=1;i<=n;i++) num[A[i]*2]--; // 减去两个相同的组合 97 for(int i=1;i<=len;i++) num[i] /= 2; // 选择的是无序的 98 for(int i=1;i<=len;i++) sum[i] = sum[i-1] + num[i]; 99 ll ans = 0; 100 for(int i=1;i<=n;i++) 101 { 102 int now = A[i]; 103 ans += sum[len] - sum[now]; // 对于每个数，加起来比它大的都是可行的 104 ans -= 1LL * (i-1) * (n-i); // 减去一个大的一个小的组合的情况 105 ans -= 1LL * (n-1); // 减去自己和任意一个组合的情况 106 ans -= 1LL * (n-i) * (n-i-1) / 2; // 减去比它大的两个组合的情况 107 // 剩下的就是两个小的加起来比A[i]大的情况 108 } 109 ll all = 1LL * n*(n-1)*(n-2) / 6; 110 printf("%.7f\n",1.0*ans/all); 111 } 112 return 0; 113 }
HDU 4609 3-idiots ——（FFT）
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。