（原创）浅析卷积

## 0. 写在前面

从粗斜体到分界线中的内容可以跳过。

本文中一些算法的简历：

| 简写 | 全称 | 中文 | ~~搞笑~~ |
|:--:|:--:|:--:|:--:|
| FWT | Fast Walsh Transformation | 快速沃尔什变换 | Fast Wonderful TLE |
| FFT | Fast Fourier Transformation | 快速傅里叶变换 | Fantastic Fluency TLE |
| NTT | Number Theory Transformation | 快速数论变换 | Natural Talented TLE|

## 1. 卷积的定义

$$ c_i= \sum_{j*k=i}{a_j\times b_k} $$

$ c $ 称为 $ a,b $ 的卷积，当 $ * $ 指不同的运算符时， $ c $ 有不同的求法，现在分类讨论。

## 2. 当 $ * $ 指 $ max/min $

可以通过简单的前/后缀和计算，以下是 $ max $ 卷积的代码：

``` cpp
#include<stdio.h>
int N,a[1000005],b[1000005],c[1000005],A[1000005],B[1000005],C[1000005];
int main(){
  scanf("%d",&N);
  for(int i=1;i<=N;i++)
    scanf("%d",a+i);
  for(int i=1;i<=N;i++)
    scanf("%d",b+i);
  for(int i=1;i<=N;i++){
    A[i]=a[i]+A[i-1];
    B[i]=b[i]+B[i-1];
  }
  for(int i=1;i<=N;i++)
    C[i]=A[i]*B[i];
  for(int i=N;i>=1;i--)
    c[i]=C[i]-C[i-1];
  for(int i=1;i<=N-1;i++)
    printf("%d ",c[i]);
  printf("%d\n",c[N]);
  return 0;
}
```

$ min $ 卷积与之类似。

**为什么它是正确的？**

讨论 $ max $ 卷积

$$ c_i=\sum_{max(j,k)=i}{a_j\times b_k} $$

即

$$ c_i=a_i\times b_i+\sum_{j=1}^{i-1}{(a_j\times b_i+b_j\times a_i)} $$

程序中

$$ A_i=\sum_{j=1}^{i}{a_j}\text{ }\text{ }\text{ }\text{ }B_i=\sum_{j=1}^{i}{b_j} $$

所以

$$ C_i=\sum_{j=1}^{i}{\sum_{k=1}^{i}{a_j\times b_k}} $$

易得

$$ c_i=C_i-C_{i-1}=\sum_{j=1}^{i}{(a_j\times b_i)}+\sum_{j=1}^{i-1}{(a_i\times b_j)} $$

与理论答案相符

---

由此我们可以总结出一点经验，求卷积的流程往往是这样：

1. 用某种变换将$ a_i,b_i $变成$ A_i,B_i $

2. $ C_i=A_i\times B_i $

3. 用其逆变换将 $ C_i $ 变成 $ c_i $ 得到答案

## 3. 当 $ * $ 指 $ \vee/\wedge $

**在这篇文章中，如果涉及到带 $lg_N$ 复杂度的卷积变换， $N=2^k(k\in\mathbb{N})$ 。实际实现时在高位补 $0$**

### 3.1 用向量表示数

一个 $ k $ 位二进制数可以表示成一个 $ k $ 维向量。例如当 $ k=3 $ 时：

| 数值 | 向量表示 | 数值 | 向量表示 |
| :----------: | :----------: | :----------: | :----------: |
| 0 | $ \begin{Bmatrix} 0 & 0 & 0 \end{Bmatrix} $ | 4 | $ \begin{Bmatrix} 1 & 0 & 0 \end{Bmatrix} $ |
| 1 | $ \begin{Bmatrix} 0 & 0 & 1 \end{Bmatrix} $ | 5 | $ \begin{Bmatrix} 1 & 0 & 1 \end{Bmatrix} $ |
| 2 | $ \begin{Bmatrix} 0 & 1 & 0 \end{Bmatrix} $ | 6 | $ \begin{Bmatrix} 1 & 1 & 0 \end{Bmatrix} $ |
| 3 | $ \begin{Bmatrix} 0 & 1 & 1 \end{Bmatrix} $ | 7 | $ \begin{Bmatrix} 1 & 1 & 1 \end{Bmatrix} $ |

### 3.2 $\vee$的实质

我们把用向量表示的数字 $ \vee $ ,例如

$$ 3 \vee 6= \begin{Bmatrix} 0 & 1 & 1 \end{Bmatrix} \vee \begin{Bmatrix} 1 & 1 & 0 \end{Bmatrix}= \begin{Bmatrix} 1 & 1 & 1 \end{Bmatrix}=7 $$

由此可知 $ \vee $ 的本质是按位 $ max $ ，所以 $ \vee $ 卷积的变换就是按位前缀和。它的逆变换其实是一个脑筋急转弯，只要把循环倒过来，把$+=$改成$-=$就可以了。代码如下：

``` cpp
#include<stdio.h>
int n,k,N,a[100005],b[100005],c[100005];
int main(){
  scanf("%d",&n);
  for(int i=0;i<n;i++)
    scanf("%d",a+i);
  for(int i=0;i<n;i++)
    scanf("%d",b+i);
  for(N=1;N<n;N<<=1,k++);
  for(int i=0;i<k;i++)
    for(int j=0;j<N;j++)
      if((j&(1<<i))==0)
        a[j+(1<<i)]+=a[j];
  for(int i=0;i<k;i++)
    for(int j=0;j<N;j++)
      if((j&(1<<i))==0)
        b[j+(1<<i)]+=b[j];
  for(int i=0;i<N;i++)
    c[i]=a[i]*b[i];
  for(int i=k-1;i>=0;i--)
    for(int j=N-1;j>=0;j--)
      if((j&(1<<i))==0)
        c[j+(1<<i)]-=c[j];
  for(int i=0;i<=N-2;i++)
    printf("%d ",c[i]);
  printf("%d\n",c[N-1]);
  return 0;
}
```

$ \wedge $ 卷积与之类似。

## 4. 当 $ * $ 指 $ \bigoplus $

### 4.1 千里之行，始于足下：N=2

由小学知识可知：
$$ <a_0,a_1> \Rightarrow <a_0+a_1,a_0-a_1> $$
$$ <A_0,A_1> \Rightarrow <\frac{A_0+A_1}{2},\frac{A_0-A_1}{2}> $$

我们将 $ a_i,b_i $ 带入
$$ C_0=A_0\times B_0=(a_0+a_1)\times (b_0+b_1)=a_0\times b_0+a_0\times b_1+a_1\times b_0+a_1\times b_1 $$
$$ C_1=A_1\times B_1=(a_0-a_1)\times (b_0-b_1)=a_0\times b_0-a_0\times b_1-a_1\times b_0+a_1\times b_1 $$
$$ c_0=\frac{C_0+C_1}{2}=a_0\times b_0+a_1\times b_1 $$
$$ c_1=\frac{C_0-C_1}{2}=a_0\times b_1+a_1\times b_0 $$

听说是正确的，于是拓展到高维
$$ A_i=\sum_{j=0}^{N}{(-1)^{bitcount(i\wedge j)\times a_j}} $$

### 4.2 一句不是废话的废话

虽然我们已经发现了正解，我们还是从向量的角度看一下。比如
$$ 3 \bigoplus 6= \begin{Bmatrix} 0 & 1 & 1 \end{Bmatrix} \bigoplus \begin{Bmatrix} 1 & 1 & 0 \end{Bmatrix}= \begin{Bmatrix} 1 & 0 & 1 \end{Bmatrix}=5 $$

再写一遍
$$ 3 \bigoplus 6= \begin{Bmatrix} 0 & 1 & 1 \end{Bmatrix} \bigoplus \begin{Bmatrix} 1 & 1 & 0 \end{Bmatrix} = \begin{Bmatrix} (0+1) \% 2 & (1+1) \% 2 & (1+0) \% 2 \end{Bmatrix} = \begin{Bmatrix} 1 & 0 & 1 \end{Bmatrix} =5 $$

发现 $ \bigoplus $ 其实是二进制无进位加法。

### 4.3 具体实现

$ \bigoplus $ 卷积的变换代码与 $ \vee/\wedge $ 卷积的代码略有不同，但 $ \vee/\wedge $ 卷积的代码也可以写成这种形式，具体题意见 [洛谷P4717【模板】快速沃尔什变换](https://www.luogu.org/problemnew/show/P4717) 。

``` cpp
#include<stdio.h>
const int p=998244353;
inline int power(int a,int k){
  int ans=1;
  for(;k;a=1LL*a*a%p,k>>=1)
    if(k&1)
      ans=1LL*ans*a%p;
  return ans;
}
int len,N,a[4][300005],b[4][300005];
void FWT(int* a,int op,int flag){
  int x,y;
  if(op==3&&flag==-1){
    for(int i=N>>1;i>0;i>>=1)
      for(int j=0;j<N;j+=i<<1)
        for(int k=0;k<i;k++){
          x=a[j+k];
          y=a[i+j+k];
          a[j+k]=1LL*(x+y)*power(2,p-2)%p;
          a[i+j+k]=(1LL*(x-y)*power(2,p-2)%p+p)%p;
        }
    return;
  }
  for(int i=1;i<N;i<<=1)
    for(int j=0;j<N;j+=i<<1)
      for(int k=0;k<i;k++){
        x=a[j+k];
        y=a[i+j+k];
        if(op==1){
          if(flag==1){
            a[j+k]=x;
            a[i+j+k]=(x+y)%p;
          }
          else{
            a[j+k]=x;
            a[i+j+k]=((y-x)%p+p)%p;
          }
        }
        if(op==2){
          if(flag==1){
            a[j+k]=(x+y)%p;
            a[i+j+k]=y;
          }
          else{
            a[j+k]=((x-y)%p+p)%p;
            a[i+j+k]=y;
          }
        }
        if(op==3){
          a[j+k]=(x+y)%p;
          a[i+j+k]=((x-y)%p+p)%p;
        }
      }
}
int main(){
  scanf("%d",&len);
  N=power(2,len);
  for(int i=0;i<N;i++){
    scanf("%d",a[1]+i);
    a[3][i]=a[2][i]=a[1][i];
  }
  for(int i=0;i<N;i++){
    scanf("%d",b[1]+i);
    b[3][i]=b[2][i]=b[1][i];
  }
  for(int j=1;j<=3;j++){
    FWT(a[j],j,1);
    FWT(b[j],j,1);
  }
  for(int i=0;i<N;i++)
    for(int j=1;j<=3;j++)
      a[j][i]=1LL*a[j][i]*b[j][i]%p;
  for(int j=1;j<=3;j++)
    FWT(a[j],j,-1);
  for(int j=1;j<=3;j++){
    for(int i=0;i<N-1;i++)
      printf("%d ",a[j][i]);
    printf("%d\n",a[j][N-1]);
  }
  return 0;
}
```

## 5. 当 $ * $ 指 $ + $

~~公式恐惧症患者请果断按下`Ctrl+W`以发起正当防卫~~

### 5.1 求多项式乘法的新方法

定义 $ N $ 次多项式 $ g,h $ ，我们可以选取 $ x_{0...2\times N} $ 带入 $ g,h $ 得到 $ G_{0...2\times N},H_{0...2\times N} $ ,将G和H逐位相乘得到 $ F $ ，最后将 $ F $ 消元得到 $ f $ , $ f=g\times h $。

不幸的是，带入多项式需要 $ O(n^2) $ 。What's worse,高斯消元需要 $ O(n^3) $ 。

### 5.2 选择带入的数

显然，瓶颈在带入的数的选择上。那我们需要带入怎样的数带入呢？

#### 5.2.1 复数

我们知道， $ x^2=-1 $ 无实数解，但我们定义 $ i^2=-1 $ 。

复数是所有能够写成 $ a+i\times b(a,b\in \mathbb{R}) $ 的数的集合，该集合记作 $ \mathbb{C} $ 。

把 `1` 看做 $ x $ 轴，把 `i` 看做 $ y $ 轴，就可以得到一个平面， ~~人称~~ 也就是复平面，复平面上的点和 $ \mathbb{C} $ 中的数一一对应。

#### 5.2.2 单位复根

单位圆指在复平面上以原点为圆心，以 $ 1 $ 为半径的圆。 $ \omega_N^k $ 指 $ 1 $ 对应的点绕原点旋转 $ \frac{k}{N} $ 圈得到的点所对应的数。特别的，我们将 $ \omega_N^k $ 称为单位复根。单位复根有以下性质。

#### 5.2.3 单位复根的性质

1. 对于任意整数 $ n\geqslant0,k\geqslant0,d\geqslant0 $ 有

$$ \omega_{N}^{k}=\omega_{Nd}^{kd} $$

2. 对于任意整数 $ n\geqslant0,k\geqslant0 $ 有

$$ \omega_{N}^{k}=\omega_{2N}^{2k} $$

3. 对于任意整数 $ n\geqslant0 $ 有

$$ \omega_{N}^{\frac{N}{2}}=-1 $$

4. 对于任意整数 $ n\geqslant0,i\geqslant0,j\geqslant0 $ 有

$$ \omega_{N}^{i+j}=\omega_{N}^{i}\times\omega_{N}^{j} $$

5. 对于任意整数 $ n\geqslant0,i\geqslant0,j\geqslant0 $ 有

$$ \omega_{N}^{ij}=(\omega_{N}^{i})^j $$

### 5.3 开始带入!

我们定义将 $ <\omega_N^0,\omega_N^1, ... ,\omega_N^{N-1}> $ 带入 $ <a_0,a_1, ... ,a_{N-1}> $ 的结果设为 $ <A_0,A_1, ... ,A_{N-1}> $ 。

根据定义得

$$ A_i=\sum_{j=0}^{N-1}a_j\times(\omega_N^i)^j $$

由性质5得

$$ A_i=\sum_{j=0}^{N-1}a_j\times\omega_N^{ij} $$

我们把 $ i $ 的定义域从 $ [0,N-1] $ 变成 $ [0,\frac{N}{2}-1] $

$$ A_i=\sum_{j=0}^{N-1}a_j\times\omega_N^{ij} $$

$$ A_{i+\frac{N}{2}}=\sum_{j=0}^{N-1}a_j\times\omega_N^{(i+\frac{N}{2})\times j} $$

奇偶分类

$$ A_i=\sum_{j=0}^{\frac{N}{2}-1}a_{2j}\times\omega_N^{2ij}+\sum_{j=0}^{\frac{N}{2}-1}a_{2j+1}\times\omega_N^{i(2j+1)} $$

$$ =\sum_{j=0}^{\frac{N}{2}-1}a_{2j}\times\omega_{\frac{N}{2}}^{ij}+\omega_N^i\sum_{j=0}^{\frac{N}{2}-1}a_{2j+1}\omega_{\frac{N}{2}}^{ij} $$

$$ A_{i+\frac{N}{2}}=\sum_{j=0}^{\frac{N}{2}-1}a_{2j}\times\omega_N^{2(i+\frac{N}{2})j}+\sum_{j=0}^{\frac{N}{2}-1}a_{2j+1}\times\omega_N^{(i+\frac{N}{2})(2j+1)} $$

$$ =\sum_{j=0}^{\frac{N}{2}-1}a_{2j}\times\omega_{\frac{N}{2}}^{ij}+\omega_N^{i+\frac{N}{2}}\sum_{j=0}^{\frac{N}{2}-1}a_{2j+1}\omega_{\frac{N}{2}}^{ij} $$

$$ =\sum_{j=0}^{\frac{N}{2}-1}a_{2j}\times\omega_{\frac{N}{2}}^{ij}-\omega_N^i\sum_{j=0}^{\frac{N}{2}-1}a_{2j+1}\omega_{\frac{N}{2}}^{ij} $$

于是我们惊喜地发现，这可以分治做。

为什么需要奇偶分类呢？也许这就是傅里叶的伟大之处吧。

### 5.4 逆变换

那么FFT的逆变换怎么写呢？

令人惊讶的是，恰有

$$ a_i=\frac{\sum_{j=0}^{N-1}A_j\times\omega_N^{-ij}}{N} $$

**为什么它是正确的？**

$$ a_i=\frac{\sum_{j=0}^{N-1}A_j\times\omega_N^{-ij}}{N} $$

$$ =\frac{\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}a_k\times\omega_N^{jk}\times\omega_N^{-ij}}{N} $$

$$ =\frac{\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}a_k\times\omega_N^{j(k-i)}}{N} $$

$$ =\frac{\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}a_k[i==k]}{N} $$

$$ =\frac{N\times a_i}{N} $$

$$ =a_i $$

---

### 5.5 蝴蝶变换

我们将 $ N=8 $ 的分治情况手动模拟一下，可以得到:

| 区间大小 | $ id_0 $ | $ id_1 $ | $ id_2 $ | $ id_3 $ | $ id_4 $ | $ id_5 $ | $ id_6 $ | $ id_7 $ |
|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
| 8 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| 4 | 0 | 2 | 4 | 6 | 1 | 3 | 5 | 7 |
| 2 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |
| 1 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |

把这个表格用二进制描述

| 区间大小 | $ id_0 $ | $ id_1 $ | $ id_2 $ | $ id_3 $ | $ id_4 $ | $ id_5 $ | $ id_6 $ | $ id_7 $ |
|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
| 8 | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
| 4 | 000 | 010 | 100 | 110 | 001 | 011 | 101 | 111 |
| 2 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
| 1 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |

我们发现表格的第一行和最后一行二进制是反转的，这样我们就发现了 $ FFT $ 的非递归写法，代码如下:

``` cpp
#include<math.h>
#include<stdio.h>
#include<algorithm>
using namespace std;
const double pi=acos(-1.0);
int n,m,res=0,N=1,len,revers[2097160];
long long ans[2097160];
int i,j,k,l;
struct node{
  double x,y;
  node(double x=0,double y=0):x(x),y(y){}
  node operator*(const node &b){
    return node(x*b.x-y*b.y,x*b.y+y*b.x);
  }
  node operator+(const node &b){
    return node(x+b.x,y+b.y);
  }
  node operator-(const node &b){
    return node(x-b.x,y-b.y);
  }
}a[2097160],b[2097160],T,t,x,y;
void FFT(node *a,double flag){
  for(i=0;i<N;i++)
    if(i<revers[i])
      swap(a[i],a[revers[i]]);
  for(j=1;j<N;j<<=1){
    T=node(cos(pi/j),flag*sin(pi/j));
    for(k=0;k<N;k+=(j<<1)){
      t=node(1,0);
      for(l=0;l<j;l++,t=t*T){
        x=a[k+l],y=t*a[k+j+l];
        a[k+l]=x+y;
        a[k+j+l]=x-y;
      }
    }
  }
}
int main(){
  scanf("%d%d",&n,&m);
  n++;
  m++;
  for(i=0;i<n;i++)    
    scanf("%lf",&a[i].x);
  for(i=0;i<m;i++)
    scanf("%lf",&b[i].x);
  for(;N<max(n,m)<<1;N<<=1,len++);
  for(i=0;i<=N;i++)
    revers[i]=(revers[i>>1]>>1)|((i&1)<<(len-1));
  FFT(a,1);
  FFT(b,1);
  for(i=0;i<=N;i++)
    a[i]=a[i]*b[i];
  FFT(a,-1);
  for(i=0;i<=N;i++)
    ans[i]+=(long long)(a[i].x/N+0.5);
  for(;!ans[N]&&N;N--);
  N++;
  for(i=0;i<n+m-2;i++)
    printf("%lld ",ans[i]);
  printf("%lld\n",ans[n+m-2]);
  return 0;
}
```

### 5.6 精度问题

把上面的代码加入模操作后提交到 [洛谷P4245【模板】任意模数NTT](https://www.luogu.org/problemnew/show/P4245) 里去发现光荣`WA`。

然后就发现 $ FFT $ 有精度问题，那么如何避免呢？

#### 5.6.1 原根

如果 $ g $ 的 $ 0 $ ~ $ \phi(p)-1 $ 在模 $ p $ 意义下正好遍历了 $ 1 $ ~ $ p-1 $ 中与 $ p $ 互质的 $ \phi(p) $ 个数，那么称 $ g $ 为 $ p $ 的原根。

当 $ p $ 为质数时，我们发现如果用 $ g^\frac{p-1}{N} $ 代替单位复根（记为 $ g_N $ ）,它满足单位复根的所有性质:

1. 对于任意整数 $ n\geqslant0,k\geqslant0,d\geqslant0 $ 有

$$ g_{N}^{k}=g_{Nd}^{kd} $$

2. 对于任意整数 $ n\geqslant0,k\geqslant0 $ 有

$$ g_{N}^{k}=g_{2N}^{2k} $$

3. 对于任意整数 $ n\geqslant0 $ 有

$$ g_{N}^{\frac{N}{2}}=-1 $$

4. 对于任意整数 $ n\geqslant0,i\geqslant0,j\geqslant0 $ 有

$$ g_{N}^{i+j}=g_{N}^{i}\times g_{N}^{j} $$

5. 对于任意整数 $ n\geqslant0,i\geqslant0,j\geqslant0 $ 有

$$ g_{N}^{ij}=(g_{N}^{i})^j $$

#### 5.6.2 能选的质数

一般情况下有三个质数可选:

- $ 469762049=7\times 2^{26}+1 $

- $ 998244353=119\times 2^{23}+1 $

- $ 1004535809=749\times 2^{21}+1 $

当 $ p $ 取上面几个质数时， $ g=3 $ ， $ p-1 $ 中有很多 $ 2 $ 的因子， $ FFT $ 中 $ N $ 又都是 $ 2 $ 的次幂，所以上面三个质数一定要记下来。

代码如下:

``` cpp
#include<stdio.h>
#include<algorithm>
using namespace std;
const long long p=998244353,g=3,invg=332748118;
int n,m,res=0,N=1,len,revers[2097160];
long long ans[2097160],a[2097160],b[2097160],T,t,x,y;
int i,j,k,l;
inline long long power(long long a,long long k,long long p){
  long long ans=1,t=a;
  for(;k;k>>=1,t=t*t%p)
    if(k&1)
      ans=ans*t%p;
  return ans;
}
void NTT(long long *a,long long flag){
  for(i=0;i<N;i++)
    if(i<revers[i])
      swap(a[i],a[revers[i]]);
  for(j=1;j<N;j<<=1){
    T=power(flag==1?g:invg,(p-1)/j/2,p);
    for(k=0;k<N;k+=(j<<1)){
      t=1;
      for(l=0;l<j;l++,t=t*T%p){
        x=a[k+l],y=t*a[k+j+l]%p;
        a[k+l]=(x+y)%p;
        a[k+j+l]=((x-y)%p+p)%p;
      }
    }
  }
}
int main(){
  scanf("%d%d",&n,&m);
  n++;
  m++;
  for(i=0;i<n;i++){
    scanf("%lld",a+i);
    a[i]=a[i]%p;
  }
  for(i=0;i<m;i++){
    scanf("%lld",b+i);
    b[i]=b[i]%p;
  }
  for(;N<max(n,m)<<1;N<<=1,len++);
  for(i=0;i<=N;i++)
    revers[i]=(revers[i>>1]>>1)|((i&1)<<(len-1));
  NTT(a,1);
  NTT(b,1);
  for(i=0;i<=N;i++)
    a[i]=a[i]*b[i]%p;
  NTT(a,-1);
  for(i=0;i<=N;i++)
    ans[i]=a[i]*power(N,p-2,p)%p;
  for(i=0;i<n+m-2;i++)
    printf("%lld ",ans[i]);
  printf("%lld\n",ans[n+m-2]);
  return 0;
}
```

### 5.7 换个角度看$\bigoplus$卷积

我们再回忆一下`4.2节`的内容，我们在做 $ \bigoplus $ 卷积时，其实可以做 $ lg_N $ 遍 $ FFT $ ，然后又因为 $ \omega_2^0=1 $ ， $ \omega_2^1=-1 $ ，就可以得到`4.1节`的结论了。

### 5.8 FFT再拓展

在洛谷上搜索`FFT`会发现还有一道 [洛谷P4245【模板】任意模数NTT](https://www.luogu.org/problemnew/show/P4245) ，怎么lòng？
我们定三个质数做模数，一般是`5.6.2节`中的三个，分别进行 $ NTT $ ，最后用 $ CRT $ 强制求出解，注意爆`long long`（对，你没看错，不是爆`int`）， $ CRT $ 时模的勤劳一点。

贴代码:

``` cpp
#include<cstdio>
#include<string.h>
#include<algorithm>
using namespace std;
int mod;
inline int power(int a,int k,int p){
  int ans;
  for(ans=1;k;k>>=1,a=1LL*a*a%p)
    if(k&1)
      ans=1LL*ans*a%p;
  return ans;
}
const int p[3]={469762049,998244353,1004535809},g=3;
class longint{
  private:
    int A,B,C;
    static longint reduce(const longint &x){
      return longint(x.A+(x.A>>31&p[0]),x.B+(x.B>>31&p[1]),x.C+(x.C>>31&p[2]));
    }
  public:
    explicit longint(){}
    explicit longint(int x){
      A=B=C=x;
    }
    explicit longint(int x,int y,int z){
      A=x;
      B=y;
      C=z;
    }
    friend longint operator+(const longint &lsh,const longint &rsh){
      return reduce(longint(lsh.A+rsh.A-p[0],lsh.B+rsh.B-p[1],lsh.C+rsh.C-p[2]));
    }
    friend longint operator-(const longint &lsh,const longint &rsh){
      return reduce(longint(lsh.A-rsh.A,lsh.B-rsh.B,lsh.C-rsh.C));
    }
    friend longint operator*(const longint &lsh,const longint &rsh){
      return longint(1LL*lsh.A*rsh.A%p[0],1LL*lsh.B*rsh.B%p[1],1LL*lsh.C*rsh.C%p[2]);
    }
    int get(){
      long long x=1LL*(B-A+p[1])%p[1]*power(p[0],p[1]-2,p[1])%p[1]*p[0]+A;
      return (1LL*(C-x%p[2]+p[2])%p[2]*power(1LL*p[0]*p[1]%p[2],p[2]-2,p[2])%p[2]*(1LL*p[0]*p[1]%mod)+x)%mod;
    }
};
const int Array_Limit=1<<18;
int N,len,revers[Array_Limit|5];
longint exp[Array_Limit|5];
void init(int n){
  len=-1;
  N=1;
  while(N<n){
    N<<=1;
    len++;
  }
  for(int i=1;i<=N;i++)
    revers[i]=revers[i>>1]>>1|(i&1)<<len;
  longint t(power(g,(p[0]-1)/N,p[0]),power(g,(p[1]-1)/N,p[1]),power(g,(p[2]-1)/N,p[2]));
  exp[0]=longint(1);
  for(longint *i=exp;i!=exp+N;i++)
    *(i+1)=*i*t;
}
void NTT(longint *a,int flag){
  for(int i=1;i<Array_Limit;i++)
    if(i<revers[i])
      swap(a[i],a[revers[i]]);
  for(int i=1;i<N;i<<=1){
    int t=N/i>>1;
    for(int j=0;j<N;j+=i<<1)
      for(int k=0;k<i;k++){
        longint x=flag==1?exp[t*k]:exp[N-t*k],y=a[i+j+k]*x;
        x=a[j+k];
        a[j+k]=x+y;
        a[i+j+k]=x-y;
      }
  }
  if(flag==-1){
    longint inv_N(power(N,p[0]-2,p[0]),power(N,p[1]-2,p[1]),power(N,p[2]-2,p[2]));
    for(longint *i=a;i!=a+N;i++)
      *i=(*i)*inv_N;
  }
}
int n,m;
longint a[Array_Limit|5],b[Array_Limit|5];
int main(){
  scanf("%d%d%d",&n,&m,&mod);
  n++;
  m++;
  for(int i=0,x;i<n;i++){
    scanf("%d",&x);
    a[i]=longint(x%mod);
  }
  for(int i=0,x;i<m;i++){
    scanf("%d",&x);
    b[i]=longint(x%mod);
  }
  init(n+m);
  NTT(a,1);
  NTT(b,1);
  for(int i=0;i<N;i++)
    a[i]=a[i]*b[i];
  NTT(a,-1);
  for(int i=0;i<n+m-2;i++)
    printf("%d ",a[i].get());
  printf("%d\n",a[n+m-2].get());
  return 0;
}
```