多项式的运算

多项式

目录

• 多项式的运算
  • 基本运算
    • 多项式加减
    • 多项式乘法
    • 多项式除法与取模
    • 多项式求导与积分
    • 多项式快速幂
  • 多项式的初等函数
    • 多项式求逆
    • 多项式ln
    • 多项式的牛顿迭代法
      • 多项式开根
      • 多项式exp
        • 多项式的三角函数
  • 多项式的除法与取模
  • 常数大到没法考的东西
    • 多项式多点求值
    • 多项式快速插值

这里主要写一下基于FFT的多项式运算，应用方面主要是生成函数，和特征多项式优化线性递推式，见那两篇博客。
感觉还是要先看一下生成函数那个博客，可以不会写代码，但是要知道生成函数是什么东西、怎么列式子，不然会不知道下面要求啥以及是干啥用的...

前置技能：[FFT优化的多项式乘法]("快速傅里叶变换")

基本运算

多项式加减

...

多项式乘法

都说了前置技能...

多项式除法与取模

用到了求逆，最后说。

多项式求导与积分

不会的话去问数学老师...

多项式快速幂

给定多项式$A(x)$，求$A^k(x) \mod{x^n}$，$k \in\mathbb{N}$。

• 直接类似正常快速幂做就好了，只是每次要暴力把$x^{n+1}$项及以上清零，复杂度$O(n\log{n}\log{k})$。
• 模$x^n$是为了能算，而且一般这个都是求生成函数用的，都说了先看那篇博客...
• 但是这个做法其实后来就不需要了，下面会说怎么做到一个$\log$，而且$k$也不仅仅局限于$\mathbb{N}$了，是多项式都可以...
• 不过对于乘法定义为：$(f \ast g)(x)=\sum_{(i+j)\mod{n}=x}f(i)\times g(j)$的循环卷积，还是只能用这个求，就是把清零的过程变为把下标模下去。

循环卷积的快速幂：

inline void mod_conv(int a[],int b[],int c[],int len){
    static int tmpa[MAXL],tmpb[MAXL],tmpc[MAXL];
    memcpy(tmpa,a,len<<2);
    memcpy(tmpb,b,len<<2);
    poly_mul(tmpa,tmpb,tmpc,len);
    memset(c,0,len<<2);
    for(int i=0;i<(len<<1);++i)
        c[i]=(c[i]+tmpc[i%len])%P;
}
inline void poly_pow(int x[],int y,int ret[],int len){
    static int tmpx[MAXL];
    memcpy(tmpx,x,len<<2);
    memset(ret,0,len<<2);
    ret[0]=1;
    for(;y;y>>=1){
        if(y&1) mod_conv(ret,tmpx,ret,len);
        mod_conv(tmpx,tmpx,tmpx,len);
    }
}

多项式的初等函数

已知一个多项式$A(x)$，我们要求的是函数$f(A(x))$，其中函数$f$是一个初等函数，但是自变量为多项式。这个操作直观意义不是很明显，主要是解决一些生成函数问题用的，所以可以理解为将多项式带入函数的泰勒展开式，譬如$e^{A(x)}$就是求：$\sum_{i=0}^{+\infty}\frac{A^i(x)}{i!}$，而并非真要算$e$的“多项式次方”。我们都知道生成函数是形式幂级数，在计算时只需要保留有用的前$n$项。而初等函数运算就是将多项式带入一个泰勒展开式，也只能得到无穷级数，而它主要是为生成函数运算服务的，所以我们也可以把它放在$\mod{x^{n+1}}$意义下求，这样就方便计算了。
注意：下文有时为了突出多项式是函数自变量，会将多项式$f(x)$写成$f$，所有大写的$F(x)$之类的则代表自变量为多项式的函数。

多项式求逆

也就是函数$F(f)=f^{-1}$，但这个显然无法快速幂，而且这个由于比较特殊，相当于乘法逆元（是后面基本上所有东西的基础），我们要单独想一种做法...

• 我们考虑从模低次意义下向上推，设当前的答案为$g$。
• 首先在$\mod{x}$意义下$g=([x^0]f)^{-1}$；
• 假设我们已经得出了$f$在$\mod{x^{\lceil\frac{n}{2}\rceil}}$意义下的逆元$g$，考虑怎么得到在$\mod{x^n}$意义下的逆元$h$：
  $$\begin{align*} fg &\equiv 1 \pmod{x^{\lceil\frac{n}{2}\rceil}} \\ fg-1 &\equiv 0 \pmod{x^{\lceil\frac{n}{2}\rceil}} \\ (fg-1)^2 &\equiv 0 \pmod{x^n} \\ f^2g^2-2fg+1 &\equiv 0 \pmod{x^n} \\ f^2g^2h-2fgh+h &\equiv 0 \pmod{x^n} \\ fg^2-2g+h &\equiv 0 \pmod{x^n} \\ h &\equiv g(2-fg) \pmod{x^n} \end{align*}$$
• 为什么那个平方操作会让$\mod{x^{\lceil\frac{n}{2}\rceil}}$变成$\mod{x^n}$？这是只有同余式为$0$的时候才成立的，因为左式平方后，多项式变为$n-1$次了，但是依然是$0$，所以改变模数不会影响。
• 我们就按上面的式子进行迭代就可以了。同时也可以发现，多项式存在逆元的充要条件是常数项有逆元。
• 复杂度是：$T(n)=T(\frac{n}{2})+O(n\log{n})$，也就是$O(n\log{n})$，不过这个复杂度分析看起来好假啊...

代码：

inline void poly_inv(int a[],int inva[],int len){
    static int tmpa[MAXL];
    memset(tmpa,0,len<<1<<2);
    memset(inva,0,len<<1<<2);
    inva[0]=inv(a[0]);
    for(int i=2;i<=len;i<<=1){
        memcpy(tmpa,a,i<<2);
        NTT(tmpa,i<<1,1),NTT(inva,i<<1,1);
        for(int j=0;j<(i<<1);++j)
            inva[j]=(2-(long long)tmpa[j]*inva[j]%P+P)*inva[j]%P;
        NTT(inva,i<<1,-1),div_len(inva,i<<1);
        memset(inva+i,0,i<<2);
    }
}

多项式ln

主要是配合后面的$\exp$来用的。
由于$\ln(x)$在$x=0$出无定义，我们不方便泰勒展开，所以用了一种同后面牛顿迭代不同的方法。

• 设$f=\ln{g}$：
  

  $$\begin{align*} f' &= \frac{g'}{g} \\ f &= \int{\frac{g'}{g}} \end{align*}$$

• 求导和积分都是$O(n)$的，求个逆再乘一下就好了，复杂度$O(n\log{n})$。

代码：

inline void poly_ln(int a[],int lna[],int len){
    static int da[MAXL],inva[MAXL],dlna[MAXL];
    memset(da,0,len<<1<<2);
    poly_diff(a,da,len);
    poly_inv(a,inva,len);
    poly_mul(da,inva,dlna,len<<1);
    poly_integ(dlna,lna,len);
}

多项式的牛顿迭代法

求函数$f(x)$的零点常用的一种方法是牛顿迭代法，其公式$x_{t+1}=x_t-\frac{f(x_t)}{f'(x_t)}$可以解释为每次拿$f(x)$在当前的近似解$x_t$处的切线与$f(x)$的交点作为新的解。这样不断迭代，通常能收敛到零点处。但是这个解释是有些片面的，真正的原理是每次在当前解$x_t$处将$f(x)$进行泰勒展开：$f(x)=\sum_{i=0}^{+\infty}\frac{f^{(i)}(x_t)}{i!}(x-x_t)^i$，带入我们要求的零点，同时省略后面的项，可以得到一个方程解出一个$x_{t+1}$。由于上面的省略操作，这个是近似的解。
对于多项式的一些运算，我们可以化为求一个自变量为多项式的函数$H(f)$的零点。考虑将牛顿迭代法推广过来：

• 设当前的解为$g$，类比多项式求逆的方法，我们从模$x$的低次项开始推。
• 在$\mod{x}$的意义下，显然$g=H([x^0]f)$；
• 假设我们得到了$\mod{x^{\lceil\frac{n}{2}\rceil}}$意义下的解$g_0$，考虑迭代出$\mod{x^n}$意义下的解$g$；
• 泰勒展开后我们会发现，对于$\sum$中$i>1$的项，其最低次数都不小于$n$，所以直接可以去掉，而求$g$的部分不会受到影响，所以多项式的牛顿迭代法求的是准确解！
  $$g \equiv g_0-\frac{H(g_0)}{H'(g_0)} \pmod{x^n}$$

多项式开根

• $H(g)=g^2-f$
• 其实多项式任意实数次方都可以像上面那样推出来，（求逆也可以，不过推出来的式子和之前的是一样的），但是太麻烦了。
• 其实模意义下常数项反而是最麻烦的，比如开方，弄不好还得求二次剩余，不过一般出题人都良心的让它是$1$。
• 复杂度分析和求逆一样，不过还套了个求逆，常数...

代码：

inline void poly_sqrt(int a[],int sqrta[],int len){
    static const int INV2=inv(2);
    static int tmpa[MAXL],invsqrta[MAXL];
    memset(tmpa,0,len<<1<<2);
    sqrta[0]=sqrt(a[0]); //通常a[0]=1
    for(int i=2;i<=len;i<<=1){
        memcpy(tmpa,a,i<<2);
        poly_inv(sqrta,invsqrta,i);
        NTT(tmpa,i<<1,1),NTT(sqrta,i<<1,1),NTT(invsqrta,i<<1,1);
        for(int j=0;j<(i<<1);++j)
            sqrta[j]=((long long)sqrta[j]*sqrta[j]%P+tmpa[j])%P*INV2%P*invsqrta[j]%P;
        NTT(sqrta,i<<1,-1),div_len(sqrta,i<<1);
        memset(sqrta+i,0,i<<2);
    }
}

多项式exp

• 就是求$\sum_{i=0}^{+\infty}\frac{A^i(x)}{i!}$，因为这个式子在生成函数中很有意义。
• $H(g)=\ln{g}-f$。但是如果是在模$P$意义下，那么$f$的常数项必须为$0$，不然由于$e$是超越数，根本没法计算$e^{[x^0]f}\mod{P}$，而且这个式子通常是算方案数一类的问题，基本都是模$P$意义下了。
• 然后就有一个黑科技了：$A^k(x)=e^{k\ln{A}}$，$k$可以是任意一个多项式，这个可以在$O(n\log{n})$内求出了，虽然常数大，但还是优于快速幂的。

代码：

inline void poly_exp(int a[],int expa[],int len){
    static int tmpa[MAXL],lnexpa[MAXL];
    memset(tmpa,0,len<<1<<2);
    memset(expa,0,len<1<<2);
    assert(a[0]==0),expa[0]=1;
    for(int i=2;i<=len;i<<=1){
        memcpy(tmpa,a,i<<2);
        poly_ln(expa,lnexpa,i);
        NTT(tmpa,i<<1,1),NTT(expa,i<<1,1),NTT(lnexpa,i<<1,1);
        for(int j=0;j<(i<<1);++j)
            expa[j]=(long long)expa[j]*(tmpa[j]+1-lnexpa[j]+P)%P;
        NTT(expa,i<<1,-1),div_len(expa,i<<1);
        memset(expa+i,0,i<<2);
    }
}

多项式的幂：

inline void poly_pow(int a[],int k[],int powa[],int len){
    static int lna[MAXL],klna[MAXL];
    poly_ln(a,lna,len);
    poly_mul(lna,k,klna,len);
    poly_exp(klna,powa,len);
}

多项式的三角函数

• 求$\sin(f)$或者$cos(f)$，其它的都可以由它们变形出来。不过好像没什么应用价值...
• 考虑到三角函数求导还是三角函数，我们要想办法把它和之前可以求的函数联系起来，比如欧拉公式：$e^{ui}=\cos(u)+i\sin(u)$；
• 再结合$\cos(-x)=\cos(x),\sin(-x)=-\sin(x)$，我们就可以加减出来了；
• 复数意义下直接FFT算，反正FFT是支持多项式系数为复数的；
• 至于说模意义下，我本来以为可以有$i \equiv \sqrt{P-1} \pmod{P}$，求个二次剩余就好了。但是看picks大佬的博客说要用二元组模拟之类的，可能是因为复数是一个二维平面上的东西，不能简单压到一维去吧，暂时没有仔细研究。

多项式的除法与取模

给定$n$次多项式$A(x)$和$m$次多项式$B(x)$（通常$n \geq m$），求满足：$A(x)=B(x)D(x)+R(x),\deg(D)=n-m,\deg(R)<m$的多项式$D(x)$和$R(x)$。显然可以类比的认为它俩分别是商和余数。

• 首先考虑求$D(x)$，因为$R(x)$等会可以暴力减出来。那么我们就应该先消除$R(x)$的影响，因为这样如果没有它式子看起来就像求逆的式子一样简单了；
• 下面会有一个不知道怎么想出来的巧妙变形：
  
  1. 将$\frac{1}{x}$带入原式：$A(\frac{1}{x})=B(\frac{1}{x})D(\frac{1}{x})+R(\frac{1}{x})$；
  2. 两边同乘$x^n$：$x^nA(\frac{1}{x})=x^mB(\frac{1}{x}) \times x^{n-m}D(\frac{1}{x})+x^{m-1}R(\frac{1}{x}) \times x^{n-m+1}$；
• 注意到每一项都有一个形如$x^{\deg(F)}F(x)$的东西，我们观察一下这是什么操作：将$F(x)$的系数翻转，即变为$F^R(x)=\sum_{i=0}^{\deg(F)}[x^{\deg(F)-i}]F(x)x^i$。
• 这样式子改写为：$A^R(x)=B^R(x)D^R(x)+x^{n-m+1}R^R(x)$；
• 然后就非常有趣，把式子放在$\mod{x^{n-m+1}}$意义下并不会影响求$D(x)$，因为$\deg(D)=n-m$，而$R^R(x)$会消失！
• 那么只需要求一个逆就好了，复杂度就是$o(n\log{n})$的，这个变形还是非常神奇的。

代码：

inline void poly_div(int a[],int n,int b[],int m,int d[]){
    static int ra[MAXL],rb[MAXL],invrb[MAXL];
    if(m>n){
        d[0]=0;
        return;
    }
    int len=poly_len(n-m);
    memset(ra,0,len<<1<<2);
    for(int i=0;i<=n-m;++i)
        ra[i]=a[n-i];
    memset(rb,0,len<<1<<2);
    for(int i=0;i<=min(n-m,m);++i)
        rb[i]=b[m-i];
    poly_inv(rb,invrb,len);
    poly_mul(ra,n-m,invrb,n-m,d,n-m);
    reverse(d,d+n-m+1);
}
inline void poly_mod(int a[],int n,int b[],int m,int r[]){
    static int d[MAXL],tmp[MAXL];
    if(m>n){
        memcpy(r,a,n+1<<2);
        return;
    }
    poly_div(a,n,b,m,d);
    poly_mul(b,m,d,n-m,tmp,m-1);
    poly_minus(a,tmp,r,m-1);
}

• 多项式取模可以用在特征多项式优化线性递推里。

常数大到没法考的东西

多项式多点求值

给一个$n$次多项式$A(x)$和$m$个点$x_1,x_2,\cdots,x_m$，求$A(x_1),A(x_2),\cdots,A(x_m)$。假设$n,m$同阶（实际上复杂度瓶颈是关于$m$的）。

• 多项式多点求值其实我们每次写FFT都会用到，但FFT是带入$n$次单位根多点求值，我们可以用它的特殊性质进行分治。
• 这里我们其实也可以那样，想办法每次减半多项式的次数以及要求值的点的个数。
• 假设这层分治我们要求$x_{\text{l}},\cdots,x_{\text{r}}$的值，下层分治我们希望分别递归处理$x_{\text{l}},\cdots,x_{\text{mid}}$和$x_{\text{mid}+1},\cdots,x_{\text{r}}$，那么我们要想办法减半递归进去的多项式的次数；
• 一个很自然的想法就是将现在的多项式中关于要递归处理的部分无关的部分消去。
• 以处理左半边为例，设当前多项式为$A(x)$，要递归到左边的多项式为$A_{\text{l}}(x)$；
• 那么有$A(x)=D(x)\prod_{i=\text{l}}^{\text{mid}}(x-x_i)+A_{\text{l}}(x)$；
• 因为前面的部分当带入要求的点值时都会变为$0$，显然是不会造成贡献，那么我们多项式取模一下就好了。
• 因为每次取模后多项式次数为$\text{r}-\text{l}$，所以次数也是减半的，复杂度是：$T(n)=2T(\frac{n}{2})+O(n\log{n})$，即$O(n\log^2{n})$。
• 实现的话，可以先分治预处理出每层要模的多项式，递归到边界时，多项式只有常数项了，这就是答案。但是由于常数巨大，可以考虑递归到小范围时就暴力求，我写的代码不卡常$50000$的数据开O2要跑快2s了。

代码：

void poly_DC(int d,int l,int r,int x[],int dc[][MAXN<<1]){
    if(l==r){
        dc[d][l<<1]=P-x[l],dc[d][l<<1|1]=1;
        return;
    }
    int mid=l+r>>1;
    poly_DC(d+1,l,mid,x,dc),poly_DC(d+1,mid+1,r,x,dc);
    poly_mul(dc[d+1]+(l<<1),mid-l+1,dc[d+1]+(mid+1<<1),r-mid,dc[d]+(l<<1),r-l+1);
}
void poly_eval(int d,int l,int r,int f[],int len,int dc[][MAXN<<1],int y[]){
    static int eval[MAXD][MAXN];
    poly_mod(f,len,dc[d]+(l<<1),r-l+1,eval[d]+l);
    if(l==r){
        y[l]=eval[d][l];
        return;
    }
    int mid=l+r>>1;
    poly_eval(d+1,l,mid,eval[d]+l,r-l,dc,y);
    poly_eval(d+1,mid+1,r,eval[d]+l,r-l,dc,y);
}
inline void poly_eval(int f[],int len,int x[],int n,int y[]){
    static int dc[MAXD][MAXN<<1];
    poly_DC(0,0,n-1,x,dc);
    poly_eval(0,0,n-1,f,len,dc,y);
}

多项式快速插值

给定$n+1$个点值$(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)$，要求插值出一个$n$次多项式来。所谓插值，就是把这个值“插到”多项式里，也就是多项式在这个点的取值是给定值。

• 多项式插值常用的是拉格朗日插值：$\sum_{i=0}^n\frac{\prod_{j \neq i}(x-x_j)}{\prod_{j \neq i}(x_i-x_j)}y_i$。这个式子的原理就是对于每一个$i$的那一项，当式子带入的是其它的$n$个值时，这里会变成$0$，否则会变成$y_i$。
• 直接求显然是$O(n^2)$的，我们考虑优化它。
• 首先考虑如何对于每个$i$求出$\prod_{j \neq i}(x_i-x_j)$，也就是将$x_i$带入多项式$A_i(x)=\frac{\prod_{j=0}^n(x-x_j)}{x-x_i}$的值；
• 考虑到分子分母此时都为$0$，我们使用洛必达法则：
  $$A_i(x_i)=\lim_{x \to x_i}\frac{\prod_{j=0}^n(x-x_j)}{x-x_i}=\lim_{x \to x_i}\left[\prod_{j=0}^n(x-x_j)\right]'$$
• 注意到最后化出的这个多项式和带入的$x_i$无关，我们直接多点求值即可。
• 令$z_i=\frac{y_i}{\prod_{j \neq i}(x_i-x_j)}$，现在就是要求$\sum_{i=0}^n z_i \prod_{j \neq i}(x-x_j)$。
• 这个可以用分治来做，就是每次将左边的乘上$\prod_{i=\text{mid}+1}^{\text{r}}(x-x_i)$，右边的同理，然后加起来就可以了。类似之前的多点求值，分治预处理出每层要用的多项式即可。
• 复杂度显然也是$O(n\log^2{n})$的，但是常数几乎是多点求值的1.5倍...，怕是也没什么题可以出了。

代码：

void poly_interp(int d,int l,int r,int val[],int dc[][MAXN<<1],int f[]){
    if(l==r){
        f[l]=val[l];
        return;
    }
    int mid=l+r>>1;
    poly_interp(d+1,l,mid,val,dc,f);
    poly_interp(d+1,mid+1,r,val,dc,f);
    static int g[MAXN],h[MAXN];
    poly_mul(f+l,mid-l,dc[d+1]+(mid+1<<1),r-mid,g,r-l);
    poly_mul(f+mid+1,r-mid-1,dc[d+1]+(l<<1),mid-l+1,h,r-l);
    poly_add(g,h,f+l,r-l);
}
inline void poly_interp(int x[],int y[],int n,int f[]){
    static int dc[MAXD][MAXN<<1],tmp[MAXN],val[MAXN];
    poly_DC(0,0,n-1,x,dc);
    poly_diff(dc[0],n,tmp);
    poly_eval(tmp,n-1,x,n,val);
    for(int i=0;i<n;++i)
        val[i]=(long long)y[i]*inv(val[i])%P;
    poly_interp(0,0,n-1,val,dc,f);
}

ps：上面讲的多项式操作复杂度都做到了理论下界，但是FFT和NTT常数那么大，跟多一两个$\log$似的也是没谁了。