（原创重整）浅析卷积

-1.先赞后看，养成习惯

本文最早发布于2019年2月17日，经格式美化后于2021年8月6日同步重发在CSDN、洛谷博客及Gitee网页，转载请注明 出处

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\ast k=i}{a}_j\times b_k$$

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

2. 当 $\ast$ 指 $max/min$

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

#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{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. 当 $\ast$ 指 $\vee /\wedge$

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

3.1 用向量表示数

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

数值	向量表示	数值	向量表示
0	$\left\{\begin{array}{ccc}0& 0& 0\end{array}\right\}$	4	$\left\{\begin{array}{ccc}1& 0& 0\end{array}\right\}$
1	$\left\{\begin{array}{ccc}0& 0& 1\end{array}\right\}$	5	$\left\{\begin{array}{ccc}1& 0& 1\end{array}\right\}$
2	$\left\{\begin{array}{ccc}0& 1& 0\end{array}\right\}$	6	$\left\{\begin{array}{ccc}1& 1& 0\end{array}\right\}$
3	$\left\{\begin{array}{ccc}0& 1& 1\end{array}\right\}$	7	$\left\{\begin{array}{ccc}1& 1& 1\end{array}\right\}$

3.2 $\vee$的实质

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

$$3\vee 6=\left\{\begin{array}{ccc}0& 1& 1\end{array}\right\}\vee \left\{\begin{array}{ccc}1& 1& 0\end{array}\right\}=\left\{\begin{array}{ccc}1& 1& 1\end{array}\right\}=7$$

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

#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. 当 $\ast$ 指 $⨁$

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⨁6=\left\{\begin{array}{ccc}0& 1& 1\end{array}\right\}⨁\left\{\begin{array}{ccc}1& 1& 0\end{array}\right\}=\left\{\begin{array}{ccc}1& 0& 1\end{array}\right\}=5$$

再写一遍
$$3⨁6=\left\{\begin{array}{ccc}0& 1& 1\end{array}\right\}⨁\left\{\begin{array}{ccc}1& 1& 0\end{array}\right\}=\left\{\begin{array}{ccc}(0+1)\%2& (1+1)\%2& (1+0)\%2\end{array}\right\}=\left\{\begin{array}{ccc}1& 0& 1\end{array}\right\}=5$$

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

4.3 具体实现

$⨁$ 卷积的变换代码与 $\vee /\wedge$ 卷积的代码略有不同，但 $\vee /\wedge$ 卷积的代码也可以写成这种形式，具体题意见 洛谷P4717【模板】快速沃尔什变换 。

#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. 当 $\ast$ 指 $+$

公式恐惧症患者请果断按下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⩾0,k⩾0,d⩾0$ 有

$${\omega }_N^k={\omega }_{Nd}^{kd}$$

2. 对于任意整数 $n⩾0,k⩾0$ 有

$${\omega }_N^k={\omega }_{2N}^{2k}$$

3. 对于任意整数 $n⩾0$ 有

$${\omega }_N^{\frac{N}{2}}=-1$$

4. 对于任意整数 $n⩾0,i⩾0,j⩾0$ 有

$${\omega }_N^{i+j}={\omega }_N^i\times {\omega }_N^j$$

5. 对于任意整数 $n⩾0,i⩾0,j⩾0$ 有

$${\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$ 的分治情况手动模拟一下，可以得到:

区间大小	$i{d}_0$	$i{d}_1$	$i{d}_2$	$i{d}_3$	$i{d}_4$	$i{d}_5$	$i{d}_6$	$i{d}_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

把这个表格用二进制描述

区间大小	$i{d}_0$	$i{d}_1$	$i{d}_2$	$i{d}_3$	$i{d}_4$	$i{d}_5$	$i{d}_6$	$i{d}_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$ 的非递归写法，代码如下:

#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 里去发现光荣WA。

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

5.6.1 原根

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

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

1. 对于任意整数 $n⩾0,k⩾0,d⩾0$ 有

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

2. 对于任意整数 $n⩾0,k⩾0$ 有

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

3. 对于任意整数 $n⩾0$ 有

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

4. 对于任意整数 $n⩾0,i⩾0,j⩾0$ 有

$$g_N^{i+j}=g_N^i\times g_N^j$$

5. 对于任意整数 $n⩾0,i⩾0,j⩾0$ 有

$$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$ 的次幂，所以上面三个质数一定要记下来。

代码如下:

#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 换个角度看 $⨁$ 卷积

我们再回忆一下4.2节的内容，我们在做 $⨁$ 卷积时，其实可以做 $l{g}_N$ 遍 $FFT$ ，然后又因为 ${\omega }_2^0=1$，${\omega }_2^1=-1$，就可以得到4.1节的结论了。

5.8 FFT再拓展

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

贴代码:

#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;
}

6. 写在最后

初稿成型时Up只有初一，如今重整时可能会有一些没发现的错误，如有欢迎指正。