bzoj 3527: [Zjoi2014]力 FFT

首先我们知道(displaystyle E_j=sum_{i<j}frac{q_i}{(i-j)^2}-sum_{i>j}frac{q_i}{(i-j)^2}),

设(displaystyle g[i]=frac{1}{i^2}),因为(g)是偶函数，所以(displaystyle E_j=sum_{i=0}^{j-1} q_i g[j-i]-sum_{i=j+1}^n q_ig[j-i])。

前面这东西很明显就是卷积，处理后面就要发挥人类智慧了，将(q_i) 翻转成新数组(q_i')，原来的式子就成了

(displaystyle E_j=sum_{i=0}^{j-1} q_i g[j-i]-sum_{i=0}^{j-1} q_i'g[j-i]),

后面这东西也变成卷积啦，用FFT处理一下就ok啦。

时间复杂度O(n log n).

#include<iostream>
#include<cstdio>
#include<cmath>
#define DB double
using namespace std;
int n, lim = 1, tmp;
const DB PI = acos(-1);
const int N = 400010;
struct lmaginary 
{
	DB x, y;
	lmaginary(double X = 0, double Y = 0) {x = X, y = Y;}
	friend lmaginary operator +(const lmaginary &a, const lmaginary &b)
	{return (lmaginary) {a.x + b.x, a.y + b.y};}
	friend lmaginary operator -(const lmaginary &a, const lmaginary &b)
	{return (lmaginary) {a.x - b.x, a.y - b.y};}
	friend lmaginary operator *(const lmaginary &a, const lmaginary &b)
	{return (lmaginary) {a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x};}
} f1[N], f2[N], g[N];
DB q[N];
int r[N];
void FFT(lmaginary *A, int lim, int opt) 
{
	for (int i = 0; i < lim; ++i)
		if (i < r[i])swap(A[i], A[r[i]]);
	for (int mid = 1; mid < lim; mid <<= 1) 
	{ //长度的一半
		lmaginary wn(cos(PI / mid), opt * sin(PI / mid));
		for (int len = mid << 1, j = 0; j < lim; j += len) 
		{
			lmaginary w((DB)1, (DB)0);
			for (int k = j; k < mid + j; ++k, w = w * wn) 
			{
				lmaginary a = A[k], b = w * A[k + mid];
				A[k] = a + b;
				A[k + mid] = a - b;
			}
		}
	}
}
void YYCH() 
{
	for (int i = 1; i <= n; ++i) 
	{
		g[i].x = (1.0 / i / i);
		f1[i].x = f2[n - i + 1].x = q[i];
	}
	while (lim < 2 * n)lim <<= 1, ++tmp;
	for (int i = 0; i < lim; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (tmp - 1));
}
int main() 
{
	cin >> n;
	for (int i = 1; i <= n; ++i)scanf("%lf", &q[i]);
	YYCH();
	FFT(f1, lim, 1); FFT(f2, lim, 1); FFT(g, lim, 1);
	for (int i = 0; i < lim; ++i)
		f1[i] = f1[i] * g[i], f2[i] = f2[i] * g[i];
	FFT(f1, lim, -1); FFT(f2, lim, -1);
	for (int i = 1; i <= n; ++i)
		printf("%f
", (f1[i].x - f2[n - i + 1].x) / lim);
	return 0;
}