题目大意
给定n个数q1,q2…qn 定义 Fj=j−1∑i=1qi×qj(i−j)2−n∑i=j+1qi×qj(i−j)2
同时 Ei=Fiqi
你的任务是对于所有的1≤i≤n, 求出 Ei
Solution
推下公式: Ej=Fjqj=j−1∑i=1qi(i−j)2−n∑i=j+1qi(i−j)2 然后如果要用FFT加速的话, 要把这个推成卷积式
卷积式就是Ck=∑k0A[i]∗B[k−i]
继续推 Ej=j∑i=1qi(i−j)2−n∑i=jqi(i−j)2 设f[i]=qi, g[i]=1i2
原式变成 Ej=j∑i=1f[i]∗g[j−i]−n∑i=jf[i]∗g[i−j] 然后设f[0]=g[0]=0 Ej=j∑i=0f[i]∗g[j−i]−n∑i=jf[i]∗g[i−j] 发现左边已经是一个卷积式子了, 继续推右边的.
将他展开, 发现 n∑i=jf[i]∗g[i−j]=f[j]∗g[0]+f[j+1]∗g[1]+⋯+f[j+(n−j)]∗g[n−j]=n−j∑i=0f[j+i]∗g[i] 然后发现f还不太行, 翻转一下.
设f′[i]=f[n−i] n−j∑i=0f[j+i]∗g[i]=n−j∑i=0f′[n−(j+i)]∗g[i]=n−j∑i=0f′[n−j−i]∗g[i] 然后设t=n−j
原式 t∑i=0f′[t−i]∗g[i] 答案就是 Ej=j∑i=0f[i]∗g[j−i]−t∑i=0f[t−i]∗g[i](t=n−j) 然后设A(x)=∑ni=0f[i],B(x)=∑ni=0g[i],C(x)=∑ni=0f′[i]
设D(x)=A(x)∗B(x),E(x)=B(x)∗C(x)
然后 Ej=D(j)−E(n−j)
就做完了
代码:
#include <bits/stdc++.h>
using namespace std;
#define int long long
const long double PI = 3.14159265358978323846264338;
const int N = 1e6 + 10;
complex<double> omega[N], omegaInv[N];
int rev[N];
double f[N], g[N], ff[N], res1[N], res2[N];
inline void init(const int n)
{
int k = 0;
while ((1 << k) < n) k++;
for (int i = 0; i < n; i++)
{
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
omega[i] = complex<double>(cos(2 * PI / n * i), sin(2 * PI / n * i));
omegaInv[i] = conj(omega[i]);
}
rev[0] = 0;
}
void fft(complex<double> *a, const int n, const complex<double> *omega)
{
for (int i = 0; i < n; i++)
{
if (i > rev[i]) continue;
swap(a[i], a[rev[i]]);
}
for (int l = 2; l <= n; l *= 2)
{
int m = l / 2;
for (complex<double> *p = a; p != a + n; p += l)
for (int i = 0; i < m; i++)
{
complex<double> t = omega[n / l * i] * p[m + i];
p[m + i] = p[i] - t;
p[i] += t;
}
}
}
inline void dft(complex<double> *a, const int n)
{
fft(a, n, omega);
}
inline void idft(complex<double> *a, const int n)
{
fft(a, n, omegaInv);
for (int i = 0; i < n; i++)
{
a[i] /= n;
}
}
inline void squared(int *a1, const int n1)
{
int n = 1;
while (n < n1 + n1 + 2) n *= 2;
static complex<double> c1[N];
for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
init(n);
dft(c1, n);
for (int i = 0; i < n; i++)
c1[i] *= c1[i];
idft(c1, n);
for (int i = 0; i <= n1 * 2; i++)
a1[i] = static_cast<int>(floor(c1[i].real() + 0.5));
}
inline void multiply(double *a1, const int n1, double *a2, const int n2, double *res)
{
int n = 1;
while (n < n1 + n2 + 2) n *= 2;
static complex<double> c1[N], c2[N];
memset(c1, 0, sizeof(c1));
memset(c2, 0, sizeof(c2));
for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
for (int i = 0; i <= n2; i++) c2[i].real(a2[i]);
init(n);
dft(c1, n);
dft(c2, n);
for (int i = 0; i < n; i++)
c1[i] *= c2[i];
idft(c1, n);
// idft(c2, n);
for (int i = 0; i <= n1 + n2; i++)
res[i] = c1[i].real();
}
signed main()
{
int n;
cin >> n;
for (int i = 1; i <= n; i++)
{
cin >> f[i];
ff[n - i] = f[i];
g[i] = (double)(1.0 / i / i);
}
// for (int i = 1; i <= n; i++)
// cin >> g[i];
multiply(f, n, g, n, res1);
// for (int i = 0; i <= n * 2; i++)
// cout << fixed << setprecision(10) << res1[i] << ' ';
// cout << endl;
// exit(0);
multiply(ff, n, g, n, res2);
for (int i = 1; i <= n; i++)
cout << fixed << setprecision(3) << res1[i] - res2[n - i] << endl;
return 0;
}