bzoj 3513 [MUTC2013]idiots FFT 生成函数

[MUTC2013]idiots

Time Limit: 20 Sec  Memory Limit: 128 MB
Submit: 806  Solved: 265
[Submit][Status][Discuss]

Description

给定n个长度分别为a_i的木棒，问随机选择3个木棒能够拼成三角形的概率。

Input

第一行T(T<=100)，表示数据组数。

接下来若干行描述T组数据，每组数据第一行是n，接下来一行有n个数表示a_i。

3≤N≤10^5,1≤a_i≤10^5

Output

T行，每行一个整数，四舍五入保留7位小数。

Sample Input

2
4
1 3 3 4
4
2 3 3 4

Sample Output

0.5000000
1.0000000

HINT

T<=20

N<=100000

Source

By sbullet

题解：先想想暴力怎么做，暴力的话可以枚举三个点对吧，没什么问题，

　　　这样是n^3的复杂度，怎么优化很关键，可以用什么数据结构还是前缀和，之类n^2logn n^2都是可以优化到的。

　　　下面我们可以枚举最大值，如果x+y<=z那么这三个数时组成不了三角形的，那么有多少呢？

　　　就是生成函数FFT优化，然后前缀和一下，就可以知道多少对二元组比z小，去重的很简单的。

　　　然后用总方案减不合法方案，就可以了。

#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iostream>

#define pi acos(-1)
#define ll long long
#define N 1000007
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    return x*f;
}

int n,mx,num,L;
int val[N],rev[N];
ll sum[N],ans,zhi;
struct comp
{
    double r,v;
    comp(){r=v=0.0;}
    comp(double x,double y){r=x,v=y;}
    void init(){r=v=0.0;}
    friend inline comp operator+(comp x,comp y){return comp(x.r+y.r,x.v+y.v);}
    friend inline comp operator-(comp x,comp y){return comp(x.r-y.r,x.v-y.v);}
    friend inline comp operator*(comp x,comp y){return comp(x.r*y.r-x.v*y.v,x.r*y.v+x.v*y.r);}
    friend inline comp operator/(comp x,int y){return comp(x.r/y,x.v/y);}
}a[N];

void FFT(comp *a,int flag)
{
    for (int i=0;i<num;i++)
        if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1;i<num;i<<=1)
    {
        comp wn=comp(cos(pi/i),flag*sin(pi/i));
        for (int j=0;j<num;j+=(i<<1))
        {
            comp w=comp(1,0);
            for (int k=0;k<i;k++,w=w*wn)
            {
                comp x=a[j+k],y=a[j+k+i]*w;
                a[j+k]=x+y,a[j+k+i]=x-y;
            }
        }
    }
    if (flag==-1) for (int i=0;i<num;i++) a[i]=a[i]/num;
}
int main()
{
    int T=read();
    while(T--)
    {
        n=read(),mx=0;
        for (int i=1;i<=n;i++)
            val[i]=read(),mx=max(mx,val[i]<<1);
        for (num=1,L=0;num<=mx;num<<=1,L++);if (L) L--;
        for (int i=0;i<num;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
        memset(sum,0,sizeof(sum));
        for (int i=0;i<num;i++) a[i].init(); 
        for (int i=1;i<=n;i++)
            a[val[i]].r+=1,sum[val[i]<<1]--;
        FFT(a,1);
        for(int i=0;i<num;i++)
            a[i]=a[i]*a[i];
        FFT(a,-1);
        for (int i=0;i<num;i++)
            sum[i]+=(ll)(a[i].r+0.5);
        for (int i=1;i<=mx;i++)
            sum[i]+=sum[i-1];
        ans=(ll)n*(n-1)*(n-2)/6,zhi=0;
        for (int i=1;i<=n;i++)
            zhi+=sum[val[i]];
        printf("%.7lf
",(1-0.5*zhi/ans));//乘0.5是因为选x,y和y,x一样的 
    }
}