HDU 4609 3-idiots ——（FFT）
阿新 • • 發佈：2017-08-06
ace 代碼 ++ per str [1] 兩個 isp 題目 
這是我接觸的第一個關於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 operatorFFT模板+ (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 }AC代碼
HDU 4609 3-idiots ——（FFT）