洛谷P3338 [ZJOI2014]力(FFT)

题意：洛谷P3338

分析：

题目给了我们这样一大坨的式子，$F_j={\sum }_{i<j}\frac{q_i{q}_j}{(i-j{)}^2}-{\sum }_{i>j}\frac{q_i{q}_j}{(i-j{)}^2}$。
然后让我们求$E_j=\frac{F_j}{q_j}$。
首先我们化简一下这个式子，让$F_i$除$q_i$得到$E_j={\sum }_{i<j}\frac{q_i}{(i-j{)}^2}-{\sum }_{i>j}\frac{q_i}{(i-j{)}^2}$，即
$E_j={\sum }_{i=0}^{j-1}\frac{q_i}{(i-j{)}^2}-{\sum }_{i=j+1}^n\frac{q_i}{(i-j{)}^2}$
这个式子暴力计算是$O(n^2)$的，无法接受，考虑$FFT$优化。

令$p(i)=\frac{1}{i^2}$，则上式化为$F_j={\sum }_{i=0}^{j-1}{q}_i\ast p(i-j)-{\sum }_{i=j+1}^n{q}_i\ast p(i-j)$

令$q_i^{\ast }$表示$q_i$的逆置，则${\sum }_{i=j+1}^n{q}_i\ast p(i-j)$可以转化为${\sum }_{i=0}^{j-1}{q}_i^{\ast }\ast p(i-j)$,
逆置有疑问的看文章末尾我画的丑图。

最后因为$i-j$会是负数，而$p(i)$是偶函数，所以最后式子转化为
$E_j={\sum }_{i=0}^{j-1}{q}_i\ast p(j-i)-{\sum }_{i=0}^{j-1}{q}_i^{\ast }\ast p(j-i)$
等号两端都是卷积直接$FFT$求即可

#include <bits/stdc++.h>
using namespace std;

typedef long long LL;
typedef unsigned long long ULL;
typedef pair<int, int> pii;
typedef pair<LL, LL> pll;

const int maxn  = 100000 + 5;
const int maxm  = 100 + 5;
const int inf   = 0x3f3f3f3f;
const LL  mod   = 1e9 + 7;//19260817
const double pi = acos(-1.0);

int n, L, len, rev[maxn << 2];

struct Complex{
	double x, y;
	Complex(double _x = 0, double _y = 0) : x(_x), y(_y){}
	Complex operator + (const Complex &a)const{
		return Complex(x + a.x, y + a.y);
	}
	Complex operator - (const Complex &a)const{
		return Complex(x - a.x, y - a.y);
	}
	Complex operator * (const Complex &a)const{
		return Complex(x * a.x - y * a.y, x * a.y + y * a.x);
	}
}a[maxn << 2], b[maxn << 2], c[maxn << 2];

void FFT(Complex a[], int len, int f){
	Complex x, y, w, wn;
	for(int i = 0; i < len; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1; i < len; i <<= 1){
		wn = Complex(cos(pi / i), f * sin(pi / i));
		for(int j = 0, p = i << 1; j < len; j += p){
			w = Complex(1, 0);
			for(int k = 0; k < i; k++, w = w * wn){
				x = a[j + k], y = w * a[j + k + i];
				a[j + k] = x + y, a[j + k + i] = x - y;
			}
		}
	}
	if(!~f) for(int i = 0; i <= len; i++) a[i].x /= 1.0 * len;
}

int main(){ 
	scanf("%d", &n);
	for(int i = 1; i <= n; i++){
		scanf("%lf", &a[i].x);
		b[n - i + 1].x = a[i].x;
		c[i].x = 1.0 / i / i;
	}
	for(len = 1; len <= n + n; len <<= 1) ++L;
	for(int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
	FFT(a, len, 1);
	FFT(b, len, 1);
	FFT(c, len, 1);
	for(int i = 1; i <= len; i++){
		a[i] = a[i] * c[i];
		b[i] = b[i] * c[i];
	}
	FFT(a, len, -1);
	FFT(b, len, -1);
	for(int i = 1; i <= n; i++) printf("%.3lf\n", a[i].x - b[n - i + 1].x);
    return 0; 
}

$FFT$入门。