[SPOJ TSUM] Triple Sums [FFT]

给一些数，从中选取3个，问他们的和。输出对于每个和，能够组成这个和的方案数。

快速傅立叶变换...第i项的系数为值为i的数出现了多少次，然后计算这个多项式的三次方，再用容斥原理去重即可。用的别人的模板..

SPOJ跑的是真慢...常数大一点都过不去...只好先DFT，然后计算，最后再进行逆变换..

另外FFT丢精度感觉丢的挺多的...本机跑样例，不除6的话，都过不去...

#include <complex>
#include <algorithm>
#include <cmath>
#include <vector>
#include <cstdio>

using namespace std;

typedef complex <long double> Complex;

const long double PI = acos(-1.);
const int N = 1 << 17;
const long double three = 3.0;

void FFT(Complex* P, int n, int oper) {
    for (int i = 1, j = 0; i < n - 1; i++) {
        for (int s = n; j ^= s >>= 1, ~j & s;);
        if (i < j) {
            swap(P[i], P[j]);
        }
    }
    Complex unit_p0;
    for (int d = 0; (1 << d) < n; d++) {
        int m = 1 << d, m2 = m * 2;
        long double p0 = PI / m * oper;
        unit_p0 = Complex(cos(p0), sin(p0));
        for (int i = 0; i < n; i += m2) {
            Complex unit = 1;
            for (int j = 0; j < m; j++) {
                Complex &P1 = P[i + j + m], &P2 = P[i + j];
                Complex t = unit * P1;
                P1 = P2 - t;
                P2 = P2 + t;
                unit = unit * unit_p0;
            }
        }
    }
}

int a[N]={0};
int b[N]={0};
int c[N]={0};
Complex A[N],B[N],AB[N];

int main() {
	int n,x,i;
	scanf("%d",&n);
	for (i=0;i<n;i++) {
		scanf("%d",&x);
		a[x+20000]++;
		b[x+x+40000]++;
		c[x+x+x+60000]++;
	}
	for (i=0;i<N;i++) {
		A[i]=a[i];
		B[i]=b[i];
	}
	FFT(A,N,1);
	FFT(B,N,1);
	for (i=0;i<N;i++) {
		Complex &tmp=A[i];
		AB[i]=tmp*(tmp*tmp-three*B[i]);
	}
	FFT(AB,N,-1);
	for (i=0;i<N;i++) {
		long long ans=((long long)(AB[i].real()/N+0.5)+2*c[i])/6;
		if (ans!=0) printf("%d : %lld\n",i-60000,ans);
	}
	return 0;
}

附上模板...

#include <complex>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;

typedef complex <double> Complex;
typedef vector <int> Polynomial;

const double PI = acos(-1.);
const int N = 1 << 17;

void FFT(Complex* P, int n, int oper) {
	for (int i = 1, j = 0; i < n - 1; i++) {
		for (int s = n; j ^= s >>= 1, ~j & s;);
		if (i < j) {
			swap(P[i], P[j]);
		}
	}
	Complex unit_p0;
	for (int d = 0; (1 << d) < n; d++) {
		int m = 1 << d, m2 = m * 2;
		double p0 = PI / m * oper;
		unit_p0 = Complex(cos(p0), sin(p0));
		for (int i = 0; i < n; i += m2) {
			Complex unit = 1;
			for (int j = 0; j < m; j++) {
				Complex &P1 = P[i + j + m], &P2 = P[i + j];
				Complex t = unit * P1;
				P1 = P2 - t;
				P2 = P2 + t;
				unit = unit * unit_p0;
			}
		}
	}
}

Complex A[N], B[N];

Polynomial operator * (const Polynomial &u, const Polynomial &v)
{
	int n=1,p=u.size(),q=v.size(),i;
	while (n<=p+q-2) n<<=1;
	for (i=0;i<n;++i) A[i]=i<p?u[i]:0;
	for (i=0;i<n;++i) B[i]=i<q?v[i]:0;
	FFT(A,n,1);
	FFT(B,n,1);
	for (i=0;i<n;++i) A[i]*=B[i];
	FFT(A,n,-1);
	Polynomial w(p+q-1);
	for (i=0;i<w.size();++i) 
		w[i]=(int)(A[i].real()/n+0.5);
	return w;
}