P3338 [ZJOI2014]力(FFT+化简)

传送门

$F_j=\sum _{i=1}^{j-1}\frac{q_i\ast q_j}{(i-j{)}^2}-\sum _{j+1}^n\frac{q_i\ast q_j}{(i-j{)}^2}$

求$E_i=\frac{F_i}{q_i}$

---

$E_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i}{(i-j{)}^2}$

我们回忆一下卷积的形式

$E_i=\sum _{j+k=i}{f}_j\ast g_k$

我们是不是可以凑成这个形式呢??

我们定义$f_i=q_i$,定义$g_i=\frac{1}{i^2}$

那么如果只看前半部分

$\sum _{i=1}^{j-1}{q}_i\ast \frac{1}{(i-j{)}^2}=\sum _{i=1}^{j-1}f[i]\ast g[j-i]$

而且我们令$f[0]=0,g[0]$,那么上下标取到$j,0$也不影响

也就是计算$\sum _{i=0}^{j}f[i]\ast g[j-i]$,很明显这就是卷积的形式

那么再来看后半部分

$\sum _{i=j+1}^n{q}_i\ast \frac{1}{(i-j{)}^2}=\sum _{i=j}^{n}f[i]\ast g[i-j]$

$\sum _{i=0}^{n-j}f[i+j]\ast g[i]$

即使这样$i+i+j!=j$

我们定义$F[i]=f[n-i]$

那么$\sum _{i=0}^{n-j}f[i+j]\ast g[i]=\sum _{i=0}^{n-j}F[n-i-j]\ast g[i]$

其实右边已经是卷积的形式了,为什么呢?因为我们令$t=n-j$

右边就是$\sum _{i=0}^{t}F[t-i]\ast g[i]$

求两次卷积相减即可

#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;
const int maxn = 4e5+10;
const double pi = acos(-1.0);
struct complex
{
	double x,y;
	complex(double xx=0,double yy=0){ x=xx; y=yy; }
}f[maxn],F[maxn],g1[maxn],g2[maxn];
int r[maxn],n;
complex operator - (complex a,complex b){ return complex(a.x-b.x,a.y-b.y); }
complex operator + (complex a,complex b){ return complex(a.x+b.x,a.y+b.y); }
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
void FFT(complex *a,int limit,int type)
{
	for(int i=0;i<limit;i++)	if( i<r[i] )	swap( a[i],a[r[i]]  );
	for(int mid=1;mid<limit;mid<<=1)
	{
		complex wn = complex( cos(pi/mid),type*sin(pi/mid) );
		for(int i=0,R=mid<<1;i<limit;i+=R)
		{
			complex w = complex(1,0);
			for(int k=0;k<mid;k++,w=w*wn)
			{
				complex x=a[i+k], y = w*a[i+k+mid];
				a[i+k] = x+y, a[i+k+mid] = x-y;
			}
		}
	}
}
void print(complex *a)
{
	for(int i=0;i<=n;i++)	cout << a[i].x << " ";
	cout << endl;
}
void mul(complex *a,complex *b,int n,int m)
{
	int limit = 1;
	while( limit<=n+m )	limit<<=1;
	for(int i=0;i<limit;i++)
		r[i] = ( r[i>>1]>>1 )|( (i&1)?limit>>1:0 );
	for(int i=0;i<=limit;i++)
	{
		if( i>n )	a[i].x=a[i].y = 0;
		if( i>m )	b[i].x=b[i].y= 0 ;
	}
	FFT(a,limit,1); FFT(b,limit,1);
	for(int i=0;i<=limit;i++)	a[i] = a[i]*b[i];
	FFT(a,limit,-1);
	for(int i=0;i<=limit;i++)	a[i].x/=limit;
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&f[i].x );
		F[n-i+1].x = f[i].x;
	}
	for(int i=1;i<=n;i++)	g1[i].x = g2[i].x = 1.0/(double)i/(double)i;
	mul(f,g1,n,n); mul(F,g2,n,n);
	for(int i=1;i<=n;i++)
		printf("%.4lf\n",f[i].x-F[n-i+1].x );
}