bzoj 3527: [Zjoi2014]力【FFT】
大力推公式，目标是转成卷积形式：\( C_i=\sum_{j=1}^{i}a_jb_{i-j} \)
首先下标从0开始存，n--
\[F_i=\frac{\sum_{j<i}\frac{q_jq_i}{(j-i)^2}-\sum_{j>i}\frac{q_jq_i}{(j-i)^2}}{q_i}
\]
\[F_i=\sum_{j<i}\frac{q_j}{(j-i)^2}-\sum_{j>i}\frac{q_j}{(j-i)^2}
\]
设
\[a_i=\sum_{j<i}\frac{q_j}{(j-i)^2}
\]
\[b_i=\sum_{j>i}\frac{q_j}{(j-i)^2}
\]
\[f_j=F_j
\]
\[g_i=\frac{1}{i^2}
\]
故
\[F_i=a_i-b_i
\]
先推ai
！注意j==i的情况下g函数为0不影响结果，所以方便起见把大于小于都换成了大于等于小于等于
\[a_i=\sum_{j<i}\frac{q_j}{(j-i)^2}
\]
\[=\sum_{j=0}^{i}\frac{q_j}{(i-j)^2}
\]
\[=\sum_{j=0}^{i}f_i*g_{i-j}
\]
于是卷积*1get
再推bi
\[b_i=\sum_{j>i}\frac{q_j}{(j-i)^2}
\]
\[=\sum_{j=i}^{n}\frac{q_j}{(j-i)^2}
\]
然后用j+i替换j
\[=\sum_{j+i=i}^{n}\frac{q_{j+i}}{(j+i-i)^2}
\]
\[=\sum_{j=0}^{n-i}\frac{q_{j+i}}{j^2}
\]
\[=\sum_{j=0}^{n-i}f_{j+i}*g_j
\]
然后用n-i-j替换j
\[=\sum_{n-i-j=0}^{n-i}f_{n-i-j+i}*g_{n-i-j}
\]
\[=\sum_{0 \leq n-i-j \leq n-i}f_{n-i-j+i}*g_{n-i-j}
\]
然后发现这样转化一下
\[0 \leq n-i-j\Rightarrow j \leq n-i
\]
\[n-i-j \leq n-i\Rightarrow 0 \leq j
\]
\[=\sum_{j=0}^{n-i}f_{n-j}*g_{n-i-j}
\]
设
\[f1_i=f_{n-i}
\]
\[=\sum_{j=0}^{n-i}f1_j*g_{n-i-j}
\]
设
\[t=n-i
\]
\[=\sum_{j=0}^{t}f1_j*g_{t-j}
\]
于是卷积*2get
然后直接上FFT即可
！！！对于g函数，应该是g[i]=1.0/i/i ！！g[i]=1.0/(i*i)会炸精度！
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=1000005;
int n,lm,bt,re[N];
struct cd
{
double r,i;
cd(double R=0,double I=0)
{
r=R,i=I;
}
cd operator + (cd &a) const
{
return cd(r+a.r,i+a.i);
}
cd operator - (cd &a) const
{
return cd(r-a.r,i-a.i);
}
cd operator * (cd &a) const
{
return cd(r*a.r-i*a.i,r*a.i+i*a.r);
}
}f[N],f1[N],g[N],a[N],b[N];
void dft(cd a[],int f)
{
for(int i=1;i<lm;i++)
if(i<re[i])
swap(a[i],a[re[i]]);
for(int i=2;i<=lm;i<<=1)
{
cd wi=cd(cos(2.0*M_PI/i),f*sin(2.0*M_PI/i));
for(int k=0;k<lm;k+=i)
{
cd w=cd(1,0),x,y;
for(int j=0;j<(i>>1);j++)
{
x=a[j+k];
y=a[j+k+(i>>1)]*w;
a[j+k]=x+y;
a[j+k+(i>>1)]=x-y;
w=w*wi;
}
}
}
if(f==-1)
for(int i=0;i<lm;i++)
a[i].r/=lm;
}
int main()
{
scanf("%d",&n);
n--;
for(int i=0;i<=n;i++)
{
scanf("%lf",&f[i].r);
f1[n-i].r=f[i].r;
}
for(int i=1;i<=n;i++)
g[i].r=1.0/i/i;
for(int i=0;;i++)
if((1<<i)>2*n)
{
bt=i;
lm=(1<<i);
break;
}
for(int i=0;i<lm;i++)
re[i]=(re[i>>1]>>1)|((i&1)<<(bt-1));
dft(f,1);
dft(f1,1);
dft(g,1);
for(int i=0;i<lm;i++)
{
a[i]=f[i]*g[i];
b[i]=f1[i]*g[i];
}
dft(a,-1);
dft(b,-1);
for(int i=0;i<=n;i++)
printf("%.3f\n",a[i].r-b[n-i].r);
return 0;
}
