多项式常用操作归纳


以下的难度顺序应该是递增的...
建议从前往后进行学习...
参考资料:

多项式求逆,多项式取模,多项式开方 学习笔记-by KsCla
多项式的多点求值与快速插值-by Miskcoo


目录
  • 多项式求逆
  • 多项式取模
  • 多项式的多点求值
  • 多项式的多点插值
      • nlog(n)^3的做法
      • nlog(n)^2的做法

多项式求逆

其实就是改变了模运算的形式,然后用到多项式上来。大概是这个形式:

\[A(x)*B(x) \equiv 1 (mod\ x^n) \]

其中A(x)和B(x)都是多项式,A(x)是输入的多项式,B(x)是我们需要求的在模x^n意义下的多项式。

可以看出B(x)的最高次项的次数一定是小于n的。

考虑如何进行求解。
设当前需要求解的是下边这个式子的B(x):

\[A(x)*B(x)\equiv 1(mod\ x^n)-----(1) \]

那么我们考虑已经求得下边这个式子的B'(x):

\[A(x)*B'(x)\equiv 1(mod\ x^{\lceil\frac{n}{2}\rceil})-----(2) \]

(以下请自行脑补上取整)
那么我们将第一个式子变化为:

\[A(x)*B(x)\equiv 1(mod\ x^{\frac{n}{2}})-----(3) \]

然后将(3)-(2)有:

\[A(x)*(B(x)-B'(x))\equiv 0(mod\ x^{\frac{n}{2}})-----(4) \]

最后将两边同时乘上B(x)(式子(1)):

\[B(x)-B'(x)\equiv 0(mod\ x^{\frac{n}{2}})-----(5) \]

这个式子就比较好了。因为要回到mod x^n的意义上来,我们考虑将这个式子进行平方。
令新的多项式\(F(x)=B(x)-B'(x)\),那么就有:\(F(x)\equiv 0 (mod\ x^{\frac{n}{2}})\),也就是说F(x)中指数

  1. 通过两个指数
  2. 通过一个指数=n/2的项相乘得到,然后系数仍然为0;
  3. 不可能通过两个指数均>=n/2的项得到,否则指数加起来就>=n了,直接就被x^n模掉了。

于是就可以得到:

\[B(x)*B(x)-2*B(x)*B'(x)+B'(x)^2\equiv 0(mod\ x^n) \]

再利用之前的套路,让两边同时乘上A(x)(式子(1)),就有:

\[B(x)-2*B'(x)+A(x)*B'(x)^2\equiv 0(mod\ x^n) \]

移项后就有:

\[B(x)\equiv 2*B'(x)-A(x)*B'(x)^2(mod\ x^n) \]

然后就可以通过不断地分治,求得子问题后,进行NTT乘法,求得B(x)。
总的时间复杂度是\(O(n \log n)\),然后还要乘上一个巨大的常数...

有了多项式求逆之后我们就有了一个比较有利的工具了。
模板的话,可以看下面多项式取模的模板,里面有多项式求逆。

多项式取模

Update:2019.5.23
终于有时间来写一下这个的推导了。
多项式取模就是多项式除法。大致形式为:\(A(x)=B(x)D(x)+R(x)\) ,其中\(A(x),B(x)\)都是已知的,分别可以看做是被除数和除数。然后\(D(x)\)可以看做是商,\(R(x)\)可以看做是余数,那么\(D(x)\)\(R(x)\)就是我们要求的。大致思路是通过倒置先求出\(D(x)\),然后再求出\(R(x)\)

假设\(A(x)\)的最高次数为\(n-1\)\(B(x)\)的最高次数为\(m-1\)。那么就可以推出\(D(x)\)的最高次数为\(n-m\),而\(R(x)\)的最高次数为\(m-2\)

直接上做法了:

首先将原始进行一个等价的变形,也就是将\(x\)替换为\(\frac{1}{x}\):

\[A(\frac{1}{x})=B(\frac{1}{x})D(\frac{1}{x})+R(\frac{1}{x}) \]

然后左右两边同时乘上\(x^{n-1}\),乘到每个多项式里面去:

\[A'(x)=B'(x)D'(x)+x^{n-m+1}R'(x) \]

仔细观察后,发现,其实\(A'(x),B'(x),D'(x),R'(x)\)分别是原函数倒过来的形式,即:

\[A(x)=\sum_{i=0}^{n-1}a_ix^i,A'(x)=\sum_{i=0}^{n-1}a_{n-i-1}x^i \]

那么再对左右两边同时模\(x^{n-m+1}\),就有:

\[A'(x)\equiv B'(x)D'(x) (mod\ x^{n-m+1}) \]

然后又因为\(D'(x)\)本来最高次数就只有\(n-m\),而模的是\(x^{n-m+1}\),那么就刚好保留了完整的多项式。也就是说我们只要能解上面的方程,就能够得到\(D(x)\)了。

然后就有:

\[A'(x) B'^{-1}(x)\equiv D'(x)(mod\ x^{n-m+1}) \]

也就是通过之前的多项式取逆的操作,把\(B'(x)\)移过去,就能够求出\(D'(x)\)了。
之后再通过将\(D'(x)\)变回\(D(x)\),就能够通过\(A(x)=B(x)D(x)+R(x)\)求出\(R(x)\)了。问题解决。
模板题:P4512 【模板】多项式除法
代码:

#include
#include
#include
#define MAXN 500000
#define MO 998244353
#define G 3
using namespace std;
int A[MAXN+5],B[MAXN+5],D[MAXN+5],R[MAXN+5];
int PowMod(int a,int b)
{
	int ret=1;
	while(b)
	{
		if(b&1)
			ret=1LL*ret*a%MO;
		a=1LL*a*a%MO;
		b>>=1;
	}
	return ret;
}
void NTT(int P[],int len,int opt)
{
	for(int i=1,j=0;i>=1,~j&s;);
		if(i

多项式的多点求值

在有了以上的两个工具之后,我们就可以进行更加复杂的操作了。
(剩下两部分的讲解可能有些偏差,如果感觉有和自己想得不一样的地方,可以参考之前给出的链接)

进入正题。

首先,多项式的多点求值主要是解决这样一个问题:给定一个多项式A(x)和n个点\(x_0,x_1,x_2,x_3,x_4...x_n\),然后要你求出\(A(x_0),A(x_1),A(x_2)...A(x_n)\)
不难想到有一个n^2的算法,就是直接对于每一个\(x_i\)枚举每一项加和得到\(A(x_i)\)
但是这显然是不够的。
我们考虑这样一个方法:
假设当前的问题是有一个x[]集合为\(x_0,x_1...x_n\),然后多项式为A(x)。
将我们需要带进去的x[]集合划分为两个部分:

\[LX[]=x_0,x_1,...x_{\lfloor\frac{n}{2}\rfloor} \]

\[RX[]=x_{\lfloor\frac{n}{2}\rfloor+1},...,x_n \]

然后再构造这样两个多项式LP(x),RP(x):

\[LP(x)=\prod_{i=1}^{\lfloor\frac{n}{2}\rfloor}(x-x_i) \]

\[RP(x)=\prod_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1}(x-x_i) \]

有了这两个多项式之后,我们再让A(x)分别模LP(x)和RP(x)(可以类比实数的运算),这里就只讨论LP(x),RP(x)的部分的处理是一模一样的。
那么就有:

\[A(x)=D(x)*LP(x)+A'(x) \]

\(x\in LX[]\)时,LP(x)由它的定义可得是等于0的。那么就直接得到了A'(x)。
再把它写成模运算的形式:

\[A(x)\equiv A'(x)(mod\ LP(x)) \]

然后就利用多项式的模运算求解就可以了。然后再讲A'(x)传给下一层,在底层(n=0)的时候就是一个常数了,就把这个位置的点值求出来了。

其实这里在参考博客中有一个问题似乎没有说清楚,就是怎么求LP(x)。虽然有定义,但是那个是无法利用的。根据我的理解,应该是先进行一次预处理,对于每一层需要用到的LP(x)和RP(x)都先用分治FFT求出来:就是在底层的时候直接返回,然后在将两个区间进行合并的时候,就对两个区间的LP(x)进行卷积(乘法),然后就求出所有可能用到的LP(x)和RP(x)了。

然后预处理的复杂度大概是\(O(n\log^2 n)\)的,求答案的复杂度大概也是\(O(n\log^2 n)\)的,但是同样也是常数巨大,无论怎么优化,都很容易T。所以小心为上。
贴代码咯~~~
版题:P5050 【模板】多项式多点求值
这道题硬是写了我一个晚上...4个小时啊...
因为担心vector会T,所以说我这里采用的是开一个巨大无比的数组,再记录一个int*指针,然后每一次都从这个大数组中静态地分配内存(常用操作啦),然后你就需要特别注意当前在大数组中的位置,注意不要超出应有的位置。最好的处理方式是每一个函数都开一些static int的数组,然后每次都先在这些数组上备份之后再搞,就不用担心越界了。
再就是求多项式的Mod的时候,需要特别注意static的数组需要进行清零,然后清零的上界也是需要进行特别注意的。否则会T。

#include
#include
#include
#define MAXN 640000
#define MAXLEN 3000000
#define MO 998244353
#define G 3
using namespace std;
int seq[MAXLEN*6+5];
int *lastp=&seq[0];
int n,m,x[MAXN*2+5],a[MAXN*2+5];
int ans[MAXN*2+5];
struct node
{
    int len,*p,*p2;
}tree[MAXN*10];
int *NewNode(int len)
{
    lastp+=len;
    return lastp-len;
}
int PowMod(int X,int Y)
{
    int ret=1;
    while(Y)
    {
        if(Y&1)
            ret=1LL*ret*X%MO;
        X=1LL*X*X%MO;
        Y>>=1;
    }
    return ret;
}
void Reverse(int X[],int len)
{
    for(int i=0;i>=1,~j&s;);
        if(i=deg)
            tmpA[i]=0;
        else
            tmpA[i]=A[i];
    for(int i=0;i=deg)
    		tmpB[i]=0;
    	else
    		tmpB[i]=B[i];
    NTT(tmpA,p,1);NTT(tmpB,p,1);
    for(int i=0;i

多项式的多点插值

多点插值其实是这样的一个问题:
给定一个点的集合:\(Point=(x_0,y_0),(x_1,y_1),(x_2,y_2)...(x_n,y_n)\),然后要求你求出一个多项式A(x),使得\(y_0=A(x_0),y_1=A(x_1),y_2=A(x_2)...y_n=A(x_n)\)

目前我看见的有两种做法。一种是\(O(nlog_n^3)\)的,还有一种是\(O(nlog_n^2)\)的。

nlog(n)^3的做法

其实有一点像多项式的多点求值的方法在进行讨论,但千万注意概念不要混淆了,特别是在看参考博客的时候。
还是先假设当前的点的集合为{\((x_0,y_0),(x_1,y_1),(x_2,y_2)...(x_n,y_n)\)},然后我们需要求的多项式为A(x)。
首先将点的集合分为两个部分:

\[LX=(x_0,y_0),(x_1,y_1),...,(x_{\lfloor\frac{n}{2}\rfloor},y_{\lfloor\frac{n}{2}\rfloor}) \]

\[RX=(x_{\lfloor\frac{n}{2}\rfloor+1},y_{\lfloor\frac{n}{2}\rfloor+1}),...,(x_{n-1},y_{n-1 }) \]

然后我们在构造一个“似曾相识”的多项式:

\[P(x)=\prod_{i=1}^{\lfloor\frac{n}{2}\rfloor}(x-x_i) \]

然后再让A(x)去模P(x)(果然和上面的很像):

\[A(x)=A_1(x)P(x)+A_2(x) \]

但是其中的A(x),A1(x),A2(x)都是不知道的。
先不管,将LX集合递归下去求A2(x)是不需要任何的前提条件的,那么我们就先递归把A2(x)求出来吧。
回溯回来之后,我们发现还有另外一半的RX点集没有使用。上面我们设出来的式子,对于LX和RX应该都是成立的,所以说就有以下的等式:

\[对于(x_i,y_i)\in RX,y_i=A_1(x_i)P(x_i)+A_2(x_i) \]

\[->A_1(x_i)=\frac{y_i-A_2(x_i)}{P(x_i)} \]

那么就利用多点求值快速求出A2(xi)和P(xi)了,然后可以求出A1(x)对应的n/2个点值,进而继续递归进行求解了。
这样子处理完了之后,A1(x),A2(x)就都被求出来了,那么再根据上面的式子,就能够直接把A(x)给求出来了。
但是先不慌,看一波时间复杂度:还需要利用多点求值!!那么就是\(O(log_n*nlog^2n)=O(nlog^3n)\),再乘上一大坨常数,这个方法几乎与O(n^2)无异了...

nlog(n)^2的做法

首先有一个结论,也就是拉格朗日插值:

\[F(x)=\sum_{i=1}^{n}\frac{\prod_{j!=i}(x-x_j)}{\prod_{j!=i}(x_i-x_j)}y_i \]

这个式子看起来非常的复杂,但其实理解之后还好(但其实还是推不出来...)。可以这样感性认知一下:

假设对于带进去的是\(x_i\),那么当k!=i的时候,分母的上面总会枚举到j使得ji,然后就是0了,根本不用看分母了;然后ki的时候,且j!=i(k),然后你就会发现分子和分母是同一个式子,那么自然就等于1了,再乘上后面的\(y_i\),自然就是\(y_1\)了。

这样子就能够理解上面的式子了。

现在考虑求分母。
设多项式\(M(x)=\prod_{i=1}^{n}(x-x_i)\),然后很容易就能够看出拉格朗日插值的式子中的第i项的分子就是\(\frac{M(x)}{x-x_i}\)。但是你会发现用M(x)来表示原式中的分母的时候,就是\(\frac{M(x_i)}{(x_i-x_i)}\),而这个式子的分子与分母都是0!!!然而利用我们学过的高中知识,也就是洛必达法则,然后就可以得到其实这就是M'(x)。这个时候,聪明的读者就会发现M'(x)在全局都是一样的,就可以先求出M'(x),然后利用多点求值快速求出答案了。

然后我们设每一项的系数就是\(val_i=\frac{y_i}{\prod_{}(x_i-x_j)}\),正式开始考虑如何求这个多项式了。再设\(LP(x)=\prod_{i=1}^{\lfloor\frac{n}{2}\rfloor}(x-x_i)\)\(RP(x)=\prod_{i=\lfloor\frac{n}{2}\rfloor+1}^{n}(x-x_i)\) ,然后原多项式(即F(x))就可以表示为:

\[F(x)=\sum_{i=1}^{\lfloor\frac{n}{2}\rfloor}(val_i*\prod_{j!=i}^{j<=\lfloor\frac{n}{2}\rfloor}(x-x_j)RP(x))+\sum_{i=\lfloor\frac{n}{2}\rfloor+1}^{n}(val_i*\prod_{j!=i}^{j>\lfloor\frac{n}{2}\rfloor}(x-x_j)LP(x)) \]

然后我们可以发现左右两边几乎形成了一个可以递归的式子(左边算的就是原式中的<=n/2的部分,右边算的就是>n/2的部分)。特别需要注意的是,需要将其与拉格朗日插值区分开来,且需要将val[i]也化到递归的下一层去,那么递归的底层就是l==r的时候,这个时候应该返回val[l]。
然后分析一波时间复杂度。大概是\(O(nlogn*logn)\)的。但也和众多的多项式运用一样,有一大坨常数...

贴代码贴代码~~~
\(O(nlog^2n)\)的做法:
代码极丑,效率中等...
其中Mod的部分没有进行批注,可以在多点插值那里去看详细的批注。
版题:P5158 【模板】多项式快速插值

#include
#include
#include
#define MAXN 1000000
#define MAXLEN 2000000
#define MO 998244353
#define G 3
using namespace std;
int n;
int x[MAXN+5],y[MAXN+5];
int M[MAXN+5];
int seq[MAXLEN*9+5];
int *ncnt=&seq[0];
int v[MAXN+5],w[2][MAXN*2+5];
struct node
{
    int *p,*p2,*p3,len;
}tree[MAXN*4];
inline int *NewNode(int len)
{
    ncnt+=len;
    return ncnt-len;
}
int PowMod(int a,int b)
{
    int ret=1;
    while(b)
    {
        if(b&1)
            ret=1LL*ret*a%MO;
        a=1LL*a*a%MO;
        b>>=1;
    }
    return ret;
}
void Prepare()
{
	int S=1;
	while(S>1);s>0;s>>=1)
	{
		w[0][s]=1LL*w[0][s<<1]*w[0][s<<1]%MO;
		w[1][s]=1LL*w[1][s<<1]*w[1][s<<1]%MO;
	}
}
void Reverse(int X[],int len)
{
    for(int i=0;i>=1,~j&s;);
        if(i'9')
    {
        c=getchar();
        if(c=='-')	f*=-1;
    }
    ret=10*ret+c-'0';
    while(true)
    {
        c=getchar();
        if(c<'0'||c>'9')	break;
        ret=10*ret+c-'0';
    }
    return ret*f;
}
void print(int x)
{
    if(x==0)
        return;
    print(x/10);
    putchar(x%10+'0');
}
int main()
{
	Prepare();//预处理NTT需要用到的"单位负根",否则会T
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
       	x[i]=read(),y[i]=read();
    Build(1,1,n);//先预处理M(x)组,同时也处理了LP(x)和RP(x)
    for(int i=0;i

各种实现的话,目前由于博主效率有限,后面会慢慢补上来的...
错误修正:2019.3.26
内容补充:2019.5.23