HDU4609 3-idiots
传送门
最近比赛遇到了要用FFT的题，就趁机加强一下。
一句话题意：在\(n\)个数中随便拿出三个数，问能组成三角形的概率。（\(1 \leqslant n \leqslant 10^ 5, 1 \leqslant a_i \leqslant10 ^ 5\)）
看到\(a_i\)的限制，就能想到开一个桶\(num[i]\)表示大小为\(i\)的数的个数。
那么\(f = num ^ 2\)，\(f(i)\)就表示拿两个数，组成大小为\(i\)的数的方案数。但这里面有不合法的和重复的，所以要将\(f(a_i+a_i)\)减\(1\)，以及将所有\(f(i)\)除以\(2\)，因为\(a_x+a_y\)和\(a_y+a_x\)是同一种选法。
上述的\(f(i)\)直接用FFT即可计算，接下来考虑怎么组成三角形。
我们可以枚举长边。
那么对于每一个长边\(a_i\)，可以和\(\sum\limits_{j>a_i} f(j)\)组成三角形，但这样算会多。
一是可能会有一个比\(a_i\)长的边，一个比\(a_i\)短的边，那么方案数就要减去\((i-1)*(n-i)\).
二是可能一个刚好也选\(a_i\)（注意，是选\(a_i\)本身，而选和\(a_i\)相同的数是合法的，而且两个短边不可能同时选到\(a_i\)，因为这种情况再刚求完\(f\)的时候已经去掉了），另一个随便选，因此要再减去\(n-1\).
最后一种就是选了两个都比\(a_i\)大的数，再减去\(C_{n - i}^2\).
总的来说这题难点不在于卷积，而是后面组合数的推导，要不重不漏的考虑到所有情况。
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<queue>
#include<assert.h>
#include<ctime>
using namespace std;
#define enter puts("")
#define space putchar(' ')
#define Mem(a, x) memset(a, x, sizeof(a))
#define In inline
#define forE(i, x, y) for(int i = head[x], y; ~i && (y = e[i].to); i = e[i].nxt)
typedef long long ll;
typedef double db;
const int INF = 0x3f3f3f3f;
const db eps = 1e-8;
const int maxn = 3e5 + 5;
const db PI = acos(-1);
In ll read()
{
ll ans = 0;
char ch = getchar(), las = ' ';
while(!isdigit(ch)) las = ch, ch = getchar();
while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar();
if(las == '-') ans = -ans;
return ans;
}
In void write(ll x)
{
if(x < 0) x = -x, putchar('-');
if(x >= 10) write(x / 10);
putchar(x % 10 + '0');
}
int n, a[maxn];
int len = 1, lim = 0, rev[maxn];
struct Comp
{
db x, y;
In Comp operator + (const Comp& oth)const {return (Comp){x + oth.x, y + oth.y};}
In Comp operator - (const Comp& oth)const {return (Comp){x - oth.x, y - oth.y};}
In Comp operator * (const Comp& oth)const {return (Comp){x * oth.x - y * oth.y, x * oth.y + y * oth.x};}
friend In void swap(Comp &a, Comp& b) {swap(a.x, b.x), swap(a.y, b.y);}
}A[maxn];
In void fft(Comp* a, int flg)
{
for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
for(int i = 1; i < len; i <<= 1)
{
Comp omg = (Comp){cos(PI / i), sin(PI / i) * flg};
for(int j = 0; j < len; j += (i << 1))
{
Comp o = (Comp){1, 0};
for(int k = 0; k < i; ++k, o = o * omg)
{
Comp tp1 = a[k + j], tp2 = o * a[k + j + i];
a[k + j] = tp1 + tp2, a[k + j + i] = tp1 - tp2;
}
}
}
}
ll num[maxn], sum[maxn];
In ll solve()
{
sort(a + 1, a + n + 1);
for(int i = 1; i <= n; ++i) num[a[i]]++;
int N = a[n] + a[n];
while(len <= N) len <<= 1, ++lim;
for(int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1));
for(int i = 0; i <= a[n]; ++i) A[i] = (Comp){num[i], 0};
for(int i = a[n] + 1; i < len; ++i) A[i] = (Comp){0, 0};
fft(A, 1);
for(int i = 0; i < len; ++i) A[i] = A[i] * A[i];
fft(A, -1);
for(int i = 0; i <= N; ++i) num[i] = A[i].x / len + 0.5;
for(int i = 1; i <= n; ++i) num[a[i] + a[i]]--;
for(int i = 1; i <= N; ++i) num[i] /= 2, sum[i] = sum[i - 1] + num[i];
ll ret = -1LL * n * (n - 1);
for(int i = 1; i <= n; ++i)
{
ret += sum[N] - sum[a[i]];
ret -= 1LL * (i - 1) * (n - i);
ret -= (1LL * (n - i) * (n - i - 1) >> 1);
}
return ret;
}
In void init()
{
Mem(num, 0);
len = 1, lim = 0;
}
int main()
{
int T = read();
while(T--)
{
init();
n = read();
for(int i = 1; i <= n; ++i) a[i] = read();
printf("%.7lf\n", 1.0 * solve() / (n - 2) / (n - 1) / n * 6);
}
return 0;
}
