多项式初步与傅里叶变换（Fourier Transform）

“虚妄与实冢的交错，逆塑变换之坤仑。”

虚数 / 复数

虚数的概念与实数相对，抽象为：在数轴中表示不出的数。

复（）数是实数和虚数的统称。

虚数单位元

在初中数学学习中，我们知道 在 的时候无意义。

那现在，我们定义 ，满足 ，任意虚数可以表示为 。

则复数可以表示为 ，此处定义 为实数域的任意实数。

---

复数的几何表示

对于任意复数 ，我们可以将实部视作 轴，虚部视作 轴，则任意复数我们都可以将其视作向量，表示为 ，其与 轴的夹角被称为复数的幅角。复数的模 。

这个平面被称为复平面。

---

复平面

轴表示实数， 轴除原点外表示虚数，四象限表示复数。

幅角：以逆时针为正方向，从 正半轴到已知复数向量的转角的有向角称为该向量的幅角。

---

复数的运算

复数的运算遵循实数性运算。

两个复数的乘积的幅角为两个复数的幅角和，两个复数的乘积的模长为两个复数的模长乘积。

---

复数意义的单位根

在复平面作单位圆，等分为 份，取其 份记为 ，称为 次单位根。表示夹角大小是 。

我们来研究其性质：

• 对于 ，当 时除外。

这些性质会在我们求解 的时候用到。

---

多项式

被表示为 ，这是一个 次多项式，如果我们把 抽象为二维平面的点，则有性质：

• 用任一 个不同点均可唯一确定一个 次多项式。

上述性质为多项式的点表示法。

点表示法的证明（选修）

我们将这 个点表示为 ，那我们可以得到 个方程。

$$\begin{cases} a_0+a_1x_1+\cdots+a_n{x_1}^n&=y_1 \\ a_0+a_1x_2+\cdots+a_n{x_2}^n&=y_2 \\ a_0+a_1x_3+\cdots+a_n{x_3}^n&=y_3 \\ &\vdots \\ a_0+a_1x_n+\cdots+a_n{x_n}^n&=y_n \\ a_0+a_1x_{n+1}+\cdots+a_n{x_{n+1}}^n&=y_{n+1} \end{cases} \nonumber$$

这是一个 元一次多项式。

这个多项式有唯一解等价于其系数矩阵的秩是满秩的，同样等价于其系数矩阵的行列式不为 。

其行列式表示为：

$$\begin{pmatrix} 1 & x_1 & {x_1}^2 & \cdots & {x_1}^n \\ 1 & x_2 & {x_2}^2 & \cdots & {x_2}^n \\ \vdots &\vdots&\vdots&\ddots&\vdots \\ 1 & x_{n+1} & {x_{n+1}}^2 & \cdots & {x_{n+1}}^n \end{pmatrix} \nonumber$$

这个行列式会在线性代数里学习到，称为范德蒙行列式。是一个 阶范德蒙行列式。

这个行列式的值为 ，所以有唯一解。

这个性质在复数范围内同样适用。

---

点表示法的运用

能够通过 个点在线性时间内求出这个多项式，运用为拉格朗日插值。

也可用于多项式乘积，对于 次多项式 ，如果我们求 ，暴力乘积消耗 得到一个 次的多项式。

而如果我们使用点表示法：

得到 个点表示为：

以及 同理。

那么我们可以在线性时间内求出 表示为：

这是否正确对应了多项式 ，我们之后再谈。但这便是 （快速傅里叶变换）的快速计算卷积的思路。

关于点表示法与系数表示法的实际操作，可以参考以下文章：

• 「洛谷日报：快速傅里叶变换」：https://www.luogu.com.cn/blog/retepviolet/kuai-su-fu-li-xie-bian-huan

---

快速傅里叶变换

观前提醒

以及多项式全家桶（ 等）并不在 提高组考察大纲内，建议省选选手食用。

根据点表示法，我们可以思考，如何用快速的时间求得 ，这样也是可以求出 的。

分治

这里我们假定 ，并有 。

我们将 按照奇偶性分类得到：

那么可以得到 。

我们将 代入，替换成 ，得到：

进行化简得到：

此时 。

那对于后半段，我们代入 得到：

化简得到：

此时 。

---

我们将 的前后两段表达式合在一起看：

$$A(x)= \begin{cases} A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^kA_2(\omega_{\frac{n}{2}}^k)&,k\in[0,\frac{n}{2}-1] \\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^kA_2(\omega_{\frac n 2}^k)&,k\in[\frac{n}{2},n-1] \end{cases} \nonumber$$

容易发现，前后半段所需要的表达式是相同的，只不过是 的区别。这意味着，我们可以通过 个表达式来求出所有的 。

这便是 的分治思想。

---

变换

但容易发现，多项式的系数表示法计算时间复杂度都是 的，所以我们考虑在点表示法的意义下进行运算。

对于点 ，我们设定 。为表示简洁，我们设 ，设一个 次多项式 。

则有 ，且有 。

我们可以来证明 。

$$\begin{align} c_k&= \sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i \nonumber \\ &=\nonumber \sum_{i=0}^{n-1}\left(\sum_{j=0}^{n-1}a_j\times(\omega_n^i)^j \right)\left(\omega_n^{-k}\right)^i \\ &=\nonumber \sum_{i=0}^{n-1}\left(\sum_{j=0}^{n-1}a_j\left(\omega_n^j\right)^i\left(\omega_n^{-k}\right)^i \right) \\ &=\nonumber \sum_{i=0}^{n-1}\left(\sum_{j=0}^{n-1}a_j\left(\omega_n^{j-k}\right)^i \right) \\ &=\nonumber \sum_{j=0}^{n-1}a_j\left(\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i \right) \end{align}$$

设定 ，我们分类讨论：

1. 则有

   即 ，因为 ，则

2. 有

那么继续化简上述式子：

$$\begin{align} c_k&=\nonumber \sum_{j=0}^{n-1}a_j\left(\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i \right) \\ &=\nonumber \sum_{j=0}^{n-1}a_j\times\sum_{i=0}^{n-1} S\left(\omega_n^{j-k}\right) \\ &=\nonumber \sum_{j=0}^{n-1}a_jn[j=k] \\ &=\nonumber na_k \end{align}$$

证毕。

---

实现

在 C++ 中，一般而言，递归常数比迭代更大。所以 一般采用迭代方式实现。

我们考虑 的分治体现：

$$\begin{align} \underline{a_0a_1a_2a_3}&\mid\underline{a_4a_5a_6a_7}& \nonumber \\ \underline{a_0a_2}\mid\underline{a_4a_6}&\mid\underline{a_1a_3}\mid\underline{a_5a_7} \nonumber \\ \nonumber \underline{a_0}\mid\underline{a_4}\mid\underline{a_2}\mid\underline{a_6}&\mid\underline{a_1}\mid\underline{a_5}\mid\underline{a_3}\mid\underline{a_7} \end{align}$$

而我们需要通过下面一行推出上面一行的值，理解为归并。所以我们要预处理出最下面一行乱序的顺序。

如果我们用二进制位表示 的下标就会发现：

$$\begin{align} \underline{000,001,010,011}&\mid\underline{100,101,110,111}& \nonumber \\ \underline{000,010}\mid\underline{100,110}&\mid\underline{001,011}\mid\underline{101,111} \nonumber \\ \nonumber \underline{000}\mid\underline{100}\mid\underline{010}\mid\underline{110}&\mid\underline{001}\mid\underline{101}\mid\underline{011}\mid\underline{111} \end{align}$$

容易发现，最后的排序就是二进制位的翻转，记为 ，其递推公式为：

其中 。

---

题目简介

题目名称：【模板】多项式乘法（FFT）
题目来源：洛谷
评测链接：https://www.luogu.com.cn/problem/P3803

形式化题意，给定一个 次多项式 和一个 次多项式 ，求其卷积。即乘积。

数据范围：

参考代码

// ----- Eternally question-----
// Problem: P3803 【模板】多项式乘法（FFT）
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3803
// Memory Limit: 500 MB
// Time Limit: 2000 ms
// Written by: Eternity
// Time: 2023-01-10 09:21:49
// ----- Endless solution-------

const int MAXN=3e6+10;
const double pi=acos(-1);
int N,M;
struct Complex
{
	double x,y;
	Complex operator+(const Complex &a) const
	{ return {x+a.x,y+a.y}; }
	Complex operator-(const Complex &a) const
	{ return {x-a.x,y-a.y}; }
	Complex operator*(const Complex &a) const
	{ return {x*a.x-y*a.y,x*a.y+y*a.x}; }
}A[MAXN],B[MAXN];
int Rev[MAXN],Bit,Tot;
inline void FFT(Complex a[],int inv)
{
	for(int i=0;i<Tot;++i)
		if(i<Rev[i]) std::swap(a[i],a[Rev[i]]);
	for(int mid=1;mid<Tot;mid<<=1)
	{
		auto w1=Complex({cos(pi/mid),inv*sin(pi/mid)});
		for(int i=0;i<Tot;i+=mid*2)
		{
			auto wk=Complex({1,0});
			for(int j=0;j<mid;++j,wk=wk*w1)
			{
				auto x=a[i+j],y=wk*a[i+j+mid];
				a[i+j]=x+y,a[i+j+mid]=x-y;
			}
		}
	}
}
int main()
{
	// freopen(".in","r",stdin);
	// freopen(".out","w",stdout);
	read(N,M);
	for(int i=0;i<=N;++i) scanf("%lf",&A[i].x);
	for(int i=0;i<=M;++i) scanf("%lf",&B[i].x);
	while((1<<Bit)<N+M+1) ++Bit;
	Tot=1<<Bit;
	for(int i=0;i<Tot;++i) Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(Bit-1));
	FFT(A,1),FFT(B,1);	//正变换
	for(int i=0;i<Tot;++i) A[i]=A[i]*B[i];
	FFT(A,-1);		//逆变换
	for(int i=0;i<=N+M;++i)
		printf("%d ",(int)(A[i].x/Tot+0.5));
	return 0;
}
/*

*/

解决的是实数域上的多项式卷积，所以在使用 C++ 进行整数域卷积时难免会造成精度缺失，所以需要四舍五入，可以用 或者 std::round(A[i].x) 解决。

---

FFT 优化高精乘

题目简介

题目名称： 升级版

题目来源：

评测链接：https://www.luogu.com.cn/problem/P1919

形式化题意：求 。

数据范围：

容易发现，我们原初 的竖式乘法已经不能满足当前这道题的要求了，所以我们考虑能不能优化成 。

对于一个 位数 ，也可以写作：

如果设定 ，那就有：

这不就是个多项式的形式吗，则我们需要求的 ，实际上就是求一个多项式 。

解决了。

AC Code

// ----- Eternally question-----
// Problem: P1919 【模板】A*B Problem 升级版（FFT 快速傅里叶变换）
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P1919
// Memory Limit: 256 MB
// Time Limit: 1500 ms
// Written by: Eternity
// Time: 2023-01-10 09:57:36
// ----- Endless solution-------

#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
	x=0;
	char ch=getchar(),t=0;
	while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
	while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1){ read(x),read(x1...); }
template<class T>
inline void write(T x)
{
	if(x<0) putchar('-'),x=-x;
	if(x>9) write(x/10);
	putchar(x%10+'0');
}
template<>
inline void write(bool x){ putchar(x?'1':'0'); }
template<>
inline void write(char c){ putchar(c); }
template<>
inline void write(char *s){ while(*s!='\0') putchar(*s++); }
template<>
inline void write(const char *s){ while(*s!='\0') putchar(*s++); }
template<class T,class ...T1>
inline void write(T x,T1 ...x1){ write(x),write(x1...); }
template<class T>
inline bool checkMax(T &x,T y){ return x<y?x=y,1:0; }
template<class T>
inline bool checkMin(T &x,T y){ return x>y?x=y,1:0; }
const int MAXN=3e6+10;
const double pi=acos(-1);
int N,M;
struct Complex
{
	double x,y;
	Complex operator+(const Complex &a) const
	{ return {x+a.x,y+a.y}; }
	Complex operator-(const Complex &a) const
	{ return {x-a.x,y-a.y}; }
	Complex operator*(const Complex &a) const
	{ return {x*a.x-y*a.y,x*a.y+y*a.x}; }
}A[MAXN],B[MAXN];
int Rev[MAXN],Bit,Tot;
inline void FFT(Complex a[],int inv)
{
	for(int i=0;i<Tot;++i)
		if(i<Rev[i]) std::swap(a[i],a[Rev[i]]);
	for(int mid=1;mid<Tot;mid<<=1)
	{
		auto w1=Complex({cos(pi/mid),inv*sin(pi/mid)});
		for(int i=0;i<Tot;i+=mid*2)
		{
			auto wk=Complex({1,0});
			for(int j=0;j<mid;++j,wk=wk*w1)
			{
				auto x=a[i+j],y=wk*a[i+j+mid];
				a[i+j]=x+y,a[i+j+mid]=x-y;
			}
		}
	}
}
char aStr[MAXN],bStr[MAXN];
int ans[MAXN];
int main()
{
	// freopen(".in","r",stdin);
	// freopen(".out","w",stdout);
	scanf("%s%s",aStr,bStr);
	N=strlen(aStr)-1,M=strlen(bStr)-1;
	for(int i=0;i<=N;++i) A[i].x=aStr[N-i]-'0';
	for(int i=0;i<=M;++i) B[i].x=bStr[M-i]-'0';
	while((1<<Bit)<N+M-1) ++Bit;
	Tot=(1<<Bit);
	for(int i=0;i<Tot;++i) Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(Bit-1));
	FFT(A,1),FFT(B,1);
	for(int i=0;i<Tot;++i) A[i]=A[i]*B[i];
	FFT(A,-1);
	int Idx=0;
	for(int i=0,x=0;i<Tot||x;++i)
	{
		x+=A[i].x/Tot+0.5;
		ans[++Idx]=x%10;x/=10;
	}
	while(Idx>2&&!ans[Idx-1]) --Idx;
	for(int i=Idx-1;i>=1;--i) write(ans[i]);
	return 0;
}
/*

*/

---

离散傅里叶（逆）变换

OI意义

法如其名，所谓离散，便是通过一系列方法来将一个复杂的多项式拆散成为很多很多很小很小的问题，就如我们在 中提到的一样，通过分治的方法解决系数表示法到点表示法的变换，便是一种离散，我们原本要考虑 个系数，但我们将其离散为 ，然后又将其离散为 ，直到很多个 ，很多对 。这就是离散。

• 离散傅里叶变换：即快速傅里叶变换的正变换。
• 离散傅里叶逆变换：即快速傅里叶变换的逆变换。

---

数学意义

观前提醒

本段与 无关，分类于离散数学，信号分析学。如有并不需要该段知识者，请跳过。

傅里叶分析方法是信号分析的最基本方法，傅里叶变换是傅里叶分析的核心，通过它把信号从时间域变换到频率域，进而研究信号的频谱结构和变化规律。

傅里叶变换集

傅里叶变换，表示能将满足一定条件的某个函数表示成三角函数（正弦和或余弦函数）或者它们的积分的线性组合。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。最初傅里叶分析是作为热过程的解析分析的工具被提出的。

设 ，其傅里叶变换 为 上的函数，定义式如下：

$$\begin{align} \hat f(\kappa)&=<f,E_\kappa> \nonumber \\&=\nonumber \int_{\mathbb T^n}f(x)e^{-2\pi i\kappa\cdot x}dx \end{align}$$

其 被称为傅里叶级数。

离散傅里叶变换是在时域和频域都呈现离散的傅里叶变换形式，将时域信号的采样变换为在离散时间傅里叶变换（）频域的采样。其变换两端（时域和频域）的序列都是有限长度的，这两组序列被认为是离散周期信号的主值序列。即使对有限长的离散信号作 ，也应当将其看作经过周期延拓成为周期信号再作变换。

傅里叶变换在虚轴上，所以它只和虚轴有关。傅里叶变换是积分形式的连续的函数内积，离散傅里叶变换是求和形式的内积。

一般情况下，可以用 快速计算 。

然后的知识就不在我可以理解的范围内了。

推荐文献：

• ：多项式与生成函数 - 快速傅里叶变换 - 离散傅里叶变换
• ：深入理解离散傅里叶变换（）
• 博客园：离散傅里叶变换（）

---

例题

怎么说呢， 在数学板块大概只是一个辅助工具，用来优化你推出来的式子，优化多项式乘法，优化高精度乘法之类的，真正的考点当然也就不会在 上，但显然，你不会 你也拿不到分。

多项式乘法

题目简介

题目名称：多项式乘法
题目来源：安徽

评测链接：https://www.luogu.com.cn/problem/P2553

形式化题意：给定一个字符串，可能是一个不定次多项式，也有可能是两个不定次方程式的卷积。如果是后者，算出该卷积并以字符串形式输出。多测。

数据范围： ，不保证升序

是个人都看得出来可以用 解决，所以难点就在于如何将字符串转换为复数了。当然，次数小于 其实暴力卷积应该也可以过。

AC Code

// ----- Eternally question-----
// Problem: P2553 [AHOI2001]多项式乘法
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P2553
// Memory Limit: 125 MB
// Time Limit: 1000 ms
// Written by: Eternity
// Time: 2023-01-10 11:04:14
// ----- Endless solution-------

#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
	x=0;
	char ch=getchar(),t=0;
	while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
	while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1){ read(x),read(x1...); }
template<class T>
inline void write(T x)
{
	if(x<0) putchar('-'),x=-x;
	if(x>9) write(x/10);
	putchar(x%10+'0');
}
template<>
inline void write(bool x){ putchar(x?'1':'0'); }
template<>
inline void write(char c){ putchar(c); }
template<>
inline void write(char *s){ while(*s!='\0') putchar(*s++); }
template<>
inline void write(const char *s){ while(*s!='\0') putchar(*s++); }
template<class T,class ...T1>
inline void write(T x,T1 ...x1){ write(x),write(x1...); }
template<class T>
inline bool checkMax(T &x,T y){ return x<y?x=y,1:0; }
template<class T>
inline bool checkMin(T &x,T y){ return x>y?x=y,1:0; }
const int MAXN=3e6+10;
const double pi=acos(-1);
int N,M;
char Str[MAXN];
struct Complex
{
	double x,y;
	void clear(){ x=y=0; }
	Complex operator+(const Complex &a) const
	{ return {x+a.x,y+a.y}; }
	Complex operator-(const Complex &a) const
	{ return {x-a.x,y-a.y}; }
	Complex operator*(const Complex &a) const
	{ return {x*a.x-y*a.y,x*a.y+y*a.x}; }
}A[MAXN],B[MAXN];
int Rev[MAXN],Bit,Tot;
inline void FFT(Complex a[],int inv)
{
	for(int i=0;i<Tot;++i)
		if(i<Rev[i]) std::swap(a[i],a[Rev[i]]);
	for(int mid=1;mid<Tot;mid<<=1)
	{
		auto w1=Complex({cos(pi/mid),inv*sin(pi/mid)});
		for(int i=0;i<Tot;i+=mid*2)
		{
			auto wk=Complex({1,0});
			for(int j=0;j<mid;++j,wk=wk*w1)
			{
				auto x=a[i+j],y=wk*a[i+j+mid];
				a[i+j]=x+y,a[i+j+mid]=x-y;
			}
		}
	}
}
inline void transform(int l,int r,Complex a[],int &n)
{
	bool s=0;
	for(int i=l,frt=0,bef=0;i<=r;++i)
	{
		if(Str[i]>='0'&&Str[i]<='9')
		{
			if(s) frt=(frt<<3)+(frt<<1)+(Str[i]^48);
			else  bef=(bef<<3)+(bef<<1)+(Str[i]^48);
		}
		else
		{
			if(Str[i]=='a') continue;
			if(Str[i]=='^') s^=1;
			if(Str[i]=='+'||Str[i]==')')
			{
				a[frt].x=bef,checkMax(n,frt);
				frt=bef=0,s^=1;
			}
		}
	}
}
inline void debug()
{
	printf("Tot:%d Bit:%d\n",Tot,Bit);
	printf("A:sum is %d\n",N);
	for(int i=0;i<=N;++i) printf("%.0lf ",A[i].x);
	printf("\nB:sum is %d\n",M);
	for(int i=0;i<=M;++i) printf("%.0lf ",B[i].x);
	puts("");
}
int main()
{
	// freopen(".in","r",stdin);
	// freopen(".out","w",stdout);
	while(~scanf("%s",Str+1))
	{
		for(int i=0;i<=Tot;++i) A[i].clear();
		for(int i=0;i<=Tot;++i) B[i].clear();
		N=M=Bit=0;
		int Len=strlen(Str+1);
		int posmid=-1;
		for(int i=1;i<=Len;++i)
			if(Str[i]=='*'){ posmid=i;break; }
		if(posmid==-1)
		{
			for(int i=2;i<Len;++i) write(Str[i]);
			puts("");continue;
		}
		transform(1,posmid-1,A,N);
		transform(posmid+1,Len,B,M);
		while((1<<Bit)<=N+M) ++Bit;
		Tot=1<<Bit;
		// debug();
		for(int i=0;i<Tot;++i) Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(Bit-1));
		FFT(A,1),FFT(B,1);
		for(int i=0;i<Tot;++i) A[i]=A[i]*B[i];
		FFT(A,-1);
		for(int i=N+M;i>=1;--i)
			printf("%da^%d+",(int)(A[i].x/Tot+0.5),i);
		write((int)(A[0].x/Tot+0.5),'\n');
	}
	return 0;
}
/*

*/

---

力

题目简介

题目名称：力
题目来源：浙江省选 第二题

评测链接 ：https://www.luogu.com.cn/problem/P3338
评测链接 ：https://loj.ac/p/2200

形式化题意：给定一个长度为 的序列 ，定义：

$$\begin{align} F(j)&=\nonumber \sum_{i=1}^{j-1}\frac{ q_i\times q_j }{ (i-j)^2 } - \sum_{i=j+1}^{n}\frac{ q_i\times q_j }{ (i-j)^2 } \\ E(i)&=\nonumber \frac{F(i)}{q_i} \end{align}$$

对于 ，求 。

数据范围：

我们可以继续推式子：

$$\begin{align} E(i)&=\nonumber \frac{F(i)}{q_i} \\&=\nonumber \sum_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum_{j=i+1}^{n}\frac{q_j}{(i-j)^2} \\&=\nonumber \sum_{j=1}^{i}\frac{q_j}{(i-j)^2}-\sum_{j=i}^{n}\frac{q_j}{(i-j)^2} \end{align}$$

设 并规定 ，代入得到：

$$\begin{align} E(x)&=\nonumber \sum_{i=1}^{x}f(i)\times g(x-i)-\sum_{i=x}^{n}f(i)\times g(i-x) \\&=\nonumber \sum_{i=0}^{x}f(i)\times g(x-i)-\sum_{i=x}^{n}f(i)\times g(i-x) \end{align}$$

化简的意义所在

由于 适用于求解多项式卷积，所以我们需要把式子写成卷积形式。得到：

$$C(x)=\sum_{i=0}^{x}A(i)\times B(x-i)$$

就可以在 的时间内求解。

对于上述式子，左半边已经是一个卷积形式了，所以我们考虑如何来优化右半边。

记 ，将式子展开得到：

所以，我们设一个 得到：

$$\begin{align} C(x)&=\nonumber \sum_{i=0}^{n-x}f'\Big(n-(x+i)\Big)\times g(i) \\&=\nonumber \sum_{i=0}^{n-x}f'(n-x-i)\times g(i) \\&=\nonumber \sum_{i=0}^{m}f'(m-i)\times g(i),m=n-x \end{align}$$

那这样的话，右半边的式子也变成了卷积。

总结一下：

$$E(x)=\sum_{i=0}^{x}f(i)\times g(x-i)-\sum_{i=0}^{m}f'(m-i)\times g(i),m=n-x$$

那推导过程就基本结束了，考虑设三个多项式分别代表：

有 ，则答案为 。

然后注意被卡精度。

AC Code

// ----- Eternally question-----
// Problem: P3338 [ZJOI2014]力
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3338
// Memory Limit: 125 MB
// Time Limit: 3000 ms
// Written by: Eternity
// Time: 2023-01-10 11:47:18
// ----- Endless solution-------

#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
	x=0;
	char ch=getchar(),t=0;
	while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
	while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1){ read(x),read(x1...); }
template<class T>
inline void write(T x)
{
	if(x<0) putchar('-'),x=-x;
	if(x>9) write(x/10);
	putchar(x%10+'0');
}
template<>
inline void write(bool x){ putchar(x?'1':'0'); }
template<>
inline void write(char c){ putchar(c); }
template<>
inline void write(char *s){ while(*s!='\0') putchar(*s++); }
template<>
inline void write(const char *s){ while(*s!='\0') putchar(*s++); }
template<class T,class ...T1>
inline void write(T x,T1 ...x1){ write(x),write(x1...); }
template<class T>
inline bool checkMax(T &x,T y){ return x<y?x=y,1:0; }
template<class T>
inline bool checkMin(T &x,T y){ return x>y?x=y,1:0; }
const int MAXN=3e5+10;
const double pi=acos(-1);
int N,M;
struct Complex
{
	double x,y;
	// Complex(){ x=y=0; }
	void clear(){ x=y=0; }
	Complex operator+(const Complex &a) const
	{ return {x+a.x,y+a.y}; }
	Complex operator-(const Complex &a) const	
	{ return {x-a.x,y-a.y}; }
	Complex operator*(const Complex &a) const
	{ return {x*a.x-y*a.y,x*a.y+y*a.x}; }
}A[MAXN],B[MAXN],C[MAXN];
int Rev[MAXN],Tot,Bit;
inline void FFT(Complex a[],int inv)
{
	for(int i=0;i<Tot;++i)
		if(i<Rev[i]) std::swap(a[i],a[Rev[i]]);
	for(int mid=1;mid<Tot;mid<<=1)
	{
		auto w1=Complex({cos(pi/mid),inv*sin(pi/mid)});
		for(int i=0;i<Tot;i+=mid*2)
		{
			auto wk=Complex({1,0});
			for(int j=0;j<mid;++j,wk=wk*w1)
			{
				auto x=a[i+j],y=a[i+j+mid]*wk;
				a[i+j]=x+y,a[i+j+mid]=x-y;
			}
		}
	}
	// if(inv==-1) for(int i=0;i<Tot;++i) a[i].x/=Tot;
}
int main()
{
	// freopen(".in","r",stdin);
	// freopen(".out","w",stdout);
	read(N);
	for(int i=1;i<=N;++i)
	{
		scanf("%lf",&A[i].x);
		C[N-i].x=A[i].x,B[i].x=(double)(1.0/i/i);
	}
	while((1<<Bit)<=(N<<1)) ++Bit;
	Tot=1<<Bit;
	for(int i=0;i<Tot;++i) Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(Bit-1));
	FFT(A,1),FFT(B,1),FFT(C,1);
	for(int i=0;i<Tot;++i) A[i]=A[i]*B[i],C[i]=C[i]*B[i];
	FFT(A,-1),FFT(C,-1);
	for(int i=1;i<=N;++i) printf("%.3lf\n",(A[i].x-C[N-i].x)/Tot);
	return 0;
}
/*

*/

---

补充讲解及多项式全家桶：

• 多项式卷积与数论变换
• 多项式分治