隐藏
「bzoj3513」「MUTC2013」idiots - 数论+FFT | Bill Yang's Blog
0%

「bzoj3513」「MUTC2013」idiots - 数论+FFT

题目大意

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


题目分析

补集转化。
将题目所求转化为$1-$不能组成三角形的概率。
概率可以通过统计数目进行计算。
而不能组成三角形的情况即两边之和大于第三边。

两边之和的方案数可以用FFT快速计算。
用Hash记录一下第三边,即可统计不合法对数了。

注意FFT因为有序会算重,因此方案数应该除以$2$,如果$i$为偶数的时候,还需要减去$\frac i2+\frac i2$选了同一条边这种情况。


代码

不知道为什么,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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<complex>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>

using namespace std;

#define cp complex<double>

inline const int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(!isdigit(x)) {
if(x=='-')bj=-1;
x=getchar();
}
while(isdigit(x)) {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int maxn=262144+5;
const double pi=acos(-1);

struct FastFourierTransform {
int n,rev[maxn];
cp omega[maxn],iomega[maxn];
void init(int n) {
this->n=n;
for(int i=0; i<n; i++) {
omega[i]=cp(cos(2*pi/n*i),sin(2*pi/n*i));
iomega[i]=conj(omega[i]);
}
int k=log2(n);
for(int i=1; i<=n; i++) {
int t=0;
for(int j=0; j<k; j++)if(i&(1<<j))t|=1<<(k-j-1);
rev[i]=t;
}
}
void transform(cp* a,cp* omega) {
for(int i=0; i<n; i++)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int len=2; len<=n; len<<=1) {
int mid=len>>1;
for(cp* p=a; p!=a+n; p+=len) {
for(int i=0; i<mid; i++) {
cp t=omega[n/len*i]*p[mid+i];
p[mid+i]=p[i]-t;
p[i]+=t;
}
}
}
}
void dft(cp* a) {
transform(a,omega);
}
void idft(cp* a) {
transform(a,iomega);
for(int i=0; i<n; i++)a[i]/=n;
}
} fft;

int t,n,Max,Hash[maxn],sum[maxn],a[maxn],b[maxn];
double up,down;

void Multiply() {
int n=1;
while(n<(Max<<1))n<<=1;
static cp c[maxn];
fill(c,c+n,0);
for(int i=0; i<Max; i++)c[i].real(a[i]);
fft.init(n);
fft.dft(c);
for(int i=0; i<n; i++)c[i]*=c[i];
fft.idft(c);
for(int i=0; i<n; i++)b[i]=round(c[i].real());
}

void Clear() {
Max=0;
memset(Hash,0,sizeof(Hash));
memset(sum,0,sizeof(sum));
memset(a,0,sizeof(a));
}

int main() {
t=Get_Int();
while(t--) {
Clear();
n=Get_Int();
for(int i=1; i<=n; i++) {
int x=Get_Int();
Hash[x]++;
a[x]++;
Max=max(Max,x);
}
for(int i=Max; i>=1; i--)sum[i]=sum[i+1]+Hash[i];
Multiply();
for(int i=0; i<=Max; i++) {
if(!(i&1))b[i]-=Hash[i>>1];
b[i]>>=1;
}
up=down=(double)n*(n-1)*(n-2)/6;
for(int i=0; i<=Max; i++)up-=(double)b[i]*sum[i];
printf("%0.7lf\n",up/down);
}
return 0;
}

姥爷们赏瓶冰阔落吧~