BZOJ 3513: [MUTC2013]idiots FFT
正难则反,首先概率就是\(\frac{合法的方案数}{总的方案数}\)=\(\frac{总的方案数-不合法的方案数}{总的方案数}\),
统计不合法的方案数只需要 两个较短的边的长度和\(\le\)较长的边,用t[i]表示长度大于等于i的木棍的数量,f[i]为长度为i的木棍的数量,g[i]表示选出两根木棍组成和为i的方案数,很明显g等于f卷上f,
注意自己不能卷自己,直接在最后让自己卷自己的方案数减去就行了。
然后\(\displaystyle \sum \frac{g[i]}{2} \times t[i]\)就是不合法的方案数了,除2是因为 先选木棍1后选木棍2 和 先选木棍2后选木棍1 是一样的。
bzoj数据比较毒瘤,NTT超时了,得用FFT才能过…
时间复杂度O(n log n).
#include<bits/stdc++.h>
#define DB double
#define LL long long
#define AK 0
#define IOI ;
using namespace std;
int T, n, x, maxx, lim;
LL ans, tot;
const int N = 400010;
const DB PI = acos(-1);
int r[N], t[N];
LL g[N];//f一根的长度为i的方案数 g两根组成i的方案数 t一根长度大于等于i的方案数
inline int read()
{
int res = 0; char ch = getchar(); bool XX = false;
for (; !isdigit(ch); ch = getchar())(ch == '-') && (XX = true);
for (; isdigit(ch); ch = getchar())res = (res << 3) + (res << 1) + (ch ^ 48);
return XX ? -res : res;
}
struct xu
{
DB x, y;
xu(DB X = 0, DB Y = 0) {x = X, y = Y;}
friend xu operator +(const xu &a, const xu &b)
{return (xu) {a.x + b.x, a.y + b.y};}
friend xu operator -(const xu &a, const xu &b)
{return (xu) {a.x - b.x, a.y - b.y};}
friend xu operator *(const xu &a, const xu &b)
{return (xu) {a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x};}
} f[N];
void FFT(xu *A, int lim, int opt)
{
for (int i = 0; i < lim; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
for (int i = 0; i < lim; ++i)
if (i < r[i])swap(A[i], A[r[i]]);
int len;
xu wn, w, x, y;
for (int mid = 1; mid < lim; mid <<= 1)
{
len = mid << 1;
wn = (xu) {cos(PI / mid), opt*sin(PI / mid)};
for (int j = 0; j < lim; j += len)
{
w = (xu) {1, 0};
for (int k = j; k < j + mid; ++k, w = w * wn)
{
x = A[k]; y = A[k + mid] * w;
A[k] = x + y; A[k + mid] = x - y;
}
}
}
}
void YYch()
{
for (int i = 0; i <= lim; ++i)f[i] = g[i] = t[i] = 0;
maxx = 0; lim = 1; ans = tot = 0;
}
inline void treAKer()
{
YYch();
cin >> n;
for (int i = 1; i <= n; ++i)
{
maxx = max(maxx, x = read());
f[x].x++; t[x]++; g[x << 1]--;
}
for (int i = maxx; i >= 1; --i)t[i] += t[i + 1];
while (lim <= (maxx << 1))lim <<= 1;
FFT(f, lim, 1);
for (int i = 0; i < lim; ++i)f[i] = f[i] * f[i];
FFT(f, lim, -1);
for (int i = 0; i < lim; ++i)g[i] += (int)(f[i].x / lim + 0.5);
tot = (LL)n * (n - 1) * (n - 2) / 6;
for (int i = 0; i < lim; ++i)ans += (g[i] >> 1) * t[i];
printf("%.7f\n", (double)(tot - ans) / tot);
}
int main()
{
cin >> T;
while (T--)treAKer();
return AK IOI;
}