「Hdu4609」3-idiots

来源和评测点

2013 Multi-University Training Contest 1
Hdu4609
Caioj1454

题目

「题目大意」
有T组数据
每组数据给出n条线段，问任意取三条，可以组成三角形的概率
「输入格式」
开头一行输入T（T<=100）
下来T组数据，每组数据第一行输入一个n（3<=n<=10^5）
第二行输入n个数，表示n条线段
线段长度（1<=n<=10^5）
「输出格式」
每组数据输出一个数p
表示可以组成三角形的概率
保留七位小数
「输入样例」
2
4
1 3 3 4
4
2 3 3 4
「输出样例」
0.5000000
1.0000000

分析

其他快速傅里叶变换题目参见：Tag-快速傅里叶变换

大致思路：
计算出两条边的和，枚举第三边表示能和它搭配的数量，经过去重后，除以总组合数就是概率

两边之和：FFT，参考另外几道简单点的题
去重1：两边之和不能是来自同一条边
去重2：两边有一条比第三边大，另一个小
去重3：两边都比第三边大
去重4：两边有任何一条=第三边
（去重2、3、4后剩下的就是两边都比第三边小的情况）

输入的是q[i]，$A(q[i])=是q[i]的数量$
NUM(i)表示两边之和=i的情况个数，用去重1
$$
NUM(i)=\frac{ A^2(i)-B(i) }{2}
$$

去重2、3、4利用排序，并用前缀和SUM获得区间和
$ ANS(i)=SUM(rn-1)-SUM(q[i]) $
$ 去重2：ANS(i)-=i\times (n-i-1) $
$ 去重3：ANS(i)-=\frac{ (n-i-2)\times (n-i-1) }{2} $
$ 去重4：ANS(i)-=n-1 $
$$
tot=\frac{ n\times (n-1)\times (n-2) }{2}
$$

答案就是
$$
\frac{ \sum_{i=0}^{n-1} ANS(i) }{tot}
$$

代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
struct complex
{
    double r,i;
    complex() { r=i=0.0; }
    complex(double rr,double ii) { r=rr;i=ii; }
    complex operator + (complex b) { return complex(r+b.r,i+b.i); }
    complex operator - (complex b) { return complex(r-b.r,i-b.i); }
    complex operator * (complex b) { return complex(r*b.r-i*b.i,r*b.i+i*b.r); }
};
const int MAXN=100000*4;
const double PI=acos(-1.0);

int R[MAXN+10];
void FFT(complex *a,int n,int op)
{
    for(int i=0;i<=n-1;i++) if(i<R[i]) swap(a[i],a[R[i]]);
    for(int l=1;l<=n/2;l*=2)
    {
        complex w1( cos(op*PI/l),sin(op*PI/l) );
        for(int j=0;j<=n-1;j+=l*2)
        {
            complex wk(1.0,0.0);
            for(int k=0;k<=l-1;k++,wk=wk*w1)
            {
                complex x=a[j+k],y=wk*a[j+l+k];
                a[j+k]=x+y;
                a[j+l+k]=x-y;
            }
        }
    }
    if(op==-1) for(int i=0;i<=n-1;i++) a[i].r=a[i].r/double(n);
}
//http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html
//http://acm.HDU.edu.cn/showproblem.php?pid=4609
//http://blog.csdn.net/rose_max/article/details/77370572
complex A[MAXN+10],A2[MAXN+10],B[MAXN+10];

int q[MAXN+10];
double NUM[MAXN+10];
double SUM[MAXN+10];
int main()
{
    int T;scanf("%d",&T);
    while(T--)
    {
        memset(A,0,sizeof(A));
        memset(B,0,sizeof(B));
        memset(SUM,0,sizeof(SUM));
        
        int n;scanf("%d",&n);int mx=0;
        for(int i=0;i<=n-1;i++)
        {
            scanf("%d",&q[i]);
            A[q[i]].r++;
            B[q[i]*2].r++;
            if(q[i]>mx) mx=q[i];
        }
        
        sort(q,q+n);
        
        int rn=1,rm=(mx)+(mx)+1,l=0;
        while(rn<rm) rn*=2,l++;
        for(int i=0;i<=rn-1;i++) R[i]=(R[i>>1]>>1) | ( (i&1)<<(l-1) );
        
        FFT(A,rn,1);
        for(int i=0;i<=rn-1;i++) A2[i]=A[i]*A[i];
        FFT(A2,rn,-1);
        
        for(int i=0;i<=rn-1;i++) NUM[i]=(A2[i].r-B[i].r)/2.0;
        for(int i=0;i<=rn-1;i++) SUM[i]=SUM[i-1]+NUM[i];
        
        double ans=0;
        for(int i=0;i<=n-1;i++)//第三边
        {
            ans+=SUM[rn-1]-SUM[q[i]];
            ans-=double(i)*double(n-i-1);//一大一小 
            ans-=double(n-i-2)*double(n-i-1)/2;//两大
            ans-=n-1;//某边=第三边 
        }
        double tot=double(n)*double(n-1)*double(n-2)/6.0;
        printf("%.7lf\n",ans/tot);
    }
}