调了好长时间qwq
深刻的教训:
一个 m 次多项式和一个 n 次多项式乘起来是一个 m+n 次多项式而不是 n 次多项式……(sb错误啊啊啊啊)
题面
给定 (n le 10^5), (q_1, cdots, q_n)
((0 < q_i < 10^9) 且为实数),
已知
[E_i = sum_{j=1}^{i-1}frac{q_j}{(i-j)^2} - sum_{j=i+1}^{n}frac{q_j}{(j-i)^2}
]
求
[E_1, E_2, cdots, E_n
]
(当你的输出与标准答案相差不超过 10^{-2}10−2 时即被认为正确。)
首先将 (E_i) 拆成两项
[A_i = sum_{j=1}^{i-1}frac{q_j}{(i-j)^2}
]
[B_i = sum_{j=i+1}^{n}frac{q_j}{(j-i)^2}
]
则 (E_i = A_i - B_i) 。
数列 (A) 可以由数列 (q) 和数列 (1, frac{1}2, frac{1}3, cdots) 卷积得到。
数列 (B) 呢?
[B_i = sum_{j=i+1}^{n}frac{q_j}{(j-i)^2}
]
[=sum_{j-i=1}^{n-i}frac{q_j}{(j-i)^2}
]
[=sum_{s=1}^{n-i}frac{q_{s+i}}{s^2}
]
山穷水尽?
看题解后,
学到许多。
一个小trick
有
[sum_{j=1}^{n-i} f[i+j] *g[j]
]
设 (f'[i] = f[n-i])
则
[sum_{j=1}^{n-i} f[i+j] *g[j]
]
[=sum_{j=1}^{n-i} f'[n-(i+j)] *g[j]
]
[=sum_{j=1}^{n-i} f'[n-i-j] *g[j]
]
令 (t = n-i)
[=sum_{j=1}^{t} f'[t-j] *g[j]
]
爽。
luogu数据AC代码
#include<bits/stdc++.h>
using namespace std;
const int maxn = 2e6+15;
const long double pi = acos(-1.0);
struct cp{
double x,y;
cp() {
x=y=0.0;
}
cp(double x, double y) {
this->x=x, this->y=y;
}
const cp operator+(const cp p)const {
return cp(x+p.x,y+p.y);
}
const cp operator-(const cp p)const {
return cp(x-p.x,y-p.y);
}
const cp operator*(const cp p)const {
return cp(x*p.x-y*p.y, x*p.y+y*p.x);
}
};
int rv[maxn];
void fft(cp *a, int n, int inv) {
for(int i=0;i<n;++i) if(i<=rv[i]) swap(a[i],a[rv[i]]);
for(int m=2;m<=n;m<<=1) {
cp w(cos(2*pi/m), inv*sin(2*pi/m));
for(int i=0;i<n;i+=m) {
cp tmp(1,0);
for(int j=0;j<(m>>1);++j) {
cp p=a[i+j], q=tmp*a[i+j+(m>>1)];
a[i+j]=p+q, a[i+j+(m>>1)]=p-q;
tmp=tmp*w;
}
}
}
if(inv==-1) {
for(int i=0;i<n;++i) a[i].x/=n;
}
}
int n,m;
long double q[maxn];
cp f[maxn], rf[maxn], g[maxn];
int main()
{
scanf("%d", &m);
for(int i=1;i<=m;++i) {
scanf("%lf", &f[i].x);
rf[m-i].x=f[i].x;
}
for(int i=1;i<=m;++i) g[i]=cp(1.0/i/i,0);
// for(int j=1;j<=m;++j) {
// double ans = 0.0;
// for(int i=1;i<=j-1;++i) ans+=q[j]*q[i]*g[j-i].x;
// for(int i=j+1;i<=m;++i) ans-=q[j]*q[i]*g[i-j].x;
// printf("%.3lf
", ans/q[j]);
// } bf algorithm
for(n=1;n<(m<<1);n<<=1);
for(int i=0;i<n;++i) rv[i]=(rv[i>>1]>>1)|((i&1)?(n>>1):0);
fft(f,n,1), fft(rf,n,1), fft(g,n,1);
for(int i=0;i<n;++i) f[i]=f[i]*g[i];
for(int i=0;i<n;++i) rf[i]=rf[i]*g[i];
fft(f,n,-1), fft(rf,n,-1);
for(int i=1;i<=m;++i) printf("%.6lf
", f[i].x-rf[m-i].x);
return 0;
}