hdu 4609 3-idiots FFT

/*
hdu 4609 3-idiots FFT
http://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 here
	int 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
",(double)cnt/tot);
	}
	return 0;
}


原文地址:https://www.cnblogs.com/riskyer/p/3347988.html