FFT乘法

fft在超大整数乘法中的应用

貌似能把时间从o(n^2) 降低到 o(n logn)

#include <stdio.h>   
#include <string.h>
#include <math.h>
#include <algorithm>
#ifdef __GNUC__
#endif // __GNUC__
using namespace std;
typedef long long ll;
const int MAXN = 400100;
const double PI = acos(-1.0);
struct mcomplex
{
	double r, i;
	mcomplex(double rr=0.0, double ii = 0.0):r(rr),i(ii){}
	mcomplex operator + (mcomplex &a)
	{
		return mcomplex(r+a.r, i+a.i);
	}
	mcomplex operator - (mcomplex &a)
	{
		return mcomplex(r-a.r, i-a.i);
	}
	mcomplex operator * (mcomplex &a)
	{
		return mcomplex(r*a.r-i*a.i, r*a.i+i*a.r);
	}
};
int a[MAXN>>2];
mcomplex xl[MAXN];
ll num[MAXN], sum[MAXN];
void change(mcomplex 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(mcomplex y[], int len, int on)
{
	change(y, len);
	for (int h =2; h <= len; h<<=1)
	{
		mcomplex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
		for (int j = 0; j< len; j+=h)
		{
			mcomplex w(1, 0);
			for (int k = j; k< j+h/2; ++k)
			{
				mcomplex u = y[k];
				mcomplex 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].r /= len;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif // ONLINE_JUDGE
	int t, n, i;
	scanf("%d", &t);
	while (t--)
	{
		scanf("%d", &n);
		memset(num, 0, sizeof num);
		for (i = 0; i< n; ++i)
		{
			scanf("%d", &a[i]);
			num[a[i]]++;
		}
		sort(a, a+n);
		int crest = a[n-1] + 1;
		int len = 1;
		while ( len < 2*crest ) len <<= 1;
		for ( i = 0; i< crest; ++i)			xl[i] = mcomplex(num[i], 0.0);
		for (i = crest ; i< len; ++i)			xl[i] = mcomplex(0, 0);
		fft(xl, len, 1);
		for (i = 0; i< len; ++i)			xl[i] = xl[i]*xl[i];
		fft(xl, len, -1);
		for (i = 0; i< len; ++i)			num[i] = (ll)(xl[i].r + 0.5);
		///.....................................
		len = 2*a[n-1];
		/// 与自身组合
		for (i = 0; i< n; ++i)			--num[a[i]<<1];
		/// 顺序
		for (i = 1; i<= len; ++i)			num[i] >>= 1;
		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[a[i]];/// 假设a[i]是最长的棒子，枚举长度大于a[i]的组合数量 
			cnt -= (ll)(n-1-i)*i; /// 一根大于，一根小于
			cnt -= n-1;	/// 另外两根中不能包括自身
			cnt -= (ll)(n-1-i)*(n-i-2)/2;	/// 两根都大于
		}
		ll tot = (ll) n*(n-1)*(n-2)/6;
		printf("%.7lf\n", (double)cnt/tot);
	}
    return 0;
}