发现自己的多项式水平菜到不行，所以写篇文章记录一下。

• 基础操作
  • FFT & NTT
  • 分治 FFT & NTT
  • 拉格朗日插值
  • 多项式求导与积分
  • 多项式求逆
  • 多项式开根
  • 多项式带余除法
  • 多项式对数函数
  • 多项式指数函数
  • 参考代码

基础操作

FFT & NTT

基础知识。

考虑如何快速计算多项式乘法，发现常规的做法是没有前途的，考虑换种方法计算。

有一种计算方法是这样的（设要计算 $A\times B$ ）：

• 将 $A$ 转化为点值表示。
• 将 $B$ 转化为点值表示。
• $O(n)$ 地计算点值表示的乘积 $C$。
• 将 $C$ 转化为数值表示。

复杂度瓶颈为点值表示与数值表示的转化，朴素做是 $O(n^2)$。

为了加速这个过程，FFT 引入了单位复根 $\omega$，NTT 引入了原根 $g$。

单位复根 $\omega$ 是 $x^n=1$ 的解集，其满足如下性质：

• $\omega^n_n=1$；
• $\omega^k_n=\omega^{2k}_{2n}$；
• $\omega^{k+n}_{2n}=-\omega^k_{2n}$。

考虑分治，对一个多项式 $f(x)$ 奇偶性分组处理，比如一个多项式：

$f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7$

考虑奇偶性分组后的状态：

$f(x)=a_0+a_2x^2+a_4x^4+a_6x^6+x(a_1+a_3x^2+a_5x^4+a_7x^6)$

建立新函数：

$G(x)=a_0+a_2x+a_4x^2+a_6x^3\\H(x)=a_1+a_3+a_5x^2+a_7x^3$

那么待求解的 $f(x)$ 可改为：

$f(x)=G(x^2)+xH(x^2)$

带入 $\omega$ 可得：

$\begin{aligned} \operatorname{DFT}(f(\omega_n^k)) &=\operatorname{DFT}(G((\omega_n^k)^2)) + \omega_n^k \times \operatorname{DFT}(H((\omega_n^k)^2))\\ &=\operatorname{DFT}(G(\omega_n^{2k})) + \omega_n^k \times \operatorname{DFT}(H(\omega_n^{2k}))\\ &=\operatorname{DFT}(G(\omega_{n/2}^k)) + \omega_n^k \times \operatorname{DFT}(H(\omega_{n/2}^k)) \end{aligned} \\ \begin{aligned} \operatorname{DFT}(f(\omega_n^{k+n/2})) &=\operatorname{DFT}(G(\omega_n^{2k+n})) + \omega_n^{k+n/2} \times \operatorname{DFT}(H(\omega_n^{2k+n}))\\ &=\operatorname{DFT}(G(\omega_n^{2k})) - \omega_n^k \times \operatorname{DFT}(H(\omega_n^{2k}))\\ &=\operatorname{DFT}(G(\omega_{n/2}^k)) - \omega_n^k \times \operatorname{DFT}(H(\omega_{n/2}^k)) \end{aligned}$

然后就可以分治做了。

诶等等，递归常数特别大呢！

为了不递归，我们考虑一开始就预处理出来每一项最后到了什么地方，这个做法是位逆序置换（又名蝴蝶置换）。

考虑找规律，发现将下标二进制表示后翻转即可，设 $R_i$ 为 $i$ 变到的地方，那么可以考虑这样求：

• 右移一位。
• 翻转。
• 再右移一位。
• 处理最高位。

那么 $R_i$ 的递推式可以轻易的写成：

$R_i=\lfloor\frac{R_{\lfloor\frac{i}{2}\rfloor}}{2}\rfloor+(i\bmod2)\times\frac{len}{2}$

以上就是如何将数值表示转化为点值表示了，所以怎样转回来呢？

直接反转后除以 $len$，或是取单位根的倒数跑即可，具体推导过程见 OI-Wiki。

模板：

void FFT(int n,cp* x,int op) {
    for(int i=0;i<n;i++)
        if(i<rk[i]) swap(x[i],x[rk[i]]);
    for(int mid=1;mid<n;mid*=2) {
        cp W(cos(Pi/mid),op*sin(Pi/mid));
        for(int R=mid*2,j=0;j<n;j+=R) {
            cp w(1,0);
            for(int k=0;k<mid;k++,w=w*W) {
                cp z1=x[j+k],z2=w*x[j+mid+k];
                x[j+k]=z1+z2,x[j+mid+k]=z1-z2;                
            }						            
        }				        
    }	    
}

NTT 呢？

首先解释 $g$，$g$ 是原根，在 $\bmod p$ 的情况下 $\forall i\not =j,g^i\not\equiv g^j\pmod p$。

我们用 $g^{\frac{\varphi(p)}{n}}$ 代替 $\omega$ 进行计算，由 $\varphi(p)=p-1$ 且 $n=2^q$，所以强制要求 $p=k\times 2^q+1$，其中 $2^q\ge n$。

较 FFT 精度高，但却对模数有着致命的要求。

时间复杂度 $O(n\log n)$。

分治 FFT & NTT

问题：求 $f_n$，其中：

$f_i=\sum^{i=1}_{j=0} f_jg_{i-j}$

考虑 CDQ 分治，对于一个区间 $[l,r]$，将其划分成 $[l,mid],[mid+1,r]$ 两段。

先处理 $[l,mid]$，再计算 $[l,mid]$ 对 $[mid+1,r]$ 的贡献，再处理 $[mid+1,r]$。

中间的计算使用 NTT/FFT 解决，时间复杂度 $O(n\log^2n)$

拉格朗日插值

问题，给定 $n$ 个点 $(x_i,y_i)$，确定一个多项式 $y=f(x)$，求出 $f(k)$ 的值。

显然的事实：

$f(x)\equiv f(y)\pmod{x-y}$

证明：$f(x)-f(y)=(a_0-a_0)+\color{red}a_1(x-y)\color{black}+a_2(x^2-y^2)+a_3(x^3-y^3)+\cdots$。

那么我们就可以列一个线性同余方程组：

$\begin{cases}f(x)\equiv y_1\pmod {x-y_1}\\f(x)\equiv y_2\pmod{x-y_2}\\\ldots\\f(x)\equiv y_n\pmod{x-y_n}\end{cases}$

直接 CRT，所有模数的积为：

$\prod_{i=1}^n (x-y_i)$

那么 $m_i=\prod _{j\not=i} (x-y_i),m_i^{-1}=\prod _{j\not=i} \frac {1}{x_i-x_j}$。

所以：

$f(x)=\sum^n_{i=1} y_i \prod_{j\not=i} \frac {x-x_j}{x_i-x_j}$

朴素实现复杂度是 $O(n^2)$ 的。

在所取点值连续的时候，即 $x_i=i$ 的时候，可以直接用前缀积优化到 $O(n)$。

在需要动态加点的时候，可以使用重心拉格朗日插值法：

回顾之前的式子：

$f(x)=\sum^n_{i=1} y_i \prod_{j\not=i} \frac {x-x_j}{x_i-x_j}$

设 $g=\prod^n_{i=1} k-x_i$，那么：

$f(x)=g\sum^n_{i=0}\prod\frac {y_i}{(k-x_i)(x_i-x_j)}$

设 $t_i=\frac {y_i}{\prod_{j\not=i} x_i-x_j}$，那么：

$f(x)=g\sum^n_{i=0}\frac {t_i}{k-x_i}$

加入一个点只需计算 $t_i$ 即可。

求自然数幂和可以直接拉格朗日插值，时间复杂度 $O(n)$。

多项式求导与积分

由于求导和积分满足可加性，所以拆成单项式分别求导或积分。

如果你学过微积分相关理论，会知道：

$A'(x)=\sum^n_{i=1}ia_ix^{i-1}\\\int A(x)dx=(\sum^n_{i=0}\frac {a_i}{i+1}x^{i+1})+c$

直接做即可，为数不多的 $O(n)$ 时间复杂度。

多项式求逆

给定 $G(x)$，求一个 $F(x)$ 使得 $F(x)G(x)\equiv0\pmod {x^n}$。

首先假设我们已经求出了一个 $f(x)$ 使得 $f(x)G(x)\equiv 0\pmod{x^{n/2}}$，考虑如何求出 $F(x)$。

因为 $F(x)G(x)\equiv 0\pmod{x^n}$，所以必然 $F(x)G(x)\equiv 0 \pmod{x^{n/2}}$，两式相减，得 $F(x)-f(x)\equiv0\pmod{x^{n/2}}$，此处 $G(x)$ 被约掉了。

两边同时平方，得：

$F^2(x)-2F(x)f(x)+f^2(x)\equiv0\color{red}{\pmod{x^{n}}}$

注意标红的部分，想一想，为什么？

考虑卷积的定义，在 $[0,n/2]$ 的部分每一项都是 $0$，而右边 $[n/2,n]$ 的部分由 $a_i=\sum^i_{j=1} a_ja_{i-j}$ 得来，容易发现，$j$ 与 $i-j$ 中必有一个小于零。

两边同时乘以 $G(x)$ 并移项，得：

$F(x)\equiv 2f(x)-G(x)f^2(x)\pmod{x^n}$

直接使用 FFT 计算即可，时间复杂度 $O(n\log n)$。

多项式开根

给定 $g(x)$，求 $f(x)$，使得：$f^2(x)\equiv g(x)\pmod{x^n}$。

类似求逆，假设已经求出了 $f_0^2(x)\equiv g(x)\pmod{x^{n/2}}$。

移项得：

$f_0^2(x)-g(x)\equiv 0\pmod{x^{n/2}}$

两边平方得：

$(f_0^2(x)-g(x))^2\equiv 0\pmod {x^n}$

使用简单的数学公式 $(x-y)^2=(x+y)^2-4xy$ 得：

$(f_0^2(x)+g(x))^2\equiv 4f_0^2(x)g(x)\pmod {x^n}$

将 $4f_0^2(x)$ 移到左边，得：

$(\frac{f_0^2(x)+g(x)}{2f_0(x)})^2\equiv g(x) \pmod{x^n}$

左右两边同时开根，得：

$\frac{f_0^2(x)+g(x)}{2f_0(x)} \equiv f(x)\pmod{x^n}$

$2^{-1}f_0(x)+2^{-1}f_0^{-1}(x)g(x)\equiv f(x)\pmod{x^n}$

直接计算即可，时间复杂度 $O(n\log n)$。

多项式带余除法

问题：给定 $n$ 次多项式 $F(x)$ 和 $m$ 次多项式 $G(x)$，求 $Q(x)$ 与 $R(x)$，满足：

• $\deg Q(x)=n-m$，$\deg R(x)<m$；
• $F(x)=Q(x)\times G(x)+R(x)$。

对于一个多项式 $f(x)$，定义 $f_R(x)=x^nf(\frac 1 x)$，等价于将多项式 $f$ 的系数翻转一次。

那么：

$F(x)=Q(x)\times G(x)+R(x)\\F(\frac 1 x)=Q(\frac 1 x)\times G(\frac 1 x)+R(\frac 1 x)\\x^nF(\frac 1 x)=x^{n-m} Q(\frac 1 x)\times x^mG(\frac 1 x)+x^{n-m+1}x^{m-1}R(\frac 1 x)\\F_R(x)=Q_R(x)\times G_R(x)+x^{n-m+1}R_R(x)\\F_R(x)\equiv Q_R(x)\times G_R(x)\pmod {x^{n-m+1}}\\Q_R(x)\equiv F_R(x)\times G_R^{-1}(x) \pmod{x^{n-m+1}}$

直接求出 $Q(x)$，那么回代可得 $R(x)$，时间复杂度 $O(n\log n)$。

多项式对数函数

求 $B(x)\equiv \ln A(x)\pmod{x^n}$。

设 $F(x)=\ln x$，那么：

$B(x)=F(A(x))\\B'(x)=F'(A(x))A'(x)\\B'(x)=\frac {A'(x)}{A(x)}\\B(x)=\int \frac {A'(x)dx}{A(x)}$

直接模拟即可，时间复杂度 $O(n\log n)$。

多项式指数函数

求 $F(x)\equiv e^{A(x) \pmod{x^n}}$。

考虑多项式牛顿迭代。

多项式牛顿迭代解决的问题是，给定函数 $G(x)$，求 $G(F(x))\equiv 0\pmod {x^n}$。

既然是迭代，那么当 $n=1$ 的时候是要单独求出的，这里不讲。

类似求逆，假设已经求出 $G(F_0(x))\equiv 0 \pmod{x^{\lceil \frac n x\rceil}}$，考虑如何扩展。

把 $G(F(x))$ 在 $F_0(x)$ 处泰勒展开：

$G(F(x))=G(F_0(x))+\frac {G'(F_0(x))}{1!}(F(x)-F_0(x))+\frac {G''(F_0(x))}{2!}(F(x)-F_0(x))+\ldots$

因为 $F(x)$ 和 $F_0(x)$ 的最后几项相同，所以 $(F(x)-F_0(x))^2$ 的最低非 $0$ 次数大于 $n$，所以可以得到：

$G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\pmod{x^n}\\G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\equiv 0\pmod{x^n}\\G(F_0(x))+G'(F_0(x))F(x)-G'(F_0(x))F_0(x)\equiv 0\pmod{x^n}\\G'(F_0(x))F(x)\equiv G'(F_0(x))F_0(x)-G(F_0(x))\pmod {x^n}\\F(x)\equiv F_0(x)-\frac {G(F_0(x))}{G'(F_0(x))} \pmod {x^n}$

那么，将 $F(x)=e^{A(x)}$ 变形得 $\ln F(x)-A(x)=0$。

设 $G(F(x))=\ln F(x)-A(x)$，那么 $G'(F(x))=\frac 1 {F(x)}$，带入牛顿迭代公式得：

$F(x)\equiv F_0(x)-\frac{\ln F_0(x)-A(x)}{F_0(x)}\pmod{x^n}\\F(x)\equiv F_0(x)(1-\ln F_0(x)+A(x))\pmod{x^n}$

直接做，时间复杂度 $O(n\log n)$。

参考代码

不想太长，所以放一份在剪贴板。