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} imes 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
", (double)(tot - ans) / tot);
}
int main() 
{
	cin >> T;
	while (T--)treAKer();
	return AK IOI;
}