Link:
Luogu https://www.luogu.com.cn/problem/P3338
Description
[forall i=1,2,cdots,n
]
,求:
[E_i~=~sumlimits_{j = 1}^{i - 1} frac{q_j}{(i - j)^2}~-~sumlimits_{j = i + 1}^{n} frac{q_j}{(i - j)^2}
]
其中
[n~in~mathbb{N}~cap~[1,10^5],
]
[forall j~in~mathbb{N}~cap~[1,n]~,~quad q_i~in~mathbb{R}~cap~(0,10^9)
]
Solution
怎么说呢,
[sum_{j = 1}^{i - 1} frac{q_j}{(i - j)^2}
]
算是个卷积吧
那右边那个反正直觉就是可以变成卷积
虽然但是,不难发现把 (i) 改成倒过来的马上就得到一个卷积了
(比如说,原来的 (i=n) 就变成 (i=1))
更加具体地,左边是
首先设 (f(j)=q_j) 以及 (g(x)=x^{-2})
[sumlimits_{j=1}^{i-1}f(j)g(i-j)
]
右边是,
[sumlimits_{t=1}^{n-i}frac{q_{i+t}}{t^2}
]
不妨设
[h(x)quad s.t.quad h(n-i-t)=q_{i+t}
]
……至于下标不够怎么处理?
补0就可以啦()
一度有过这样的神奇想法
“为什么输出答案的时候,左边和右边两个卷积不用先求和呢?”
……没救了。卷积啊
经典错误了
就是那个
lim
。。因为卷积所以至少不小于 (2n)
然后因为
for(0; <lim; )
所以实际上 lim
不能比 (2n+1) 小。
Code
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<cctype>
#include<cmath>
using namespace std;
const int MAXN = 4e5 + 10;
const double tpi = acos(-1.0);
int n, lim = 1, loglim = 0;
int rev[MAXN];
struct complex
{
double real, imag;
complex(const double &aug1 = 0, const double &aug2 = 0)
{
real = aug1;
imag = aug2;
}
complex operator * (complex b)
{
return complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real);
}
complex operator + (complex b)
{
return complex(real + b.real, imag + b.imag);
}
complex operator - (complex b)
{
return complex(real - b.real, imag - b.imag);
}
void operator *= (complex b)
{
*this = *this * b;
}
void operator += (complex b)
{
real += b.real;
imag += b.imag;
}
void operator -= (complex b)
{
real -= b.real;
imag -= b.imag;
}
void operator *= (const double &x)
{
real *= x;
imag *= x;
}
} f[MAXN], g[MAXN], h[MAXN], Wn[MAXN][2];
void fft(complex *f, int opt)
{
static complex t;
for (int i = 0; i < lim; ++i)
{
if (rev[i] > i) swap(f[i], f[rev[i]]);
}
for (int n, stats, half_n = 1; half_n < lim; half_n = n)
{
n = half_n << 1;
stats = lim / n;
for (int pos = 0; pos < lim; pos += n)
{
for (int k = 0; k < half_n; ++k)
{
t = Wn[stats * k][opt] * f[pos + half_n + k];
f[pos + half_n + k] = f[pos + k] - t;
f[pos + k] += t;
}
}
}
if (opt == 1)
{
for (int i = 0; i < lim; ++i)
f[i].real /= lim;
}
}
void init()
{
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
{
scanf("%lf", &f[i].real);
h[n - i].real = f[i].real;
g[i].real = 1.0 / i / i;
}
while (lim <= (n << 1))
{
lim <<= 1;
++loglim;
}
for (int i = 0; i < lim; ++i)
{
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (loglim - 1));
}
const double modpi = 2.0 * tpi / lim;
double theta;
for (int i = 0; i < lim; ++i)
{
theta = modpi * i;
Wn[i][0].real = cos(theta);
Wn[i][1].real = +Wn[i][0].real;
Wn[i][0].imag = sin(theta);
Wn[i][1].imag = -Wn[i][0].imag;
}
}
void work()
{
fft(f, 0);
fft(g, 0);
fft(h, 0);
for (int i = 0; i < lim; ++i)
{
f[i] *= g[i];
h[i] *= g[i];
}
fft(f, 1);
fft(h, 1);
}
int main()
{
init();
work();
for (int i = 1; i <= n; ++i)
{
printf("%.3f
", f[i].real - h[n - i].real);
}
system("pause");
return 0;
}