【算法讲15：多项式卷积】FFT | NTT

模板题

• 【模板】多项式乘法（FFT） | 洛谷 P3803
  求两个多项式的乘积

多项式的表示法

• 系数表示：$n$ 个系数确定一个 $n-1$ 次多项式
  $$(a_0,a_1,a_2,\cdots ,a_{n-1})$$
• 点值表示：$n$ 个点确定一个 $n-1$ 刺多项式
  $$\begin{gathered} (x_0,x_1,x_1,\cdots ,x_{n-1}) \\ (y_0,y_1,y_1,\cdots ,y_{n-1}) \end{gathered}$$
  其中满足
  $$y_i=\sum _{j=0}^{n-1}{a}_j{x}_i^j$$

多项式的乘法与卷积

多项式乘法

• 有两个多项式 $A(x)，B(x)$，两个多项式相乘结果为 $C(x)$。
• 新的多项式系数可以表示为：
  $$c_k=\sum _{i+j=k}{a}_i{b}_j=\sum _{i=0}^k{a}_i{b}_{k-i}$$
• 多项式乘法得到的结果又可以叫做这两个多项式的卷积，可以用 $FFT$ 进行快速计算。

$\lceil$快速傅里叶变换$\rfloor$$FFT$

• 如果暴力求解两个多项式的卷积，时间复杂度为 $O(nm)$
  利用 $FFT$ 可以把时间复杂度降到 $O(n\log n)$
  具体分为两个过程：
• 把 $A(x)、B(x)$ 的系数表示转化成点值表示，然后用 $O(n)$ 的方法求得 $C(x)$ 的点值表示，然后还原成 $C(x)$ 的系数表示。
• 那么第一步，我们想要找一个快速的方法把 $A(x)、B(x)$ 的点值表示求出来。

单位根

• $FFT$ 的方法是找到了一组特殊的点代到多项式当中去，可以快速找到一组点的数值。那组特殊的点就是复数中的单位根。
• 一个复数可以在复平面上用一个向量来表示。在单位元上的向量都可以成为单位根。
  我们把一个单位元分成 $n$ 份，我们就得到了 $n$ 个单位根，其中的第 $k$ 个可以即为 ${\omega }_n^k$，并且 ${\omega }_n^0=1$
• 欧拉公式：$e^{i\theta }=cos\theta +i\cdot \sin \theta$
  于是，
  $${\omega }_n^k=e^{2\pi \frac{k}{n}i}$$
• 性质：
  ${\omega }_n^k={\omega }_{2n}^{2k}$
  ${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
  ${\omega }_n^k\cdot {\omega }_n^l={\omega }_n^{k+l}$

离散傅里叶变换 $DFT$

• 考虑 $n=2^b$ 项的多项式 $A(x)$ 系数表示为 $(a_0,a_1,\cdots ,a_{n-1})$
• 我们现在想导入一组单位根 $({\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1})$，快速算出 $(A({\omega }_n^0),A({\omega }_n^1),\cdots ,A({\omega }_n^{n-1}))$，这个过程成为 $DFT$
• 首先，我们把多项式按奇偶性分类：
  $$A(x)=(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+x(a_1+a_3{x}^2+\cdots +a_{n-1}{x}^{n-2})$$
  现在，我们设两个多项式：
  $$\begin{gathered} A0(x)=a_0+a_2{x}^1+\cdots +a_{n-2}{x}^{\frac{n-2}{2}} \\ A1(x)=a_1+a_3{x}^1+\cdots +a_{n-1}{x}^{\frac{n-2}{2}} \end{gathered}$$
  于是，原来的多项式就变为：
  $$A(x)=A0(x^2)+xA1(x^2)$$

$DFT$ 的详细讨论

• 当 $0\le k\le \frac{n}{2}-1$ 时，
  $$A({\omega }_n^k)=A0(({\omega }_n^k{)}^2)+{\omega }_n^k\cdot A1(({\omega }_n^k{)}^2)$$
  化简得到：
  $$A({\omega }_n^k)=A0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\cdot A1({\omega }_{\frac{n}{2}}^k)$$
  可以看到，我们分为了两个子问题，可以递归求解。
• 那么 $\frac{n}{2}\le k+\frac{n}{2}\le n-1$，我们带入得到：
  $$A({\omega }_n^{k+\frac{n}{2}})=A0(({\omega }_n^{k+\frac{n}{2}}{)}^2)+{\omega }_n^{k+\frac{n}{2}}\cdot A1(({\omega }_n^{k+\frac{n}{2}}{)}^2)$$
  化简得到：
  $$A({\omega }_n^{k+\frac{n}{2}})=A0({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\cdot A1({\omega }_{\frac{n}{2}}^k)$$
• 时间复杂度：$T(n)=2T(\frac{n}{2})+O(n)=O(n\log n)$

离散傅里叶逆变换 $IDFT$

• 将 $DFT$ 得到的点值转化为系数，这一过程就叫做 $IDFT$
• 假设我们得到了 $A(x)$ 的点值，设为 $(d_0,d_1,\cdots ,d_{n-1})$，我们要获得其多项式 $F(x)$
  设 $c(k)$ 为 $x={\omega }_n^{-k}$ 时的点值。
  $$c_k=\sum _{i=0}^{n-1}{d}_i\cdot ({\omega }_n^{-k}{)}^i=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j\cdot ({\omega }_n^i{)}^j\cdot ({\omega }_n^{-k}{)}^i=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=\sum _{j=0}^{n-1}{a}_j\cdot n\cdot [j=k]=n\cdot a_k$$

减小常数 (蝴蝶变换)

• 我们一遍遍递归下去，在最后一层发现一个规律：原序列和置换后的序列的对应二进制位刚好是翻转的。比如 $1\to 8$，$14\to 7$
• 代码：
  $2.77s$

const int MAX = 1e7+50;

const double PI = acos(-1);
struct comp{
    double x,y;
    comp(double xx = 0,double yy = 0){x = xx,y = yy;}
};
comp operator +(comp a,comp b){return comp(a.x+b.x , a.y+b.y);}
comp operator -(comp a,comp b){return comp(a.x-b.x , a.y-b.y);}
comp operator *(comp a,comp b){return comp(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}

int lim = 1;
int l,r[MAX];
void FFT(comp *A,int type){
    for(int i = 0;i < lim;++i)
        if(i < r[i])swap(A[i],A[r[i]]);
    for(int mid = 1;mid < lim;mid<<=1){
        comp omg(cos(PI/mid),type * sin(PI/mid));
        for(int R = mid<<1,j = 0;j < lim;j += R){
            comp w(1,0);
            for(int k = 0;k < mid;++k,w = w * omg){
                comp x = A[j+k],y = w * A[j + mid + k];
                A[j + k] = x + y;
                A[j + mid + k] = x - y;
            }
        }
    }
}
comp aa[MAX],bb[MAX];
int main()
{
    int N,M;N = read();M = read();
    for(int i = 0;i <= N;++i)aa[i].x = read();
    for(int i = 0;i <= M;++i)bb[i].x = read();
    while(lim <= N + M)lim<<=1,l++;
    for(int i = 0;i < lim;++i)
        r[i] = (r[i>>1]>>1) | ((i&1)<<(l-1));
    FFT(aa,1);
    FFT(bb,1);
    for(int i = 0;i <= lim;++i)aa[i] = aa[i] * bb[i];
    FFT(aa,-1);
    for(int i = 0;i <= N + M;++i)
        printf("%d ",(int)(aa[i].x/lim+0.5));
    return 0;
}

$\lceil$快速数论变换$\rfloor$ (NTT)

• 求的是多项式乘法的系数取模之后的结果。这里模数应该要是一个质数。
• 我们希望带入一组幂次以便于更好取模运算。

原根

• 如果 $g^i\mathrm{mod}p(1\le i\le p-1)$ 的值互不相同，那么 $g$ 称为 $p$ 的原根。
  并非所有的数都有原根，且如果一个数存在原根，其个数也不唯一。
• 事实上，原根存在的重要条件是：$m=2,4,p^a,2p^a(a>0)$ 且 $a$ 为奇素数。
  如果这个数有原根，那么一定有 $\phi (\phi (m))$ 个原根。
• 原根没有快速求法，只能通过暴力枚举判断，或者查表。
  常见的 $NTT$ 模数是 $998244353$，其原根为 $3$

同样的性质

• 加入模数为 $p$，我们发现构造一下原根，得到与单位根有同样的性质的数，即：
  $${\omega }_n^i=g^{\frac{i(p-1)}{n}}\mathrm{mod}p$$
  直接将 $FFT$ 中的原根替换为 $NTT$ 中的原根。
• 代码
  $1.64s$

const int MAX = 1e7+50;

ll qpow(ll a,ll n){/* */ll res = 1LL;while(n){if(n&1)res=res*a%MOD;a=a*a%MOD;n>>=1;}return res;}

const int P = 998244353,G = 3,Gi = 332748118;
int lim = 1;
int l,r[MAX];
void NTT(ll *A,int type){
    for(int i = 0;i < lim;++i)
        if(i < r[i])swap(A[i],A[r[i]]);
    for(int mid = 1;mid < lim;mid<<=1){
        ll omg = qpow(type == 1 ? G : Gi , (P - 1) / (mid << 1));
        for(int j = 0;j < lim;j += (mid << 1)){
            ll w = 1;
            for(int k = 0;k < mid;++k,w = w * omg % P){
                int x = A[j+k],y = w * A[j + mid + k] % P;
                A[j + k] = (x + y) % P;
                A[j + mid + k] = (x - y + P) % P;
            }
        }
    }
}
ll aa[MAX],bb[MAX];
int main()
{
    int N,M;N = read();M = read();
    for(int i = 0;i <= N;++i)aa[i] = (read() + P) % P;
    for(int i = 0;i <= M;++i)bb[i] = (read() + P) % P;
    while(lim <= N + M)lim<<=1,l++;
    for(int i = 0;i < lim;++i)
        r[i] = (r[i>>1]>>1) | ((i&1)<<(l-1));
    NTT(aa,1);
    NTT(bb,1);
    for(int i = 0;i < lim;++i)aa[i] = (aa[i] * bb[i]) % P;
    NTT(aa,-1);
    ll iv = qpow(lim,P - 2);
    for(int i = 0;i <= N + M;++i)
        printf("%d ",(aa[i] * iv) % P);
    return 0;
}

参考与查阅

• 【SDUACM-暑期专题div1】FFT、NTT
• 洛谷改题的 $dalao$ 们的题解