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);
	}
}