有关多项式处理的各种算法总结

有关多项式处理的各种算法总结

---

前言

最近频繁接触了各种各样的生成函数，以及各种十(sang)分(xin)有(bing)趣(kuang)的多项式相关算法。
由于这部分内容太过于丧病外加内容较多，因此决定写篇总结来辅助复习。

阅读本文所需前置技能：快速数论变换

以下内容中代码片涉及的多项式卷积(NTT)及各种函数使用代码如下:

typedef long long ll;
const ll md=998244353;
inline int add(int a,int b){a+=b;if(a>=md)return a-md;return a;}
inline int mul(int a,int b){return (ll)a*b%md;}
inline void initrev(int n)
{
    for(int i=0;i<n;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
}
inline void NTT(int *a,int n,bool f)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int h=2;h<=n;h<<=1)
    {
        int w=qpow(3,(md-1)/h);
        if(f)w=qpow(w,md-2);
        for(int j=0;j<n;j+=h)
        {
            int wn=1;
            for(int k=j;k<j+(h>>1);k++)
            {
                int x=a[k],y=mul(wn,a[k+(h>>1)]);
                a[k]=add(x,y);
                a[k+(h>>1)]=add(x-y,md);
                wn=mul(wn,w);
            }
        }
    }
    if(f)
        for(int i=0,invn=inv[n];i<n;i++)
            a[i]=mul(a[i],invn);
}

inline void mult(int *a,int n,int *b,int m,int *c)
{
    static int x[N],y[N],l;
    for(l=1;l<=(n+m-1);l<<=1);
    initrev(l);
    memcpy(x,a,sizeof(a[0])*n);
    memcpy(y,b,sizeof(b[0])*m);
    fill(x+n,x+l,0);fill(y+m,y+l,0);
    NTT(x,l,0);NTT(y,l,0);
    for(int i=0;i<l;i++)
        c[i]=mul(x[i],y[i]);
    NTT(c,l,1);
}

同时，绝大多数代码中传入的参数$n$为多项式长度，要求为一个不小于所需答案长度的一个$2$的次幂。

---

1.一些定义的简要介绍

非零多项式$G(x)$的次数为其最高项的次数，通常记为$deg(G(x))$。本文本着尽量直观的原则将避免使用此符号。

形如$[x^n]G(x)$的符号代表多项式$G(x)$的$x^n$一项的系数。

多项式之间的乘积叫做卷积，可以理解为初中学习的因式分解后展开的方法。
举个栗子就明白了:

$$(1+2x+3x^2)*(9+2x)=9+20x+31x^2+6x^3$$

多项式通常以生成函数的形式出现。这里简单介绍一下生成函数。

所谓生成函数，就是把一个数列用一个函数的各项系数来表示……
比如，这里以著名的斐波那契数列做例子:
假设$fib(1)=1,fib(2)=1,fib(n)=fib(n-1)+fib(n-2)$
得到数列:

$$Fib=\{1,1,2,3,5,8,13\cdots\}$$
令其生成函数为 $F(x)$，得到的函数如下:
$$F(x)=1+x^1+2x^2+3x^3+5x^4+8x^5+13x^6+\cdots$$
相信大家这就大致能明白生成函数是个什么意思了~
不如直接去看百科

对于一个无穷长的数列的生成函数，理论上会是个无穷长的多项式，而因此难以表示出来。
因此定义$mod\ x^n$意义下的多项式为仅保留前$n$项的的原多项式。

以上关于生成函数的内容和剩下的内容貌似没有什么太大关系。

---

2.牛顿迭代

问题:

给出多项式$G(x)$，求解多项式 $F(x)$，满足：

$$G(F(x)) \equiv 0\ (mod\ x^n)$$

需要给出$F(x)\ mod\ x^n$意义下的答案，即前$n$个系数。

---

解决方法(没看懂或不想看可直接到下面记忆结论，下面的内容同理):
假设咱们现在有

$$G(F_0(x)) \equiv 0\ (mod\ x^{2^t})$$
要求的是

$$G(F(x)) \equiv 0\ (mod\ x^{2^{t+1}})$$
也就是利用倍增的思想，将长度加倍。
注意如果求解的是特定长度的答案，可以考虑直接求出向上取最近的 $2$的次幂长度的答案，并保留所需的长度即可。

对$G(F(x))$在 $F_0(x)$处来一发泰勒展开:

$$G(F(x))=G(F_0(x))+\frac{G'(F_0(x))}{1!}(F(x)-F_0(x))^1+\frac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+\cdots$$
可以发现，在 $mod\ x^{2^{t+1}}$意义下， $F(x)-F_0(x)$的最低非 $0$次项大于$2^{t+1}$，于是从第三项开始，后面的项贡献均为 $0$。
那么得到一个简单的式子:
$$G(F(x))=G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\ (mod\ x^{2^{t+1}})$$
又由咱们的目标知 $G(F(x)) \equiv 0\ (mod\ x^{2^{t+1}})$
那么最终移项得到:
$$F(x) \equiv F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))}\ (mod\ {x^{2^{t+1}}})$$
这个形式正是咱们想要的~

---

3.多项式求逆

提出给定多项式$F(x)$，求多项式$G(x)$，使得下式成立:

$$G(x)*F(x)-1 \equiv 0\ (mod\ x^n)$$

---

解决方法:
考虑使用倍增，边界为$G_0=Inv(F_0)$，即直接求逆元。
考虑使用牛顿迭代公式，令$G(x)$为变量，$G_t(x)$代表$mod\ x^{2^{t}}$意义下的$G(x)$。
假设已经有了$G_t(x)$，要求$G_{t+1}(x)$，将$G_t(x)*F(x)-1$带入牛顿迭代公式，得：

$$\begin{align} G_{t+1}(x) & \equiv G_t(x) - \frac{G_t(x)*F(x)-1}{F(x)}\ (mod\ {x^{2^{t+1}}}) \\ & \equiv G_t(x) - (G_t(x)*F(x)-1)*G_t(x)\ (mod\ {x^{2^{t+1}}}) \\ & \equiv 2*G_t(x)-F(x)*G_t(x)*G_t(x)\ (mod\ {x^{2^{t+1}}}) \\ \end{align}$$
最后的式子便是结果！

代码实现(个人习惯较为明显，仅供参考):

inline void cinv(int *a,int *b,int n)
{
    static int C[N];
    if(n==0)
    {
        b[0]=qpow(a[0],md-2);
        return;
    }

    cinv(a,b,n>>1);
    memcpy(C,a,sizeof(a[0])*n);
    memset(C+n,0,sizeof(a[0])*n);
    n<<=1;initrev(n);
    NTT(C,n,0);NTT(b,n,0);
    for(int i=0;i<n;i++)
        b[i]=add(mul(2,b[i])-mul(C[i],mul(b[i],b[i])),md);
    NTT(b,n,1);n>>=1;
    memset(b+n,0,sizeof(b[0])*n);
}

---

4.多项式求$\ln$

问题:
给出多项式$F(x)$，求多项式$G(x)$，满足

$$\ln(F(x))-G(x) \equiv 0 \ (mod\ x^n)$$

---

解决方法:
考虑对两边同时求导:

$$G'(x) \equiv \frac{F'(x)}{F(x)}$$
那么有
$$G(x) \equiv \int \frac{F'(x)}{F(x)} dx$$
考虑多项式如何求导……
每一项同时左移一位并乘上原来下标即可~
积分同理，右移一位并把下标乘回去即可~
配合多项式求逆食用，复杂度为 $O(nlogn)$~

代码实现:

inline void cln(int *f,int *g,int n)
{
    static int D[N],A[N],m;
    memset(D,0,sizeof(D));
    memset(A,0,sizeof(A));

    for(int i=0;i<n;i++)
        D[i]=mul(f[i+1],(i+1));
    D[n]=0;

    for(m=1;m<=n;m<<=1);
    cinv(f,A,n);initrev(m);
    NTT(A,m,0);NTT(D,m,0);
    for(int i=0;i<m;i++)
        A[i]=mul(A[i],D[i]);
    NTT(A,m,1);

    for(int i=1;i<=n;i++)
        g[i]=mul(A[i-1],inv[i]);
    g[0]=0;
}

---

5.多项式开方

问题:
给出多项式$F(x)$，求多项式$G(x)$，满足

$$G(x)^2-F(x) \equiv 0\ (mod\ x^n)$$

---

解决方法:
考虑和多项式求逆一样使用倍增，初值$G_0=\sqrt{F(0)}\ (mod\ x)$。
同样将$G^2(x)-F(x)$带入牛顿迭代公式:

$$\begin{align} G_{t+1}(x) & = G_t(x)-\frac{G_t(x)^2-F(x)}{2G_t(x)}\ (mod\ {x^{2^{t+1}}}) \\ & = \frac{2G_t(x)^2-G_t(x)^2+F(x)}{2G_t(x)}\ (mod\ {x^{2^{t+1}}}) \\ & = \frac{G_t(x)^2+F(x)}{2G_t(x)}\ (mod\ {x^{2^{t+1}}}) \\ \end{align}$$

配合多项式求逆食用，分析可知复杂度为$O(nlogn)$

代码实现:

void csqrt(int *a,int *b,int n)
{
    static int C[N],D[N];
    if(n==1)
    {
        b[0]=sqrt(a[0]);
        return;
    }

    csqrt(a,b,n>>1);
    memset(D,0,sizeof(D));

    cinv(b,D,n);
    memcpy(C,a,sizeof(a[0])*n);
    memset(C+n,0,sizeof(C[0])*n);
    n<<=1;NTT(C,n,0);NTT(b,n,0);NTT(D,n,0);
    for(int i=0;i<n;i++)
        b[i]=mul(inv2,add(b[i],mul(C[i],D[i])));
    NTT(b,n,1);n>>=1;
    memset(b+n,0,sizeof(b[0])*n);
}

---

6.多项式求指数函数(exp)

问题:
给出多项式$F(x)$，求多项式$G(x)$，满足

$$e^{F(x)}-G(x) \equiv 0\ (mod\ x^n)$$

---

解决方法:
依旧考虑使用倍增和牛顿迭代。
然而发现一个问题:对于指数函数，无法知道常数项。
那么考虑通过同时取$\ln$来变形:

$$F(x)-\ln(G(x)) \equiv 0\ (mod\ x^n)$$
然后把这个式子带入牛顿迭代式:
$$\begin{align} G_{t+1}(x) & = G_t(x)-\frac{F(x)-\ln(G_t(x))}{-\frac{1}{G_t(x)}}\ (mod\ {x^{2^{t+1}}}) \\ & = G_t(x)+G_t(x)(F(x)-\ln(G_t(x)))\ (mod\ {x^{2^{t+1}}}) \\ \end{align}$$
配合多项式求逆和多项式求 $\ln$共同食用，复杂度仍为 $O(nlogn)$。

代码实现:

inline void cexp(int *a,int *b,int n)
{
    static int D[N];
    if(n==1)
    {
        b[0]=1;
        return;
    }

    cexp(a,b,n>>1);
    memset(D,0,sizeof(D));
    cln(b,D,n);
    for(int i=0;i<n;i++)
        D[i]=add(a[i]-D[i],md);
    D[0]=add(1,D[0]);
    n<<=1;initrev(n);
    NTT(D,n,0);NTT(b,n,0);
    for(int i=0;i<n;i++)
        b[i]=mul(D[i],b[i]);
    NTT(b,n,1);n>>=1;
    memset(b+n,0,sizeof(b[0])*n);
}

---

7.多项式快速幂

问题:
给出多项式$F(x)$和正整数$k$，求多项式$G(x)$，满足

$$F(x)^k-G(x) \equiv 0\ (mod\ x^n)$$

---

解决方法:
对等式求$\ln$得:

$$k\ln(F(x))-\ln(G(x)) \equiv 0\ (mod\ x^n)$$
因此
$$G(x)=exp(k\ln(F(x)))$$
配合多项式求 $\ln$和多项式求 $exp$共同食用，复杂度依旧为 $O(nlogn)$。

代码实现:

inline void cpow(int *f,int *g,int n,int k)
{
    static int C[N];
    memset(C,0,sizeof(C));

    cln(f,C,m);
    for(int i=0;i<n;i++)
        C[i]=mul(C[i],k);
    cexp(C,g,m);
}

---

8.拉格朗日反演

问题:
给出多项式$F(x)$，求多项式$G(x)$的$n$次项系数，满足

$$F(G(x)) = x$$

---

解决方法:
采用拉格朗日反演:

$$[x^n]G(x)=\frac{1}{n}[x^{n-1}]\frac{1}{F(x)^n}$$
本蒟蒻不会证明0_0
但是看完了上面的内容后显然还是会用的~
配合多项式求逆和多项式快速幂食用即可。

代码实现:

inline int lagrange(int *f,int n,int x)
{
    static int C[N],D[N];
    memset(C,0,sizeof(C));
    memset(D,0,sizeof(D));
    cinv(f,C,n);
    cpow(C,D,n,x);
    return mul(D[x-1],qpow(x,md-2));
}

---

9.多项式除法和取模

问题:
给出$n$次多项式$F(x)$和$m$次多项式$D(x)$，求多项式$G(x)$和$P(x)$，满足

$$F(x) = D(x)*G(x)+P(x)$$
且
$$deg(P(x))<m$$

---

解决方法:
定义运算$A^R(x)$，满足对于$deg(A(x))=n$的多项式$A(x)$：

$$A^R(x)=x^nA(\frac{1}{x})$$
或者说
$$[x^i]A^R(x)=[x^{n-i-1}]A(x)$$
再举个例子
$$A(x)=1+2x+3x^2+4x^3$$
$$A^R(x)=4+3x+2x^2+1$$
也就是系数反转的意思~

那么现在来考虑将原多项式反转一下:

$$x^nF(\frac{1}{x})=x^mD(\frac{1}{x})x^{n-m}G(\frac{1}{x})+x^{n-m+1}x^{m-1}R(\frac{1}{x})$$
注意此处默认 $deg(G(x))=n-m$， $deg(R(x))=m-1$，高次项补 $0$。

此时，$x^{n-m+1}x^{m-1}R(\frac{1}{x})$最低次项为 $x^{n-m+1}x^{m-1}(\frac{1}{x})^{m-1}=x^{n-m+1}$，而 $deg(G(x))=n-m$。
则放到 $mod\ n-m+1$意义下， $R(x)$的影响被完全抵消，而不会影响 $G(x)$的结果，更不会影响已知的 $F(x)$和 $D(x)$。

那么有

$$F^R(x)=D^R(x)G^R(x)\ (mod\ x^{n-m+1})$$

于是来一发多项式求逆+reverse函数即可求出$G(x)$，进而求出$R(x)$。复杂度$O(nlogn)$。

代码实现:

    inline void cdiv(int *a,int n,int *b,int m,int *ret)
    {
        static int c[N],d[N],e[N],l;
        if(n<m)return;
        for(l=1;l<=n-m+1;l<<=1);
        reverse_copy(a,a+n,c);
        reverse_copy(b,b+m,d);
        for(int i=m;i<n+m+1;i++)//不能用fill(),否则RE!!!
            d[i]=0;
        memset(e,0,sizeof(e[0])*(l*2));
        cinv(d,e,l);
        mult(e,n-m+1,c,n-m+1,d);
        reverse_copy(d,d+n-m+1,ret);
    }

    inline void cmod(int *a,int n,int *b,int m,int *ret)
    {
        static int c[N];
        cdiv(a,n,b,m,c);
        mult(c,n-m+1,b,m,c);
        for(int i=0;i<m-1;i++)
            ret[i]=add(a[i]-c[i],md);
    }

---

10.多项式的多点求值

问题:
给出多项式$F(x)$和点集$X=\{x_0,x_1,\cdots,x_{n-1}\}$，求

$$Y=\{F(x_0),F(x_1),\cdots,F(x_{n-1})\}$$

---

解决办法:
考虑分治。
如何把问题变成两个规模减半的子问题?
首先把所求的值划成两半:

$$\begin{align} X_0 & = \{x_0,x_1,\cdots,x_{\lfloor \frac{n}{2} \rfloor}\} \\ X_1 & = \{x_{\lfloor \frac{n}{2} \rfloor+1},x_{\lfloor \frac{n}{2} \rfloor+2},\cdots,x_{n-1}\} \end{align}$$
此时，答案集合分别为 $Y_0$和 $Y_1$。

然后系数也要划成两半。
考虑构造两个多项式:

$$\begin{align} G_0(x) & = \prod\limits_{i=0}^{\lfloor \frac{n}{2} \rfloor} (x-x_i)\\ G_1(x) & = \prod\limits_{i=\lfloor \frac{n}{2} \rfloor+1}^{n-1} (x-x_i)\\ \end{align}$$
可以发现，当带入的点 $x_i \in X_0$，那么 $G_1(x_i)=0$。同理当 $x_i \in X_1$，那么 $G_0(x_i)=0$。

令点集$X_0$在知道$Y_0$的情况下可以解出的多项式为$F_0(x)$，$X_1$在知道$Y_1$的情况下可以解出的多项式为$F_1(x)$。
可以发现，根据定义，$F_0(x)$和$F_1(x)$正是咱们需要的子问题系数，因为它们分别满足了带入$X_0$和$X_1$能解出正确的$Y_0$和$Y_1$。

现在可以用如下方式表示$F(x)$：

$$\begin{align} F(x) & = D_0(x)*G_0(x)+F_0(x)\\ F(x) & = D_1(x)*G_1(x)+F_1(x)\\ \end{align}$$
其中 $D_0(x)$和 $D_1(x)$为任意多项式。
可以发现在上式中，若带入的值 $x_i \in X_0$，则 $D_0(x)*G_0(x)=0$，此时 $F(x)=F_0(x)$。
同理，若带入的值 $x_i \in X_1$，则 $D_1(x)*G_1(x)=0$，此时 $F(x)=F_1(x)$。

也就是说，如下关系成立:

$$\begin{align} F(x) & \equiv F_0(x)\ (mod\ G_0(x))\\ F(x) & \equiv F_1(x)\ (mod\ G_1(x))\\ \end{align}$$

这样得到的子问题系数$F_0(x)$和$F_1(x)$满足次数不超过递归区间长度，同时可以求出正确的结果，因此这样做是正确的。
最后递归到底层时，余下来的系数便是答案!

配合多项式取模食用，复杂度$O(nlog^2n)$。
注意，使用前需要先使用分治预处理出每一层分治区间的$G(x)$，预处理复杂度同样为$O(nlog^2n)$。

代码实现:

inline void pre(int dep,int l,int r,int *a,int b[M][N])
{
    if(l==r-1)
    {
        b[dep][l<<1]=add(md,-a[i]),b[dep][l<<1|1]=1;
        return;
    }
    int mid=l+r>>1;
    pre(dep+1,l,mid,a,b);
    pre(dep+1,mid,r,a,b);
    mult(b[dep+1]+(l<<1),mid-l+1,b[dep+1]+(mid<<1),r-mid+1,b[dep]+(l<<1));
}

inline void calc(int dep,int l,int r,int *a,int b[M][N])
{
    static int c[N];
    if(l==r-1)return;int mid=l+r>>1;
    memcpy(c,a+l,sizeof(a[0])*(r-l));
    cmod(c,r-l,b[dep+1]+(l<<1),mid-l+1,a+l);
    cmod(c,r-l,b[dep+1]+(mid<<1),r-mid+1,a+mid);
    calc(dep+1,l,mid,a,b);
    calc(dep+1,mid,r,a,b);
}

---

不出意外还是未完待续…..
（因为还没学快速插值）

---

写到这里才惊讶地发现咱最近都学的是些什么毒瘤东西……
不过有用就是坠吼的~

一些比较好的练手题

BZOJ 4555 求和 （多项式求逆）
BZOJ 3456 城市规划 （多项式求$\ln$）
BZOJ 3625 小朋友和二叉树 （多项式开方）
BZOJ 3684 大朋友和多叉树 （拉格朗日反演）标程
COGS 2189 帕秋莉的超级多项式 （使用到的方法为第2-7条）标程