[ZJOI2014]力 FFT
题面
题解:
\[F_j = \sum_{i < j}\frac{q_iq_j}{(i - j)^2} - \sum_{i > j}{\frac{q_iq_j}{(i - j)^2}}
\]
\[E_j = \sum_{i < j}\frac{q_i}{(i - j)^2} - \sum_{i > j}{\frac{q_i}{(i - j)^2}}
\]
对式子的2个部分分别计算。
令\(S_i = i^2\)
\[\sum_{i < j}\frac{q_i}{(i - j)^2} = \sum_{i < j}q_i S_{j - i}
\]
看上去就是卷积形式,FFT计算即可。
对于后半部分,将序列翻转,\(i > j\)就变成\(i < j\)了,而\(S\)可以看做距离,所以不会变,直接计算就好了.
计算完之后需要将序列翻转回来
#include<bits/stdc++.h>
using namespace std;
#define R register int
#define ld double
#define LL long long
#define AC 310000
const double pi = acos(-1);
int n, lim = 1, len;
int Next[AC];
ld q[AC], f[AC];
struct node{
ld x, y;
node(ld xx = 0, ld yy = 0){x = xx, y = yy;}
}a[AC], b[AC], s[AC];
node operator * (node x, node y){return node(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);}
node operator + (node x, node y){return node(x.x + y.x, x.y + y.y);}
node operator - (node x, node y){return node(x.x - y.x, x.y - y.y);}
inline int read()
{
int x = 0;char c = getchar();
while(c > '9' || c < '0') c = getchar();
while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
return x;
}
void FFT(node *A, int opt)
{
for(R i = 0; i < lim; i ++)
if(i < Next[i]) swap(A[i], A[Next[i]]);
for(R i = 1; i < lim ; i <<= 1)
{
node W(cos(pi / i), opt * sin(pi / i));
for(R r = i << 1, j = 0; j < lim; j += r)
{
node w(1, 0);
for(R k = 0; k < i ; k ++, w = w * W)
{
node x = A[j + k], y = w * A[j + k + i];
A[j + k] = x + y, A[j + k + i] = x - y;
}
}
}
}
void pre()
{
n = read() - 1;
for(R i = 0; i <= n; i ++) scanf("%lf", &q[i]);
while(lim <= n + n) lim <<= 1, ++ len;
for(R i = 0; i <= lim; i ++)
Next[i] = (Next[i >> 1] >> 1) | ((i & 1) << (len - 1));
}
void work()
{
for(R i = 0; i <= n; i ++)
{
if(i != 0) s[i].x = 1.0 / i / i;
a[i].x = q[i], b[n - i].x = q[i];
}
FFT(s, 1);
FFT(a, 1);
for(R i = 0; i < lim; i ++) a[i] = a[i] * s[i];
FFT(a, -1);
for(R i = 0; i < lim; i ++) f[i] = a[i].x / lim;
FFT(b, 1);
for(R i = 0; i < lim; i ++) b[i] = b[i] * s[i];
FFT(b, -1);
for(R i = 0; i <= n; i ++)
if(i < n - i) swap(b[i], b[n - i]);
for(R i = 0; i < lim; i ++) f[i] -= b[i].x / lim;
for(R i = 0; i <= n; i ++) printf("%.3lf\n", f[i]);
}
int main()
{
//freopen("in.in", "r", stdin);
pre();
work();
//fclose(stdin);
return 0;
}