多项式的各种操作

版权声明:神犇们抄蒟蒻的总结有意思吗? https://blog.csdn.net/semiwaker/article/details/73251486

多项式的各种操作

By SemiWaker

目录

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

多项式逆元

定义

A(x)A1(x)1(modxn)

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

做法

暴力做法

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

非暴力:

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

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

A(x)B(x)1(modxn2)

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

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

A(x)A1(x)1(modxn2)

把两个式子合在一起得到
A(x)(B(x)A1(x))0(modxn2)

显然A(x)不为0,所以后面的为0。
B(x)A1(x)0(modxn2)

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

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

(B(x)A1(x))20(modxn)

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

接着把平方拆开。

B2(x)+A2(x)2B(x)A1(x)0(modxn)

A2(x)会很不好处理,我们整体乘一个A(x)。(显然A(x)非0,所以不会有问题。)
得到:
B2(x)A(x)+A1(x)2B(x)0(modxn)

这样我们就凑出了现在的逆元和之前的逆元的关系了。
A1(x)B(x)(2B(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)D1(x)1(modxn)
2、没有模,不需要模除法的结果长度也是确定的。
我们知道逆元怎么做,那么我们试着用逆元来做除法。

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

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

AR(x)=xn1A(1x)=i=0nani1xi

也就是,我们把 1x 带入 A(x),然后乘上xn1
这个操作的意义是什么呢?
反转系数。
原来的某一项aixi,变为aixi×xn1=aixni1
也就是系数反转了。

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

AR(x)=xn1A(1x)=xm1B(1x)xnmD(1x)+xnm+1xm2R(1x)

也就是:
AR(x)=BR(x)DR(x)+xnm+1RR(x)

这时再模xnm+1
AR(x)BR(x)DR(x)(modxnm+1)

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

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

注意点

1、除法没有模,不要写多了。
2、要保证mn,否则没有意义。
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)(xx0)+f′′(x0)(xx0)22+

为了简单,只保留线性部分。
f(x)=f(x0)+f(x0)(xx0)=0

然后得到
x=x0f(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=ai1i

记得将空出来的一项设为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)(1lnf0(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

展开阅读全文

没有更多推荐了,返回首页