3-idiots HDU - 4609 FFT模板

给N个木棍问任选三个可组合出三角形的概率

以前写的，整理整理存个档

//#pragma comment(linker, "/STACK:1024000000,1024000000")
//#include<bits/stdc++.h>
#include<stdio.h>
#include<algorithm>
#include<queue>
#include<string.h>
#include<iostream>
#include<math.h>
#include<set>
#include<map>
#include<vector>
#include<iomanip>
using namespace std;
#define ll long long
#define pb push_back
#define FOR(a) for(int i=1;i<=a;i++)
const double PI = acos(-1.0);

struct complex
{
double r, i; //real and image
complex(double _r = 0, double _i = 0)
{
r = _r; i = _i;
}
complex operator +(const complex &b)
{
return complex(r + b.r, i + b.i);
}
complex operator -(const complex &b)
{
return complex(r - b.r, i - b.i);
}
complex operator *(const complex &b)
{
return complex(r*b.r - i*b.i, r*b.i + i*b.r);
}
};
void change(complex y[], int len)
{
int i, j, k;
for (i = 1, j = len / 2; i < len - 1; i++)
{
if (i < j)swap(y[i], y[j]);
k = len / 2;
while (j >= k)
{
j -= k;
k /= 2;
}
if (j < k)j += k;
}
}
void fft(complex y[], int len, int on)
{
change(y, len);
for (int h = 2; h <= len; h <<= 1)
{
complex wn(cos(-on * 2 * PI / h), sin(-on * 2 * PI / h));
for (int j = 0; j < len; j += h)
{
complex w(1, 0);
for (int k = j; k < j + h / 2; k++)
{
complex u = y[k];
complex t = w*y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w*wn;
}
}
}
if (on == -1)
for (int i = 0; i < len; i++)
y[i].r /= len;
}

const int maxn = 4e5 + 6;
complex x1[maxn];
int arr[maxn >> 2];
ll num[maxn]; ll sum[maxn];

int main() {
int T, n;
scanf("%d", &T);
while (T--) {
memset(num, 0, sizeof num);
scanf("%d", &n);
for (int i = 0; i<n; i++) {
scanf("%d", &arr[i]);
num[arr[i]]++;
}
sort(arr, arr + n);
int len1 = arr[n - 1] + 1; //数域范围
int len = 1;
while (len<2 * len1)len <<= 1; //扩充成最小2的幂次，即多项式项数
for (int i = 0; i<len1; i++)x1[i] = complex(num[i], 0);
for (int i = len1; i<len; i++)x1[i] = complex(0, 0);
fft(x1, len, 1); //应用2n阶的fft计算出A(x)的点值表达
for (int i = 0; i<len; i++)x1[i] = x1[i] * x1[i]; //逐点相乘,O(n)
fft(x1, len, -1); //对2n个点值应用fft，计算其逆dft
for (int i = 0; i<len; i++) {
num[i] = (long long)(x1[i].r + 0.5);
}
len = 2 * arr[n - 1]; //从2的幂次缩回原来那么大
/************************************************************************************************************/

//减掉相同组合
for (int i = 0; i<n; i++)num[arr[i] + arr[i]]--;
//选择无序，地位等价
for (int i = 1; i <= len; i++)num[i] /= 2;
sum[0] = 0;
for (int i = 1; i <= len; i++)sum[i] = sum[i - 1] + num[i];

ll cnt = 0;
for (int i = 0; i<n; i++) {
cnt += sum[len] - sum[arr[i]];
cnt -= (long long)(n - 1 - i)*i;
cnt -= n - 1;
cnt -= (long long)(n - 1 - i)*(n - i - 2) / 2;
}
ll tot = 1ll * n*(n - 1)*(n - 2) / 6;
printf("%.7lf\n", (double)cnt / tot);
}
}