bzoj 3513: [MUTC2013]idiots FFT

时间：2021-09-25 02:10:21

bzoj 3513: [MUTC2013]idiots FFT

链接

bzoj

思路

参考了学姐TRTTG的题解

统计合法方案,最后除以总方案。

合法方案要不好统计，统计不合法方案。

\(a+b<=c\)的个数

f[i]是i出现的个数

g[i]表示a+b=i的个数，a<=b

这个可以fft加速到\(nlogn\)统计.

具体的，fft算出ff的卷积，减去自己自己的贡献，然后/2就是了g[i]。

不合法方案数就是:\(\sum f[i]*g[i]\)

最终答案是\(ans=\frac{C_n^3-tot}{C_n^3}\)

代码

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int N=4e5+7;

const double Pi=acos(-1.0);

int read() {

	int x=0,f=1;char s=getchar();

	for(;s>'9'||s<'0';s=getchar()) if(s=='-') f=-1;

	for(;s>='0'&&s<='9';s=getchar()) x=x*10+s-'0';

	return x*f;

}

int n,m,r[N],limit=1,l;

struct Complex {

	double x,y;

	Complex(double xx=0,double yy=0) {x=xx,y=yy;}

}a[N],b[N];

Complex operator + (Complex a,Complex b) {return Complex(a.x+b.x,a.y+b.y);}

Complex operator - (Complex a,Complex b) {return Complex(a.x-b.x,a.y-b.y);}

Complex operator * (Complex a,Complex b) {return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

void fft(Complex *a,int type) {

	for(int i=0;i<=limit;++i)

		if(i<r[i]) swap(a[i],a[r[i]]);

	for(int mid=1;mid<limit;mid<<=1) {

		Complex Wn(cos(Pi/mid),type*sin(Pi/mid));

		for(int i=0;i<limit;i+=mid<<1) {

			Complex w(1,0);

			for(int j=0;j<mid;++j,w=w*Wn) {

				Complex x=a[i+j],y=w*a[i+j+mid];

				a[i+j]=x+y;

				a[i+j+mid]=x-y;

			}

		}

	}

}

ll f[N],g[N],sum[N];

void solve() {

	memset(a,0,sizeof(a));

	memset(f,0,sizeof(f));

	n=read();

	for(int i=1;i<=n;++i) f[read()]++;

	for(int i=0;i<=limit;++i) a[i].x=f[i];

	fft(a,1);

	for(int i=0;i<=limit;++i) a[i]=a[i]*a[i];

	fft(a,-1);

	for(int i=0;i<=limit;++i)

		g[i]=(int)(a[i].x/limit+0.5);

	for(int i=0;i<=limit;++i) g[i*2]-=f[i];

	for(int i=0;i<=limit;++i) g[i]>>=1;

	for(int i=0;i<=limit;++i) g[i]+=g[i-1];

	ll tot=0;

	for(int i=0;i<=limit;++i) tot+=1LL*f[i]*g[i];

	ll ans=1LL*n*(n-1)*(n-2)/6;

	 printf("%.7lf\n",(double)(ans-tot)/ans);

}

int main() {

	while(limit<=200000) limit<<=1,l++;

	for(int i=0;i<=limit;++i)

		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

	for(int T=read();T;T--) solve();

	return 0;

}

相关文章