[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{ ext{总方案数}- ext{不合法方案数}}{ ext{总方案数}})。

代码：

#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
",1.0*up/down);		
	}
	return 0;
}