这是我接触的第一个关于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 }
FFT模板
关于这个模板有几点需要注意的：
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代码
时间： 2024-07-29 07:06:43