hdu 4609 3-idiots 【FFT快(gui)速傅立叶变换】

标签：style   blog   os   io   for   2014   art   cti   

FFT实现起来挺复杂的，开始用vector，结果发现空间超了，换成数组还是超，删掉了几个后又超时了

sin cos 函数调用次数太多了，改成乘法，还是超时

最后把FFT里的除法运算和模运算优化了一下，终于过了，排的老后面

坑，3843MS，一看时间最少的只有671MS，我都怀疑这是不是同一个算法。。为毛差距这么大

#pragma comment(linker, "/STACK:102400000,102400000")
#include<iostream>
#include<vector>
#include<algorithm>
#include<cstdio>
#include<queue>
#include<stack>
#include<string>
#include<map>
#include<set>
#include<cmath>
#include<cassert>
#include<cstring>
#include<iomanip>
#include<ctime>
using namespace std;
#ifdef _WIN32
typedef __int64 i64;
#define out64 "%I64d\n"
#define in64 "%I64d"
#else
typedef long long i64;
#define out64 "%lld\n"
#define in64 "%lld"
#endif
/************ for topcoder by zz1215 *******************/
#define foreach(c,itr)  for(__typeof((c).begin()) itr=(c).begin();itr!=(c).end();itr++)
#define FOR(i,a,b)      for( int i = (a) ; i <= (b) ; i ++)
#define FF(i,a)         for( int i = 0 ; i < (a) ; i ++)
#define FFD(i,a,b)      for( int i = (a) ; i >= (b) ; i --)
#define S64(a)          scanf(in64,&a)
#define SS(a)           scanf("%d",&a)
#define LL(a)           ((a)<<1)
#define RR(a)           (((a)<<1)+1)
#define pb              push_back
#define pf              push_front
#define X               first
#define Y               second
#define CL(Q)           while(!Q.empty())Q.pop()
#define MM(name,what)   memset(name,what,sizeof(name))
#define MC(a,b)        memcpy(a,b,sizeof(b))
#define MAX(a,b)        ((a)>(b)?(a):(b))
#define MIN(a,b)        ((a)<(b)?(a):(b))
#define read            freopen("in.txt","r",stdin)
#define write           freopen("out.txt","w",stdout)

const int inf = 0x3f3f3f3f;
const i64 inf64 = 0x3f3f3f3f3f3f3f3fLL;
const double oo = 10e9;
const double eps = 10e-9;
const double pi = acos(-1.0);

const int maxn = 401111;

struct complex {
	double a, b;
	complex(double _a = 0.0, double _b = 0.0) :a(_a), b(_b){}
	complex operator + (const complex &c) const {
		return complex(a + c.a, b + c.b);
	}
	complex operator - (const complex &c) const {
		return complex(a - c.a, b - c.b);
	}
	complex operator * (const complex &c) const {
		return complex(a*c.a - b*c.b, a*c.b + b*c.a);
	}
};

complex w[maxn];
int id[maxn];  // change id

complex fac[maxn];  // vector a
complex fbc[maxn];  // vector b
complex fabc[maxn]; // vector a*b

void make_id(int id[],int len)
{
	for (int step = len / 2; step >= 1; step >>= 1) {
		for (int it = len / step / 2; it < len / step; it++){
			id[it] = id[it - len / step / 2] + step;
		}
	}
}

void dft(complex a[], complex tc[] , complex w[], int id[],int len) // 变换函数 a ，临时函数 tc， 复数矩阵 w ，长度 len 
{
	for (int i = 0; i < len; i++){
		tc[i] = a[id[i]];
	}
	int temp,tw;
	for (int step = 1; step <= len / 2; step <<= 1){
		for (int x = 0; x < len; x += step *2){
			temp = len / (step * 2);
			tw = 0;
			for (int it = x; it < x + step; it++){
				a[it] = tc[it] + tc[it + step] * w[tw];
				tw += temp;
			}
			for (int it = x + step; it < x + step * 2; it++){
				a[it] = tc[it - step] + tc[it] * w[tw];
				tw += temp;
			}
		}
		for (int i = 0; i < len; i++){
			tc[i] = a[i];
		}
	}
}

void convolution(complex a[], complex b[],complex c[],int len)   //卷积
{
	for (int i = 0; i < len; i++){
		c[i] = a[i] * b[i];
	}
}

void divide(complex c[], int len)
{
	for (int i = 0; i < len; i++){
		c[i].a /= len;
	}
}


void fft(complex a[], complex b[], complex c[],int id[],complex w[], int len)
{
	make_id(id,len);
	double ang = double(2.0*pi) / double(len); 
	w[1] = complex(cos(ang), sin(ang));
	w[0] = complex(1, 0);
	for (int x = 2; x < len; x++){
		w[x] = w[1] * w[x - 1];
	}
	dft(a, c, w, id, len);
	dft(b, c, w, id, len);
	convolution(a, b, c, len);

	w[1] = complex(cos(-ang), sin(-ang));
	for (int x = 2; x < len; x++){
		w[x] = w[1] * w[x - 1];
	}
	dft(c, a, w, id, len);
	divide(c, len);
}

void tocomplex(int a[],complex ca[],int len)  
{
	for (int i = 0; i < len; i++){
		ca[i] = complex(a[i], 0);
	}
}

void toint(i64 a[], complex ac[], int len)
{
	for (int i = 0; i < len; i++){
		a[i] = i64(ac[i].a + 0.5);
	}
}

int n;
int branche[maxn];
int fx[maxn];
i64 res[maxn];

double start()
{
	int temp = 0;
	for (int i = 1; i <= n; i++){
		temp = max(temp, branche[i]);
	}
	temp += 10;
		
	int len = 1;
	while (len < temp){
		len <<= 1;
	}
	len <<= 1;

	for (int i = 0; i < len; i++){
		fx[i] = 0;
	}
	for (int i = 1; i <= n; i++){
		fx[branche[i]]++;
	}

	tocomplex(fx, fac, len);
	tocomplex(fx, fbc, len);
	fft(fac, fbc, fabc, id, w, len);
	toint(res, fabc, len);

	for (int i = 1; i <= n; i++){
		res[2 * branche[i]] --;
	}
	for (int i = 0; i <len; i++){
		res[i] >>= 1;
	}
	for (int i = 1; i < len; i++){
		res[i] += res[i - 1];
	}
	i64 num = 0;
	for (int i = 1; i <= n; i++){
		num += res[branche[i]];
	}
	i64 tot = (i64)n*(n - 1)*(n - 2);
	tot /= 6;
	double ans = double(num) / double(tot);
	ans = 1.0 - ans;
	return ans;
}


int main()
{
	int T;
	cin >> T;
	while (T--){
		cin >> n;
		for (int i = 1; i <= n; i++){
			//cin >> branche[i];
			SS(branche[i]);
		}
		printf("%.7lf\n", start());
	}
	return 0;
}

hdu 4609 3-idiots 【FFT快(gui)速傅立叶变换】,布布扣,bubuko.com

hdu 4609 3-idiots 【FFT快(gui)速傅立叶变换】

标签：style   blog   os   io   for   2014   art   cti   

原文地址：http://blog.csdn.net/zz_1215/article/details/38341381