3-idiots
题意
t组测试，每组测试给出n条木棒 n<1e5
然后输入n条木棒的长度ai,问等概率取出三条木棒能构成三角形的概率
题解
参考题解A 参考题解B
首先将读入信息从每个木棒长度 变为 每种长度几根木棒
例如lenth={1,3,3,4} 变为 num={0,1,0,2,1}
即长度为0的有0根,1的有1根,…..
那么如果从它们中任选两根组成三角形的两边和 那么可以得到
num={0, 0, 1, 0, 4, 2, 4, 4, 1}
(其中有AB互换得到的和以及同一根木棍自己和自己相加得到的和,之后会减去)
假如用生成函数来描述 那么 {0,1,0,2,1} 记为$f(x)=0x^0+1x^1+0x^2+2x^3+1x^4$
而两两匹配得到的{0, 0, 1, 0, 4, 2, 4, 4, 1}
实际上就是$f(x)×f(x)=0x^0+0x^1+1x^2+0x^3+4x^4+2x^5+4x^6+4x^7+1x^8$
其中次幂i表示长度,系数ai表示有多少根木棍=该长度i
于是我们可以利用FFT求出木棍两两组合得到新的长度各有多少种
另外需要容斥排除掉一些情况，在代码的注释中会有解释
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
| #include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define endl '\n'
#define Turnoff std::ios::sync_with_stdio(false)
const ll Max=3e5+5;
const double Pi=acos(-1);
/*
*/
struct Complex{
double r,i;
Complex(){r=0,i=0;}
Complex(double real,double imag):r(real),i(imag){};
};
Complex operator + (Complex a, Complex b){return Complex(a.r+b.r,a.i+b.i);}
Complex operator - (Complex a, Complex b){return Complex(a.r-b.r,a.i-b.i);}
Complex operator * (Complex a, Complex b){return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
ll rev[Max],len=0,lim=1,ans[Max];
Complex a[Max],b[Max];///原多项式系数
ll sum[Max],lenth[Max],num[Max];
///预处理位逆变换
inline void getrev(int lena,int lenb){
int limit=1,lenth=0;
while(limit<=lena+lenb)limit<<=1,lenth++;
for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lenth-1));
lim=limit,len=lenth;
return;
}
inline void FFT(Complex *A,int opt){
for(int i=0;i<lim;i++)if(i<rev[i]) swap(A[i],A[rev[i]]);///求出要迭代的序列
for(int mid=1;mid<lim;mid<<=1){///待合并区间长度的一半
Complex Wn=Complex(cos(Pi/mid),opt*sin(Pi/mid));
for(int R=mid<<1,i=0;i<lim;i+=R){///R是区间的右端点，i表示前已经到哪个位置了
Complex W=Complex(1,0);
for(int k=0;k<mid;k++){///枚举左半部分
Complex x=A[i+k],y=W*A[i+mid+k];///蝴蝶操作
A[i+k]=x+y;
A[i+mid+k]=x-y;
W=W*Wn;
}
}
}
if(opt==-1)for(int i=0;i<lim;i++)A[i].r=(A[i].r/lim+0.5);
}
int main() {
//Turnoff;
int t;
scanf("%d",&t);
while(t--){
int n;scanf("%d",&n);
ll maxlen=0,tot=1ll*n*(n-1)*(n-2)/6;
memset(a,0,sizeof a);
memset(sum,0,sizeof sum);
//memset(num,0,sizeof num);
for(int i=0;i<n;i++){
scanf("%lld",lenth+i);
a[lenth[i]].r++;
//num[lenth[i]]++;
maxlen=max(maxlen,lenth[i]);
}
getrev(maxlen,maxlen);FFT(a,1);
maxlen*=2;
for(int i=0;i<lim;i++)a[i]=a[i]*a[i];
FFT(a,-1);
for(int i=0;i<lim;i++)sum[i]=a[i].r;
for(int i=0;i<n;i++)sum[lenth[i]*2]--;
///同一根木棒自己加自己 不合法减去
for(int i=0;i<=maxlen;i++)sum[i]/=2;
/// 两根木棒AB互换后BA重复统计 不合法/2
for(int i=1;i<=maxlen;i++)sum[i]+=sum[i-1];
/// 做前缀和处理
ll ans=0;
sort(lenth,lenth+n);
///排序构成三角形的第三边(最长边)长度
///枚举最长边
for(int i=0;i<n;i++){
ans+=sum[maxlen]-sum[lenth[i]];
///根据定理两边之和大于第三边得到两边和大于枚举边长度lenth[i]的总数
///我们假定lenth[i]最长
ans-=1ll*i*(n-i-1);
///两边中有一个大于lenth[i]一个小于lenth[i]
ans-=1ll*(n-i-1)*(n-i-2)/2;
///两边都大于lenth[i]
ans-=1ll*(n-1);
///另外两边中有一个就是lenth[i]本身
}
//cout<<ans<<" "<<tot<<endl;
printf("%.7f\n",1.0*ans/tot);
}
}
|