bzoj3527: [Zjoi2014]力 卷积+FFT

时间：2022-01-09 08:14:50

先写个简要题解：本来去桂林前就想速成一下FFT的，结果一直没有速成成功，然后这几天断断续续看了下，感觉可以写一个简单一点的题了，于是就拿这个题来写，之前式子看着别人的题解都不太推的对，然后早上6点多推了一个多小时终于发现了一个很巧妙的方法，首先问题的关键在于后半个式子，因为显然前半个式子很容易想到卷积的形式，那么直接FFT就好了，但是后半部分不好考虑，一般肯定是通过类似换元的做法化到后来得出结论，到中间有一步就有点难度，那个地方我一直卡。后来突然想到，既然前半部分i<j时那么好处理，那么i>j的情况我把i和j分别用(n-i)和(n-j)代入不就转化过去了，然后直接就会发现他是卷积后的第（n-i)项，所以后面一半需要反转a数组和反转结果数组，这里要注意下标，第二部分，第0项对应第n-1项。

$E_i=\sum_{j<i} \frac{q_j}{(j-i)^2}-\sum_{j>i} \frac{q_j}{(j-i)^2}$
$E_i=\sum_{j<i} \frac{q_j}{(j-i)^2}-\sum_{n-j<n-i} \frac{q_{n-j}}{((n-j)-(n-i))^2}$
$E_i=\sum_{j<i} \frac{q_j}{(i-j)^2}-\sum_{n-j<n-i} \frac{q_{n-j}}{(j-i)^2}$

$E_i=f[j]\ast g[i-j]-f[n-j]\ast g[j-i]=c[i]-c^{'}[n-i]$

代码：

#include<bits/stdc++.h>

#include<complex>

#define db double

#define ll long long

#define mp make_pair

#define fi first

#define pb push_back

#define se second

#define rep(i,a,b)for(int i=a;i<=b;i++)

using namespace std;

const double pi=acos(-1);

const int maxn=5e5+7;

db spt(db x){return x*x;}

int r[maxn];

complex<double>a[maxn],b[maxn],c[maxn],d[maxn];

db ans[maxn];

int M,N,n,l=0;

db p[maxn];

void FFT(complex<double> f[],int op)

{

	for(int i=0;i<N;i++)if(i<r[i])swap(f[i],f[r[i]]);

	for(int i=1;i<N;i<<=1)

	{

		complex<double >w(cos(pi/i),op*sin(pi/i));

		for(int p=(i<<1),j=0;j<N;j+=p)

		{

			complex<double>W(1,0);

			for(int k=0;k<i;k++,W*=w)

			{

				complex<double>x=f[j+k],y=W*f[j+i+k];

				f[j+k]=x+y;f[j+k+i]=x-y;

			}

		}

	}

}

int main()

{

	scanf("%d",&n);

	for(int i=1;i<=n;i++)

	{

		scanf("%lf",&p[i]);

	}

	N=M=n-1;

	for(int i=0;i<n;i++){a[i]=p[i+1];if(i>=1)b[i]=1.0/spt(i);}

	M+=N;

	for(N=1;N<=M;N<<=1)l++;

	for(int i=0;i<N;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

	FFT(a,1);FFT(b,1);

	for(int i=0;i<=N;i++)a[i]=a[i]*b[i];

	FFT(a,-1);

	for(int i=0;i<n;i++)ans[i]+=a[i].real()*1.0/N;

	for(int i=0;i<n;i++){c[i]=p[n-i];if(i>=1)d[i]=1/spt(i);}

	FFT(c,1);FFT(d,1);

	for(int i=0;i<=N;i++)c[i]=c[i]*d[i];

	FFT(c,-1);

	for(int i=0;i<n;i++)ans[i]-=c[n-i-1].real()*1.0/N;

	for(int i=0;i<n;i++)printf("%.3lf\n",ans[i]);

	return 0;

}

　　PS:终于用上latex了，好不适应，，还不会用。

相关文章