HDU - 4609 3-idiots （FFT入门）

3-idiots

Time Limit: 10000/5000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 7656    Accepted Submission(s): 2672

 

Problem 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根木棍中挑选一根，问这三个傻子挑选出来的木棍能组成三角形的概率是多少。

题目思路：题目的问题可以转化成，求在n条线段中任意选出3条能组成三角形的概率是多少。在看题解之前，对于我这种不会FFT的菜鸡来说，这题根本无从下手。在学习了别人的做法之后也算理解了如何做。

我们知道三条线段能组成三角形的条件为两条短的边长之和大于长的边，那么我们现在只要求出任意两边的边长之和，然后再枚举最长边，即可求出能组成的三角形的个数了。

我们可以用一个数组num[i]来表示边长为 i 的边的数量。我们知道对于多项式乘法，两项相乘之后得到的数的下标正好是这两项的下标之和。

如果我们令数组num代表一个多项式，使得nu数组和num数组相乘，这样我们就能求出任意两边的边长之和的数量了。但是直接对于num做多项式乘法的话复杂度是O(n^2)的，所以这个时候我们要借助FFT的帮助，在O(nlogn)的时间内求出num数组和num数组的卷积。

卷积的第 i 项就代表着两边长度和为 i 的数量。但是由于题目中的要求每根木棍只能挑选一次，所以要减掉选了两根相同的木棍的情况，由于选取木棍没有顺序关系，所以还得将方案数除2。接着就能借助一个前缀和sum求出两条边长大于 a[i] 的方案数了。

对a[i]就能O(n)的枚举a[i]作为最长边的所有情况了，对于a[i]作为最长边的情况，我们还要减掉选出来的一条边的边长大于a[i]，一条边等于a[i]的情况，以及两条边都大于a[i]情况。最后再除以所有方案数即可得到答案。

具体实现看代码：

#include <bits/stdc++.h>
#define fi first
#define se second
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
#define lowbit(x) x&-x
#define pb push_back
#define MP make_pair
#define clr(a) memset(a,0,sizeof(a))
#define _INF(a) memset(a,0x3f,sizeof(a))
#define FIN freopen("in.txt","r",stdin)
#define fuck(x) cout<<"["<<x<<"]"<<endl
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int>pii;
//head
const int MX=1e5+7;

int n,a[MX],_;
ll num[MX],sum[MX];
const double pi = acos(-1.0);
int len,mx;
ll res[MX<<2];
struct Complex {
    double r,i;
    Complex(double r=0,double i=0):r(r),i(i) {};
    Complex operator+(const Complex &rhs) {return Complex(r + rhs.r,i + rhs.i);}
    Complex operator-(const Complex &rhs) {return Complex(r - rhs.r,i - rhs.i);}
    Complex operator*(const Complex &rhs) {return Complex(r*rhs.r - i*rhs.i,i*rhs.r + r*rhs.i);}
} va[MX<<2],vb[MX<<2];
void rader(Complex F[],int len) { //len = 2^M,reverse F[i] with  F[j] j为i二进制反转
    int j = len >> 1;
    for(int i = 1; i < len - 1; ++i) {
        if(i < j) swap(F[i],F[j]);  // reverse
        int k = len>>1;
        while(j>=k) {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}
void FFT(Complex F[],int len,int t) {
    rader(F,len);
    for(int h=2; h<=len; h<<=1) {
        Complex wn(cos(-t*2*pi/h),sin(-t*2*pi/h));
        for(int j=0; j<len; j+=h) {
            Complex E(1,0); //旋转因子
            for(int k=j; k<j+h/2; ++k) {
                Complex u = F[k];
                Complex v = E*F[k+h/2];
                F[k] = u+v;
                F[k+h/2] = u-v;
                E=E*wn;
            }
        }
    }
    if(t==-1)   //IDFT
        for(int i=0; i<len; ++i)
            F[i].r/=len;
}
void Conv(Complex a[],Complex b[],int len) { //求卷积
    FFT(a,len,1);
    FFT(b,len,1);
    for(int i=0; i<len; ++i) a[i] = a[i]*b[i];
    FFT(a,len,-1);
}
void work() {
    Conv(va,vb,len);
    for(int i=0; i<len; ++i)res[i]=va[i].r + 0.5;
}

int main(){
	//FIN;
	for(scanf("%d",&_);_;_--){
		clr(num);
		scanf("%d",&n);
		int len1=0;
		for(int i=0;i<n;i++){
			scanf("%d",&a[i]);
			num[a[i]]++;
			len1=max(len1,a[i]);
		}
		sort(a,a+n);
		len1++;
		len=1;
		while(len<2*len1) len<<=1;
		for(int i=0;i<len;i++){
			if(i<len1) va[i]=vb[i]=Complex(num[i],0);
			else va[i]=vb[i]=Complex(0,0);
		}
		work();
		len=a[n-1]*2;
		for(int i=0;i<n;i++) res[a[i]+a[i]]--;
		for(int i=1;i<=len;i++) res[i]/=2;
		for(int i=1;i<=len;i++) sum[i]=sum[i-1]+res[i];
		ll ans=0,tot=(ll)n*(n-1)*(n-2)/6;
		for(int i=0;i<n;i++){
			//统计a[i]为三角形的最长边的所有情况
			ans+=(sum[len]-sum[a[i]]);//在剩下的n-1条边中选出两条边的长度和大于a[i]
			ans-=(ll)(n-i-1)*i;//减掉一条边大于a[i],一条边小于a[i]的情况
			ans-=(ll)(n-i-1)*(n-i-2)/2;//减掉两条边都大于a[i]的情况
			ans-=(n-1);//减掉一条边为a[i]的情况
		}
		printf("%.7f\n",1.0*ans/tot);
	}
	return 0;
}