1. 程式人生 > >bzoj3513 [MUTC2013]idiots FFT

bzoj3513 [MUTC2013]idiots FFT

Description


給定n個長度分別為a_i的木棒,問隨機選擇3個木棒能夠拼成三角形的概率。

第一行T(T<=100),表示資料組數。
接下來若干行描述T組資料,每組資料第一行是n,接下來一行有n個數表示a_i。
3≤N≤105,1≤a_i≤105

T<=20
N<=100000

Solution


考慮算不能拼成三角形的概率。我們設f[x]表示兩條邊組成x的方案數,統計比x大的邊有多少直接算就可以了
求f[]可以FFT

感覺壓一波時間就可以進第一頁了呀。。

Code


#include <stdio.h>
#include
<string.h>
#include <algorithm> #include <math.h> #define rep(i,st,ed) for (register int i=st;i<=ed;++i) #define drp(i,st,ed) for (register int i=st;i>=ed;--i) #define fill(x,t) memset(x,t,sizeof(x)) typedef long long LL; const int N=444489; const double pi=3.1415926535897932
; struct com { double real,imag; com operator +(const com &b) const {return (com) {real+b.real,imag+b.imag};} com operator -(const com &b) const {return (com) {real-b.real,imag-b.imag};} com operator *(const com &b) const {return (com) {real*b.real-imag*b.imag,real*b.imag+imag*b.real}
;} com operator /(const double b) const {return (com) {real/b,imag/b};} } a[N]; int rev[N],c[N]; int read() { int x=0,v=1; char ch=getchar(); for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):(v),ch=getchar()); for (;ch<='9'&&ch>='0';x=x*10+ch-'0',ch=getchar()); return x*v; } void FFT(com *a,int len,double f) { for (register int i=0;i<len;i++) if (i<rev[i]) std:: swap(a[i],a[rev[i]]); for (register int i=1;i<len;i<<=1) { com wn=(com) {cos(pi/i),f*sin(pi/i)}; for (register int j=0;j<len;j+=i*2) { com w=(com) {1,0}; for (register int k=0;k<i;k++) { register int tmp=j+k; com u=a[tmp],v=a[tmp+i]*w; a[tmp]=u+v; a[tmp+i]=u-v; w=w*wn; } } } if (f==-1) for (register int i=0;i<len;i++) a[i]=a[i]/len; } int main(void) { freopen("data.in","r",stdin); freopen("myp.out","w",stdout); for (int T=read();T--;fill(c,0)) { int n=read(),mx=0; rep(i,1,n) { int x=read(); c[x]+=1; mx=std:: max(mx,x); } int len,lg; for (len=1,lg=0;len<=mx*2;len<<=1,lg++); for (register int i=0;i<len;++i) rev[i]=(rev[i/2]/2)|((i&1)<<(lg-1)),a[i].real=c[i],a[i].imag=0; FFT(a,len,1); for (register int i=0;i<len;++i) a[i]=a[i]*a[i]; FFT(a,len,-1); for (register int i=0;i<len;i+=2) a[i].real-=1.0*c[i>>1]; for (register int i=0;i<len;++i) a[i]=a[i]/2.0; LL ans=0; drp(i,mx,0) c[i]+=c[i+1]; rep(i,0,len-1) ans+=(LL)(0.5+a[i].real)*c[i]; LL tmp=1LL*n*(n-1)/2*(n-2)/3; printf("%.7lf\n", 1.0*(tmp-ans)/tmp); } return 0; }