多项式的各种操作

By SemiWaker

目录

• 多项式乘法点这里
• 多项式逆元
• 多项式除法
• 多项式牛顿迭代法
• 多项式对数
• 多项式指数函数
• 多项式幂函数
• 多项式三角函数
• 多项式多点求值
• 多项式快速插值
  

  ---

多项式逆元

定义

A(x)A−1(x)≡1(modxn)

• A(x)：要求逆元的多项式
• A−1(x)：多项式 A(x) 的逆元
• (modxn)：在模xn的意义下。xn是一个多项式，不是一个数。
  模xn的意义类似数的除法：设a=bd+r，且r<b，那么有a≡r(modb)。在多项式中，只要r的界比b的界小就好了。
  多项式除法要用多项式逆元，但是多项式逆元的定义又是多项式除法，这里先介绍多项式逆元，所以除法部分简单略过
  xn比较特殊，因为模xn实际上就是舍弃掉≥n次的所有项，只剩下 0 到 n-1次的项。
  （原因显然∑N−1i=0aixi=xn∑N−n−1i=0aixi+∑n−1i=0aixi）
  为什么要模xn呢？
  因为如果不模的话，除了只有常数项的多项式以外，其他多项式的逆元都是无穷级数。
  例：1+2x+3x2。为了方便我们只写系数数列，从次数小的开始。
  可以写成 1 1 1 。
  设这个数列是ai，然后设逆元的系数数列是bi，最后设它们乘起来是ci。
  根据定义，ci满足c0=1，∀i>0,ci=0。
  然后我们考虑这个的逆元会是什么。
  显然b0=1。
  然后c1=a0b1+a1b0，可以解出b1=−1
  然后c2=a0b2+a1b1+a2b0，可以解出b2=0
  依次类推b3=1，b4=−1……
  发现无论如何都有非0项存在，所以是无穷级数。
  模xn之后就可以解决这个问题，我们就可以只求我们需要的前n项。

做法

暴力做法

上面举例为什么一定要模xn的时候，已经给出了暴力做法了。
直接递推，时间复杂度O(n2)

非暴力：

基本思路：倍增
只有常数项的多项式求逆元是比较简单的事，即普通的求逆元。
然后考虑怎么让已经求出的多项式翻倍，也就是倍增的思路。

首先，我们设长度为⌈n2⌉的已经求出的逆元为B(x)。
注意，是向上取整。
向上取整的好处是：翻倍后保证长度≥n。

A(x)B(x)≡1(modx⌈n2⌉)

然后直接设当前求的逆元为A−1(x)

A(x)A−1(x)≡1(modxn)

注意到，前n项为(1 0 0 0 0 ……），那么前⌈n2⌉项也是(1 0 0 0 ……）。
所以当前的逆元折半后还是逆元。

A(x)A−1(x)≡1(modx⌈n2⌉)

把两个式子合在一起得到

A(x)(B(x)−A−1(x))≡0(modx⌈n2⌉)

显然A(x)不为0，所以后面的为0。

B(x)−A−1(x)≡0(modx⌈n2⌉)

我们现在有了之前的逆元和当前的逆元的关系，但是这时在模⌈n2⌉的情况下的，考虑怎样才能变成在模n的情况下表示之前逆元和当前逆元的关系。

有一个很简单的方法：把它平方了。

(B(x)−A−1(x))2≡0(modxn)

为什么一个在模⌈n2⌉的情况下为0的多项式，平方之后模xn也是0呢？
因为模⌈n2⌉的情况下为0的意思是前⌈n2⌉项都是0。
平方之后，显然前⌈n2⌉项也还是0。
然后我们考虑后面的那些项。
因为是卷积，所以an=∑aian−i。其中 i 和 n-i 必然有一个小于等于⌈n2⌉，那么每一项至少有一个数是0，所以全部都是0。
那么就说明，平方之后的多项式的前n项必然都是0，所以我们可以用平方来倍增。

接着把平方拆开。

B2(x)+A−2(x)−2B(x)A−1(x)≡0(modxn)

A−2(x)会很不好处理，我们整体乘一个A(x)。（显然A(x)非0，所以不会有问题。）
得到：

B2(x)A(x)+A−1(x)−2B(x)≡0(modxn)

这样我们就凑出了现在的逆元和之前的逆元的关系了。

A−1(x)≡B(x)(2−B(x)∗A(x))(modxn)

注意到B(x)只有⌈n2⌉这么多项，所以后面的式子应该最多2n项。
所以把B(x)和A(x)用FFT得到2n项的点值然后计算即可。

时间复杂度

T(n)=T(n2)+O(nlogn)=O(nlogn)

常数大概是FFT的6倍。（每一层3个FFT，加起来要遍历两次。）

代码

void PolyRev(int n,Complex *A,Complex *B){    static Complex tmp[MAXN*4];    if (n==1) {B[0]=1.0/A[0];return;}    PolyRev((n+1)>>1,A,B);    int N=1;while (N<n<<1) N<<=1;    for (int i=0;i<n;++i) tmp[i]=A[i];    for (int i=n;i<N;++i) tmp[i]=0;    FFT(N,tmp,1);FFT(N,B,1);    for (int i=0;i<N;++i) B[i]=B[i]*(2.0-tmp[i]*B[i]);    FFT(N,B,-1);    for (int i=n;i<N;++i) B[i]=0;}

NTT版类似，这里用复数只是方便理解

注意点

• 是的，n不一定是2的幂，因为每次向上取整，所以n是多少都可以。
• 多项式是否有逆元，根据上面的过程，只取决于常数项是否有逆元。
  

  ---

多项式除法

定义

A(x)=B(x)D(x)+R(x)

• A(x)：被除数多项式，界为n。
• B(x)：除数多项式，界为m，显然有m<n。
• D(x)：商多项式，界为n-m+1。
• R(x)：余数多项式，界小于m。
  类似普通的除法中余数要小于除数才能保证结果唯一，虽然多项式无法比较大小，但是R(x)的界比B(x)的界小，就可以保证结果唯一了。

多项式取模

A(x)≡R(x)(modB(x))

即上面除法式子的模的形式，是等价的。

做法

思路：转换
我们观察除法和逆元的差别：
1、有余数，所以B(x)D−1(x)≠1(modxn)
2、没有模，不需要模除法的结果长度也是确定的。
我们知道逆元怎么做，那么我们试着用逆元来做除法。

先考虑求出商，求出商之后余数是可以带进去求出来的。
第二个差别很容易解决，我们只需要强行模xn−m+1就可以保证除法的结果不变。
那么有余数的问题怎么办呢？
能不能在模的时候把余数去掉？
观察每一个多项式的项数：
A(x) n
B(x) m
D(x) n-m+1
R(x) <=m-1
也就是说，在模xn−m+1的情况下，会抛弃掉前m-1项。
但是，R(x) 只有后m-1项，所以直接模无法去掉余数。

考虑怎么让后m-1项移动到前m-1项。
定义反转操作：

AR(x)=xn−1A(1x)=∑i=0nan−i−1xi

也就是，我们把 1x 带入 A(x)，然后乘上xn−1。
这个操作的意义是什么呢？
反转系数。
原来的某一项aixi，变为aix−i×xn−1=aixn−i−1。
也就是系数反转了。

我们试着用系数反转来去掉余数。

AR(x)=xn−1A(1x)=xm−1B(1x)xn−mD(1x)+xn−m+1xm−2R(1x)

也就是：

AR(x)=BR(x)DR(x)+xn−m+1RR(x)

这时再模xn−m+1

AR(x)≡BR(x)DR(x)(modxn−m+1)

然后问题解决。
具体步骤如下：
反转A(x)、B(x)，求BR(x)在模xn−m+1意义下的逆元，然后乘上AR(x)得到DR(x)，再反转过来就得到了商D(x)。
如果要求R(x)，代入R(x)=A(x)−B(x)D(x)即可。

需要的只有逆元和乘法而已。
时间复杂度O(nlogn)

注意点

1、除法没有模，不要写多了。
2、要保证m≤n，否则没有意义。
3、所有计算长度不会超过2n，所以FFT的长度是2n。

程序

void PolyDiv(int n,Complex *A,int m,Complex *B,Complex *D,Complex *R){    static Complex A0[MAXN*4],B0[MAXN*4],B1[MAXN*4];//D    int N=1;while (N<n<<1) N<<=1;    for (int i=0;i<N;++i) A0[i]=B0[i]=B1[i]=D[i]=R[i]=0.0;    for (int i=0;i<n;++i) A0[n-i-1]=A[i];    for (int i=0;i<m;++i) B0[m-i-1]=B[i];    PolyRev(n-m+1,B0,B1);    FFT(N,A0,1);FFT(N,B1,1);    for (int i=0;i<N;++i) D[i]=A0[i]*B1[i];    FFT(N,D,-1);    for (int i=n-m+1;i<N;++i) D[i]=0.0;    for (int i=0,j=n-m;i<j;++i,--j) swap(D[i],D[j]);//R    for (int i=0;i<m;++i) B0[i]=B[i];    for (int i=0;i<n-m+1;++i) B1[i]=D[i];    for (int i=n-m+1;i<N;++i) B1[i]=0.0;    FFT(N,B0,1);FFT(N,B1,1);    for (int i=0;i<N;++i) R[i]=B0[i]*B1[i];    FFT(N,R,-1);    for (int i=0;i<n;++i) R[i]=A[i]-R[i];}

---

多项式牛顿迭代法

定义

普通的牛顿迭代法

问题：求函数f(x)的零点。
首先随便选一个点x0。
然后把f(x)泰勒展开。

f(x)=f(x0)+f′(x0)(x−x0)+f′′(x0)(x−x0)22+⋯

为了简单，只保留线性部分。

f(x)=f(x0)+f′(x0)(x−x0)=0

然后得到

x=x0−f(x0)f′(x0)

当然，这个不是准确解。
只需要用x代替x0，迭代多次后就可以得到近似解。

多项式牛顿迭代法

问题：求函数g(f(x))的零点多项式法f(x)
解释：g(f(x))是一个关于多项式f(x)的函数，其中f(x)是变量。
即求一个多项式f(x)使得

g(f(x))≡0(modxn)

举例

g(f(x))=3f2(x)−2f(x)+1

这个函数的零点即方程3f2(x)−2f(x)+1=0的解，

g(f(x))=f(x)−Ak(x)

这个函数的零点即Ak(x)，也就是把A(x)的k次方求出来。

g(f(x))=lnf(x)−A(x)

这个函数的零点为eA(x)，即A(x)的exp。
简单来说，多项式牛顿迭代法就是用来解关于多项式的方程的。

做法

思路：倍增
当n=1的时候，只有常数项，可以直接算出。
常数项的解可能有多个，代入哪一个决定了最后算出哪个多项式
然后考虑长度⌈n2⌉的多项式f0(x)和长度为n的f(x)的关系。
我们可以把f0(x)当成普通的牛顿迭代法中的x0去泰勒展开，然后观察会怎样：

g(f(x))=g(f0(x))+g′(f0(x))(f(x)−f0(x))+g′′(f0(x))2(f(x)−f0(x))2+⋯

注意到：f0(x)和f(x)的前⌈n2⌉项必然是相同的。
那么，f(x)−f0(x)的前⌈n2⌉项也会是0。
所以(f(x)−f0(x))2以及更高次方的前n项都是0。
所以，对于多项式来说，只保留线性部分得到的是准确解。
然后我们列出对于多项式的牛顿迭代公式：

f(x)=f0(x)−g(f0(x))g′(f0(x))

代入求值即可。

每一层需要逆元和乘法，所以时间复杂度:

T(n)=T(n2)+O(nlogn)=O(nlogn)

常数的话，每一层需要至少3次FFT和一次逆元，相当于9次FFT，所以总共大约18次FFT，即常数至少是FFT的18倍。
几乎等于多了一个logn

注意点

1、多项式函数求导的时候，常数多项式应该视为常数去求导，不应该把常数多项式直接求导。
例如g(f(x))=f2(x)−A(x)
求导得：g′(f(x))=2f(x)
而不是2f(x)−A′(x)
2、FFT的长度是2n。
3、精准度取决于常数项，之后的操作都是完全没有误差的。（如果无视FFT的话。）
4、多个零点的函数，求出哪个零点取决于常数项代入哪个零点。零点的数量也是取决于常数项。

---

多项式对数

定义

对于一个多项式A(x)，求：

lnA(x)(modxn)

即一个n项的多项式。
说实话，多项式对数到底是什么我并不清楚，但是知道功能和对数一样就可以了。

做法

直接通过数学的方法算出。
即

lnA(x)=∫(lnA(x))′=∫A′(x)A(x)

为啥这里又直接把A(x)求导了呢？因为这里的A(x)不是常数。
然后步骤就很简单了：
求出A(x)在模xn意义下的逆元，再求出A′(x)，乘起来之后再积分。

多项式的求导和积分都是可以O(n)完成的。
简单来说：
假设A(x)的系数数列为ai，求导之后的系数数列为bi，积分之后的为ci。
则有：

bi=(i+1)ai+1

ci=ai−1i

记得将空出来的一项设为0就好了。

时间复杂度O(nlogn)

程序

void PolyLn(int n,Complex *A,Complex *B){    static Complex A0[MAXN*4],A1[MAXN*4];    int N=1;while (N<n<<1) N<<=1;    for (int i=0;i<N;++i) A0[i]=0.0,A1[i]=0.0;    PolyRev(n,A,A0);    for (int i=0;i<n-2;++i) A1[i]=A[i+1]*double(i+1);A1[n-1]=0.0;    FFT(N,A0,1);FFT(N,A1,1);    for (int i=0;i<N;++i) A0[i]*=A1[i];    FFT(N,A0,-1);    for (int i=0;i<n-2;++i) B[i+1]=A0[i]/double(i+1);B[0]=0.0;    for (int i=n;i<N;++i) B[i]=0.0;}

---

多项式指数函数

定义

对于一个多项式A(x)，求：

eA(x)(modxn)

即exp(A(x))。

做法

牛顿迭代
首先数学变换

f(x)≡eA(x)(modxn)

两边取ln。

lnf(x)−A(x)≡0(modxn)

定义多项式函数g(f(x))为：

g(f(x))≡lnf(x)−A(x)(modxn)

然后带入牛顿迭代公式

f(x)=f0(x)−lnf0(x)−A(x)1f0(x)

化简得：

f(x)=f0(x)(1−lnf0(x)+A(x))

然后按照牛顿迭代的方法去做就好了。
时间复杂度O(nlogn)

程序

void PolyExp(int n,Complex *A,Complex *B){    if (n==1) {B[0]=exp(A[0]);return;}    static Complex A0[MAXN*4],B0[MAXN*4];    PolyExp((n+1)>>1,A,B);    int N=1;while (N<n<<1) N<<=1;    for (int i=0;i<N;++i) B0[i]=0.0;    PolyLn(n,B,B0);    for (int i=0;i<n;++i) A0[i]=A[i];for (int i=n;i<N;++i) A0[i]=0.0;    FFT(N,A0,1);FFT(N,B,1);FFT(N,B0,1);    for (int i=0;i<N;++i) B[i]=B[i]*(1.0-B0[i]+A0[i]);    FFT(N,B,-1);    for (int i=n;i<N;++i) B[i]=0.0;}

---

多项式幂函数

定义

对于一个多项式A(x)和一个数k，求：

Ak(x)(modxn)

做法

数学变换

Ak(x)≡eklnA(x)(modxn)

只需要先求lnA(x)，再乘上k，再求exp即可。
实际，乘上一个多项式也是可行的。所以用这种方法也可以求出多项式的多项式次方。
时间复杂度O(nlogn)。
当然，常数太大，如果是正整数的话直接快速幂也差不多了。

程序

void PolyPow(int n,Complex *A,Complex *B,double k){    static Complex A0[MAXN*4];    PolyLn(n,A,A0);    for (int i=0;i<n;++i) A0[i]*=k;    PolyExp(n,A0,B);}

---

多项式三角函数

定义

对于一个多项式A(x)，求：

sin(A(x))(modxn)

cos(A(x))(modxn)

求出这两个三角函数之后就可以求出所有三角函数了。
虽然意义不明

做法

求导是没有什么用的，求导之后还是三角函数。
所以需要找一个三角函数和其他初等函数的关系式。
欧拉公式

eix=cos(x)+isin(x)

把A(x)带入x得到

eiA(x)=cos(A(x))+isin(A(x))

这样一次exp就可以求出cos(A(x))和sin(A(x))。
唯一需要注意的问题就是乘了i之后怎么做exp。
首先，多项式的系数要用复数表示。
然后考虑FFT，实际上FFT的过程完全没有变化，直接代入即可。
接着考虑逆元，没有什么区别，常数项的逆元也是1a0，只不过是复数的除法而已。
然后考虑ln，完全没区别。
最后考虑exp，有区别的地方在于常数项，要代入普通的欧拉公式来求。
这时求出来的显然是无理数，所以NTT就不要想了。

---

多项式多点求值

多项式快速插值

留坑待填

---

BySemiWaker