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小，去重的很简单的。

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

 1 #include<cstring>
 2 #include<cmath>
 3 #include<algorithm>
 4 #include<cstdio>
 5 #include<iostream>
 6 
 7 #define pi acos(-1)
 8 #define ll long long
 9 #define N 1000007
10 using namespace std;
11 inline int read()
12 {
13     int x=0,f=1;char ch=getchar();
14     while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
15     while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
16     return x*f;
17 }
18 
19 int n,mx,num,L;
20 int val[N],rev[N];
21 ll sum[N],ans,zhi;
22 struct comp
23 {
24     double r,v;
25     comp(){r=v=0.0;}
26     comp(double x,double y){r=x,v=y;}
27     void init(){r=v=0.0;}
28     friend inline comp operator+(comp x,comp y){return comp(x.r+y.r,x.v+y.v);}
29     friend inline comp operator-(comp x,comp y){return comp(x.r-y.r,x.v-y.v);}
30     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);}
31     friend inline comp operator/(comp x,int y){return comp(x.r/y,x.v/y);}
32 }a[N];
33 
34 void FFT(comp *a,int flag)
35 {
36     for (int i=0;i<num;i++)
37         if (i<rev[i]) swap(a[i],a[rev[i]]);
38     for (int i=1;i<num;i<<=1)
39     {
40         comp wn=comp(cos(pi/i),flag*sin(pi/i));
41         for (int j=0;j<num;j+=(i<<1))
42         {
43             comp w=comp(1,0);
44             for (int k=0;k<i;k++,w=w*wn)
45             {
46                 comp x=a[j+k],y=a[j+k+i]*w;
47                 a[j+k]=x+y,a[j+k+i]=x-y;
48             }
49         }
50     }
51     if (flag==-1) for (int i=0;i<num;i++) a[i]=a[i]/num;
52 }
53 int main()
54 {
55     int T=read();
56     while(T--)
57     {
58         n=read(),mx=0;
59         for (int i=1;i<=n;i++)
60             val[i]=read(),mx=max(mx,val[i]<<1);
61         for (num=1,L=0;num<=mx;num<<=1,L++);if (L) L--;
62         for (int i=0;i<num;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
63         memset(sum,0,sizeof(sum));
64         for (int i=0;i<num;i++) a[i].init(); 
65         for (int i=1;i<=n;i++)
66             a[val[i]].r+=1,sum[val[i]<<1]--;
67         FFT(a,1);
68         for(int i=0;i<num;i++)
69             a[i]=a[i]*a[i];
70         FFT(a,-1);
71         for (int i=0;i<num;i++)
72             sum[i]+=(ll)(a[i].r+0.5);
73         for (int i=1;i<=mx;i++)
74             sum[i]+=sum[i-1];
75         ans=(ll)n*(n-1)*(n-2)/6,zhi=0;
76         for (int i=1;i<=n;i++)
77             zhi+=sum[val[i]];
78         printf("%.7lf
",(1-0.5*zhi/ans));//乘0.5是因为选x,y和y,x一样的 
79     }
80 }