HDU4609 3-idiots（母函数 + FFT）
题目
Source
http://acm.hdu.edu.cn/showproblem.php?pid=4609
Description
King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king's forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn't pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.
Input
An integer T(T≤100) will exist in the first line of input, indicating the number of test cases.
Each test case begins with the number of branches N(3≤N≤105).
The following line contains N integers a_i (1≤a_i≤105), which denotes the length of each branch, respectively.
Output
Output the probability that their branches can form a triangle, in accuracy of 7 decimal places.
Sample Input
2
4
1 3 3 4
4
2 3 3 4
Sample Output
0.5000000
1.0000000
分析
题目大概说，有n条边，长度在[1,100000]，从中任选三条边，问组成三角形的概率是多少。
这题思路好绕。
详见：http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html
大概这样：
- 首先求出两两组合成长度和位i的方案数tot[i]
- 然后从小到大枚举三角形最长的边len，那么现在要确定另外两条边合法组合的方案数
这个方案数为：
Σtot[i]（i>len）- 最长边与其他边组合的方案数 - 长度小于最长边与长度大于等于最长边的其他边组合的方案数 - 长度大于等于最长边的其他边之间组合的方案数
Σtot[i]可以通过前缀和差分求得，其他部分简单的组合计算。
而tot[i]如何求？
构造一个多项式，$\sum_ic_ix^i$，表示长度$i$的边有$c_i$个，而$\sum_ic_ix^i \times \sum_ic_ix^i$的结果$\sum_iC_ix^i$就表示有先后顺序且放回地选出两条边长度和为$i$的方案数为$C_i$
这其实是母函数吧。。然后利用FFT去快速求得这个多项式的值。
不过这个多项式求出来还不是最终的，再减去自身与自身组合，最后除以2即是要的tot[i]。
代码
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define MAXN 277777
const double PI=acos(-1.0);
struct Complex{
double real,imag;
Complex(double _real,double _imag):real(_real),imag(_imag){}
Complex(){}
Complex operator+(const Complex &cp) const{
return Complex(real+cp.real,imag+cp.imag);
}
Complex operator-(const Complex &cp) const{
return Complex(real-cp.real,imag-cp.imag);
}
Complex operator*(const Complex &cp) const{
return Complex(real*cp.real-imag*cp.imag,real*cp.imag+cp.real*imag);
}
void setValue(double _real=0,double _imag=0){
real=_real; imag=_imag;
}
};
int len;
Complex wn[MAXN],wn_anti[MAXN];
void FFT(Complex y[],int op){
for(int i=1,j=len>>1,k; i<len-1; ++i){
if(i<j) swap(y[i],y[j]);
k=len>>1;
while(j>=k){
j-=k;
k>>=1;
}
if(j<k) j+=k;
}
for(int h=2; h<=len; h<<=1){
Complex Wn=(op==1?wn[h]:wn_anti[h]);
for(int i=0; i<len; i+=h){
Complex W(1,0);
for(int j=i; j<i+(h>>1); ++j){
Complex u=y[j],t=W*y[j+(h>>1)];
y[j]=u+t;
y[j+(h>>1)]=u-t;
W=W*Wn;
}
}
}
if(op==-1){
for(int i=0; i<len; ++i) y[i].real/=len;
}
}
void Convolution(Complex A[],Complex B[],int n){
for(len=1; len<(n<<1); len<<=1);
for(int i=n; i<len; ++i){
A[i].setValue();
B[i].setValue();
}
for(int i=0; i<=len; ++i){
wn[i].setValue(cos(2.0*PI/i),sin(2.0*PI/i));
wn_anti[i].setValue(wn[i].real,-wn[i].imag);
}
FFT(A,1); FFT(B,1);
for(int i=0; i<len; ++i){
A[i]=A[i]*B[i];
}
FFT(A,-1);
}
Complex A[MAXN],B[MAXN];
int a[111111],cnt[111111];
long long tot[MAXN];
int main(){
int t,n;
scanf("%d",&t);
while(t--){
scanf("%d",&n);
memset(cnt,0,sizeof(cnt));
int maxa=0;
for(int i=0; i<n; ++i){
scanf("%d",a+i);
++cnt[a[i]];
maxa=max(maxa,a[i]);
}
for(int i=0; i<=maxa; ++i){
A[i].setValue(cnt[i]);
B[i].setValue(cnt[i]);
}
Convolution(A,B,maxa+1);
for(int i=0; i<len; ++i){
tot[i]=(long long)(A[i].real+0.5);
}
for(int i=0; i<n; ++i){
--tot[a[i]+a[i]];
}
for(int i=0; i<len; ++i){
tot[i]>>=1;
}
for(int i=1; i<len; ++i){
tot[i]+=tot[i-1];
}
sort(a,a+n);
long long res=0;
for(int i=0; i<n; ++i){
long long tmp=tot[len-1]-tot[a[i]];
tmp-=n-1; // 本身与其它的组合
tmp-=(long long)i*(n-i-1); // 小于与大于等于的组合
if(n-i-1>1) tmp-=(long long)(n-i-1)*(n-i-2)/2; // 大于等于之间的组合
res+=tmp;
}
printf("%.7f\n",res*1.0/((long long)n*(n-1)*(n-2)/6));
}
return 0;
}
