hdu 4609 3-idiots [fft 生成函数 计数]

hdu 4609 3-idiots

题意：

给出(A_i)，问随机选择一个三元子集，选择的数字构成三角形的三边长的概率。

---

一开始一直想直接做....

先生成函数求选两个的方案(注意要减去两次选择同一个的，然后/2)，然后统计三角形个数。

枚举三角形最长边，求(i+j>k,i<k,j<k)的方案数。后两个条件减去不合法的。

不合法很好统计 (i ge k ightarrow i+j > k)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N = (1<<18) + 5;
const double PI = acos(-1.0);
inline int read() {
    char c=getchar(); int x=0,f=1;
    while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
    return x*f;
}

struct meow{
    double x, y;
    meow(double a=0, double b=0):x(a), y(b){}
};
meow operator +(meow a, meow b) {return meow(a.x+b.x, a.y+b.y);}
meow operator -(meow a, meow b) {return meow(a.x-b.x, a.y-b.y);}
meow operator *(meow a, meow b) {return meow(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);}
meow conj(meow a) {return meow(a.x, -a.y);}
typedef meow cd;

namespace fft {
	int n, rev[N];
	cd omega[N], omegaInv[N];
	void init(int lim) {
		n = 1; int k = 0; while(n < lim) n <<= 1, k++;
		for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));

		for(int i=0; i<n; i++) {
			omega[i] = cd(cos(2*PI/n*i), sin(2*PI/n*i));
			omegaInv[i] = conj(omega[i]);
		}
	}

	void dft(cd *a, int flag) {
		cd *w = flag == 1 ? omega : omegaInv;
		for(int i=0; i<n; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
		for(int l=2; l<=n; l<<=1) {
			int m = l>>1;
			for(cd *p = a; p != a+n; p += l) 
				for(int k=0; k<m; k++) {
					cd t = w[n/l*k] * p[k+m];
					p[k+m] = p[k] - t;
					p[k] = p[k] + t;
				}
		}
		if(flag == -1) for(int i=0; i<n; i++) a[i].x /= n;
	}
}

int n, q[N]; ll s[N];
cd a[N];
int main() {
	freopen("in", "r", stdin);
	int T = read();
	while(T--) {
		n = read(); int m = 0;
		for(int i=1; i<=n; i++) q[i] = read(), m = max(m, q[i]);
		fft::init(m+m+1);
		memset(a, 0, sizeof(a));
		for(int i=1; i<=n; i++) a[q[i]].x ++;

		fft::dft(a, 1);
		for(int i=0; i<fft::n; i++) a[i] = a[i] * a[i];
		fft::dft(a, -1);

		for(int i=1; i <= m<<1; i++) s[i] = (ll) floor(a[i].x + 0.5);
		for(int i=1; i <= n; i++) s[q[i]<<1] --;
		for(int i=1; i <= m<<1; i++) s[i] = s[i-1] + (s[i]>>1);
		sort(q+1, q+1+n);
		ll ans = 0, tot = (ll) n * (n-1) * (n-2) / 6;
		for(int i=1; i<=n; i++) {
			//ll t = s[m<<1] - s[q[i]] - (ll) (n-i) * (i-1) - (n-1) - (ll) (n-i) * (n-i-1) / 2;
			ll t = s[m<<1] - s[q[i]] - (ll) (n-i+1) * (n-1) + (ll) (n-i+1) * (n-i) / 2;
			ans += t;
		}
		printf("%.7lf
", (double) ans / tot);
	}
}