HDU-4609 3-idiots FFT

　　题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4609　

　　题意：给n个棒子，求任意组成3根能形成三角形的概率。

　　用FFT求出任意两根棒子组合成新的长度，这种长度的组合有多少种。

　　用cnt[i]表示和的长度为i的棒子组合有多少种。

　　(1) 首先去掉cnt里重复的部分，一根棒子不能与自己组合，所以有cnt[a[i] * 2] --

　　(2) 任意两个棒子组合的顺序我们不需要考虑，而实际上他们算了两次，所以有cnt[i] /=2

　　(3) 枚举所有棒子，假设这根棒子a[i]是组成三角形其中最长的那根，我们首先求出另外两根的长度和比它的长的总种数，即sigma(cnt[j]) (a[i]<j<=crest*2，crest是最长的棒子长度)

　　(4) 另外两根中，这根被枚举的棒子不能再出现了，所以要减去n-1

　　(5) 另外两根中，可能有一根长度大于a[i]，另一根小于a[i]，所以要减去(n-i) * (i-1)

　　(6) 另外两根中，可能两根长度都大于a[i]，再减去(n-i) * (n-i-1) / 2

  1 //STATUS:C++_AC_2140MS_7128KB
  2 #include <functional>
  3 #include <algorithm>
  4 #include <iostream>
  5 //#include <ext/rope>
  6 #include <fstream>
  7 #include <sstream>
  8 #include <iomanip>
  9 #include <numeric>
 10 #include <cstring>
 11 #include <cassert>
 12 #include <cstdio>
 13 #include <string>
 14 #include <vector>
 15 #include <bitset>
 16 #include <queue>
 17 #include <stack>
 18 #include <cmath>
 19 #include <ctime>
 20 #include <list>
 21 #include <set>
 22 #include <map>
 23 using namespace std;
 24 //#pragma comment(linker,"/STACK:102400000,102400000")
 25 //using namespace __gnu_cxx;
 26 //define
 27 #define pii pair<int,int>
 28 #define mem(a,b) memset(a,b,sizeof(a))
 29 #define lson l,mid,rt<<1
 30 #define rson mid+1,r,rt<<1|1
 31 #define PI acos(-1.0)
 32 //typedef
 33 typedef __int64 LL;
 34 typedef unsigned __int64 ULL;
 35 //const
 36 const int N=400010;
 37 const int INF=0x3f3f3f3f;
 38 const int MOD=10007,STA=8000010;
 39 const LL LNF=1LL<<55;
 40 const double EPS=1e-4;
 41 const double OO=1e30;
 42 const int dx[4]={-1,0,1,0};
 43 const int dy[4]={0,1,0,-1};
 44 const int day[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};
 45 //Daily Use ...
 46 inline int sign(double x){return (x>EPS)-(x<-EPS);}
 47 template<class T> T gcd(T a,T b){return b?gcd(b,a%b):a;}
 48 template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
 49 template<class T> inline T lcm(T a,T b,T d){return a/d*b;}
 50 template<class T> inline T Min(T a,T b){return a<b?a:b;}
 51 template<class T> inline T Max(T a,T b){return a>b?a:b;}
 52 template<class T> inline T Min(T a,T b,T c){return min(min(a, b),c);}
 53 template<class T> inline T Max(T a,T b,T c){return max(max(a, b),c);}
 54 template<class T> inline T Min(T a,T b,T c,T d){return min(min(a, b),min(c,d));}
 55 template<class T> inline T Max(T a,T b,T c,T d){return max(max(a, b),max(c,d));}
 56 //End
 57 //复数结构体
 58 struct complex
 59 {
 60     double r,i;
 61     complex(double _r = 0.0,double _i = 0.0)
 62     {
 63         r = _r; i = _i;
 64     }
 65     complex operator +(const complex &b)
 66     {
 67         return complex(r+b.r,i+b.i);
 68     }
 69     complex operator -(const complex &b)
 70     {
 71         return complex(r-b.r,i-b.i);
 72     }
 73     complex operator *(const complex &b)
 74     {
 75         return complex(r*b.r-i*b.i,r*b.i+i*b.r);
 76     }
 77 };
 78 /*
 79  * 进行FFT和IFFT前的反转变换。
 80  * 位置i和 （i二进制反转后位置）互换
 81  * len必须去2的幂
 82  */
 83 void change(complex y[],int len)
 84 {
 85     int i,j,k;
 86     for(i = 1, j = len/2;i < len-1; i++)
 87     {
 88         if(i < j)swap(y[i],y[j]);
 89         //交换互为小标反转的元素，i<j保证交换一次
 90         //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的
 91         k = len/2;
 92         while( j >= k)
 93         {
 94             j -= k;
 95             k /= 2;
 96         }
 97         if(j < k) j += k;
 98     }
 99 }
100 /*
101  * 做FFT
102  * len必须为2^k形式，
103  * on==1时是DFT，on==-1时是IDFT
104  */
105 void FFT(complex y[],int len,int on)
106 {
107     change(y,len);
108     for(int h = 2; h <= len; h <<= 1)
109     {
110         complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
111         for(int j = 0;j < len;j+=h)
112         {
113             complex w(1,0);
114             for(int k = j;k < j+h/2;k++)
115             {
116                 complex u = y[k];
117                 complex t = w*y[k+h/2];
118                 y[k] = u+t;
119                 y[k+h/2] = u-t;
120                 w = w*wn;
121             }
122         }
123     }
124     if(on == -1)
125         for(int i = 0;i < len;i++)
126             y[i].r /= len;
127 }
128 
129 LL sum[N/2],tot[N/2];
130 int cnt[N/4],l[N/4];
131 complex a[N],b[N];
132 int T,n;
133 
134 int main(){
135  //   freopen("in.txt","r",stdin);
136     int i,j,len,hig;
137     LL ans;
138     scanf("%d",&T);
139     while(T--)
140     {
141         scanf("%d",&n);
142         mem(cnt,0);hig=0;
143         for(i=1;i<=n;i++){
144             scanf("%d",&l[i]);
145             cnt[l[i]]++;
146             hig=Max(hig,l[i]);
147         }
148         sort(l+1,l+n+1);
149         hig++;
150         for(len=1;len<(hig<<1);len<<=1);
151         for(i=0;i<hig;i++)
152             a[i]=b[i]=complex(cnt[i],0);
153         for(;i<=len;i++)
154             a[i]=b[i]=complex(0,0);
155 
156         FFT(a,len,1);
157         FFT(b,len,1);
158         for(i=0;i<=len;i++)
159             a[i]=a[i]*b[i];
160         FFT(a,len,-1);
161         len=(hig<<1)-1;
162         for(i=0;i<=len;i++)
163             tot[i]=(LL)(a[i].r+0.5);
164         for(i=1;i<=n;i++)
165             tot[l[i]<<1]--;
166         for(i=1;i<=len;i++)tot[i]>>=1;
167         sum[0]=0;
168         for(i=1;i<=len;i++)
169             sum[i]=tot[i]+sum[i-1];
170         ans=0;
171         for(i=1;i<=n;i++){
172             ans+=sum[len]-sum[l[i]];
173             ans-=(LL)(n-1)+(LL)(i-1)*(n-i)+(LL)(n-i)*(n-i-1)/2;
174         }
175 
176         printf("%.7lf
",(double)ans/n/(n-1)/(n-2)*6);
177     }
178     return 0;
179 }