读书人

hdu 4609 三-idiots FFT

发布时间: 2013-11-08 17:52:01 作者: rapoo

hdu 4609 3-idiots FFT

/*hdu 4609 3-idiots FFThttp://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html*/#pragma warning(disable : 4786)#pragma comment(linker, "/STACK:102400000,102400000")#include <iostream>#include <string>#include <cstring>#include <cstdlib>#include <cstdio>#include <cmath>#include <algorithm>#include <functional>#include <map>#include <set>#include <vector>#include <stack>#include <queue>//priority_queue#include <bitset>#include <complex>#include <utility>/***make_pair()**/using namespace std;const double eps=1e-8;const double PI=acos(-1.0);const int inf=0x7fffffff;template<typename T> inline T MIN(T a,T b){return a<b?a:b;}template<typename T> inline T MAX(T a,T b){return a>b?a:b;}template<typename T> inline T SQR(T a){return a*a;}inline int dbcmp(double a){return a>eps?(1):(a<(-eps)?(-1):0);}typedef __int64 LL;int n;const int size=400040;int data[size/4];LL num[size],sum[size];complex<double> x1[size];void change(complex<double> y[],int len){    int i,j,k;    for(i = 1, j = len/2;i < len-1;i++)    {        if(i < j)swap(y[i],y[j]);        k = len/2;        while( j >= k)        {            j -= k;            k /= 2;        }        if(j < k)j += k;    }}void fft(complex<double> y[],int len,int on){    change(y,len);    for(int h = 2;h <= len;h <<= 1)    {        complex<double> wn(cos(-on*2*PI/h),sin(-on*2*PI/h));        for(int j = 0;j < len;j += h)        {            complex<double> w(1,0);            for(int k = j;k < j+h/2;k++)            {                complex<double> u = y[k];                complex<double> t = w*y[k+h/2];                y[k] = u+t;                y[k+h/2] = u-t;                w = w*wn;            }        }    }    if(on == -1)        for(int i = 0;i < len;i++)            y[i].real(y[i].real()/len);}int main() {// your code goes hereint t,i;cin>>t;while(t--){cin>>n;memset(num,0,sizeof(num));for(i=0;i<n;++i){cin>>data[i];num[data[i]]++;}sort(data,data+n);int len1=data[n-1]+1;int len=1;while(len<2*len1) len<<=1;for(i=0;i<len1;++i){x1[i]=complex<double>(num[i],0);}for(;i<len;++i){x1[i]=complex<double>(0,0);}fft(x1,len,1);for(i=0;i<len;++i){x1[i]=x1[i]*x1[i];}fft(x1,len,-1);for(i=0;i<len;++i){num[i]=(LL)(x1[i].real()+0.5);}len=2*data[n-1];for(i=0;i<n;++i){num[data[i]+data[i]]--;}for(i=1;i<=len;++i){num[i]/=2;}sum[0]=0;for(i=1;i<=len;++i){sum[i]=sum[i-1]+num[i];}LL cnt=0;for(i=0;i<n;++i){cnt+=sum[len]-sum[data[i]];//减掉一个取大一个取小的情况,违背了data[i]是最长边的规则cnt-=(LL)(n-1-i)*i;//减掉一个取本身,一个取其他数的情况,否则data[i]就被取了两次cnt-=(LL)(n-1);//减掉取两个比他大的数的情况,违背了data[i]是最长边的规则cnt-=(LL)(n-1-i)*(n-2-i)/2;}LL tot = (LL)n*(n-1)*(n-2)/6;printf("%.7lf\n",(double)cnt/tot);}return 0;}

读书人网 >编程

热点推荐