HDU4609 计数问题+FFT

时间:2023-03-09 01:50:53
HDU4609 计数问题+FFT

题目大意:

给出n(1e5)条线段(长度为整数,<=1e5),求任取3条线段能组成一个三角形的概率。

用cnt[i]记录长度为i的线段的个数,通过卷积可以计算出两条线段长度之和为i的方案数sum[i]:先用FFT计算出cnt[i]的卷积sum[i],为取两条线段长度和为i的排列数(包括自己和自己),去掉自己和自己的方案数,再对所有sum[i]除以2即为所求方案数。

之后对所有线段a[i]有大到小排列,考虑第i条线段是三角形最长边的情况(长度相同则将编号大的视为更长,就没有长度相同的情况了。实际计算时无影响)。首先是sum[i + 1] + sum[i + 2] + ... + sum[len - 1],即其他两边之和要大于他的方案数。然后去掉:1.另两边有一条比他大,有一条比他小。 2.两条都比他长。 3.有一条是他自己。这三种情况。

要注意有些地方爆int。还有len的计算要注意边界(一直wa,wa了半天)。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std; const double pi = acos(-); struct Complex{
double real, imag;
Complex (){}
Complex (double _real, double _imag){
real = _real;
imag = _imag;
}
}; Complex operator + (Complex x, Complex y){
return Complex(x.real + y.real, x.imag + y.imag);
} Complex operator - (Complex x, Complex y){
return Complex(x.real - y.real, x.imag - y.imag);
} Complex operator * (Complex x, Complex y){
return Complex(x.real * y.real - x.imag * y.imag, x.real * y.imag + x.imag * y.real);
} Complex operator / (Complex x, double y){
return Complex(x.real / y, x.imag / y);
} int reverse(int x, int len){
int t = ;
for (int i = ; i < len; i <<= ){
t <<= ;
if (x & i) t |= ;
}
return t;
} Complex A[];
void FFT(Complex *a, int n, int DFT){
for (int i = ; i < n; ++i) A[reverse(i, n)] = a[i];
for (int i = ; i <= n; i <<= ){
Complex wn = Complex(cos( * pi / i), DFT * sin( * pi / i));
for (int j = ; j < n; j += i){
Complex w = Complex(, );
for (int k = ; k < (i >> ); ++k){
Complex x = A[j + k];
Complex y = w * A[j + k + (i >> )];
A[j + k] = x + y;
A[j + k + (i >> )] = x - y;
w = w * wn;
}
}
}
if (DFT == -) for (int i = ; i < n; ++i) A[i] = A[i] / n;
for (int i = ; i < n; ++i) a[i] = A[i];
} int T, n;
int a[];
int cnt[];
Complex B[];
long long sum[]; int main(){ scanf("%d", &T);
while (T--){
scanf("%d", &n);
int maxL = ;
memset(cnt, , sizeof(cnt));
memset(sum, , sizeof(sum));
for (int i = ; i < n; ++i){
scanf("%d", a + i);
if (a[i] > maxL) maxL = a[i];
++cnt[a[i]];
}
int len = ;
while (len <= maxL * ) len <<= ;
for (int i = ; i <= maxL; ++i) B[i] = Complex(cnt[i], );
for (int i = maxL + ; i < len; ++i) B[i] = Complex(, );
FFT(B, len, );
for (int i = ; i < len; ++i) B[i] = B[i] * B[i];
FFT(B, len, -);
for (int i = ; i < len; ++i) sum[i] = (long long)(B[i].real + 0.5);
for (int i = ; i < n; ++i) --sum[a[i] + a[i]];
for (int i = ; i < len; ++i) sum[i] >>= ;
for (int i = len - ; i >= ; --i) sum[i] += sum[i + ];
sort(a, a + n);
long long ans = ;
for (int i = ; i < n; ++i){
long long tmp = sum[a[i] + ];
tmp -= ((long long)n - i - ) * i;
tmp -= ((long long)n - i - ) * (n - i - ) / 2LL;
tmp -= n - ;
ans += tmp;
}
double Ans = (double)ans * 6.0 / n / (n - ) / (n - );
printf("%.7f\n", Ans);
} return ;
}