[BZOJ3513: [MUTC2013]idiots]

BZOJ3513: [MUTC2013]idiots

FFT只是一个工具，重点还是你如何运用。

我们设一个函数\(a(x)\)表示长度为\(x\)的共有\(a(x)\)根木棍；设一个\(f(x)\)表示选出\(2\)根木棍长度和为\(x\)共有\(f(x)\)种方法。显然，\(f(x)=\sum_{y=0}^x a(y)a(x-y)\)。换句话说，\(f=a^2\)。这就是使用FFT的部分。

当然咯，我们还要重复计算的部分。例如，当\(x\)为偶数时，\(f(x)\)就包含了两次都选同一根木棍的不合法情形。因此，这时，\(f(x)\)应该减去\(a(\dfrac{x}{2})\)。

同时，第一次选\(A\)木棍而第二次选\(B\)木棍与第一次选\(B\)木棍而第二次选\(A\)木棍两者本质相同。因此，我们应该把所有的\(f(x)\)都除以\(2\)。

我们设一个\(g(x)\)表示长度\(\geq x\)的木棍共有\(g(x)\)根。则所有不合法的方案数为\(\sum f(x)g(x)\)。

则答案为\(\dfrac{\text{总方案数}-\text{不合法方案数}}{\text{总方案数}}\)。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const double pi=acos(-1);
const double eps=1e-8;
const int lg=18,lim=(1<<lg);
int T,n,rev[lim+5],cnt[lim+5],t[lim+5];
struct cp{
	double x,y;
	cp(double u=0,double v=0){x=u,y=v;}
	friend cp operator +(const cp &u,const cp &v){return cp(u.x+v.x,u.y+v.y);}
	friend cp operator -(const cp &u,const cp &v){return cp(u.x-v.x,u.y-v.y);}
	friend cp operator *(const cp &u,const cp &v){return cp(u.x*v.x-u.y*v.y,u.y*v.x+u.x*v.y);}
}f[lim+5];
void FFT(cp *a,int tp){
    for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int md=1;md<lim;md<<=1){
        cp rt=cp(cos(pi/md),tp*sin(pi/md));
        for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
            cp w=cp(1,0);
            for(int i=0;i<md;i++,w=w*rt){
                cp x=a[pos+i],y=w*a[pos+md+i];
                a[pos+i]=x+y;
                a[pos+md+i]=x-y;
            }
        }
    }
}
ll s[lim+5],up,down;
int main(){
	scanf("%d",&T);
    for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	while(T--){
		scanf("%d",&n),memset(cnt,0,sizeof(cnt)); 
		for(int i=0,x;i<n;i++)scanf("%d",&x),cnt[x]++;
		for(int i=0;i<lim;i++)f[i]=cp(cnt[i],0);
		for(int i=lim-1;i>=0;i--)t[i]=t[i+1]+cnt[i];
		FFT(f,1);
		for(int i=0;i<lim;i++)f[i]=f[i]*f[i];
		FFT(f,-1);
		for(int i=0;i<lim;i++)s[i]=(ll)(f[i].x/lim+0.5);
		for(int i=0;i<lim;i++){
			if(!(i&1))s[i]-=cnt[i>>1];
			s[i]>>=1;
		}
		up=down=(1ll*n*(n-1)*(n-2)/6);
		for(int i=0;i<lim;i++)up-=1ll*s[i]*t[i];
		printf("%.7lf\n",1.0*up/down);		
	}
	return 0;
}