浅析卷积

· · 个人记录

-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*k=i}{a_j\times b_k} ## 2. 当 $*$ 指 $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;
}
**为什么它是正确的？** 讨论 $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
  1. C_i=A_i\times B_i

  2. 用其逆变换将 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 卷积的变换就是按位前缀和。它的逆变换其实是一个脑筋急转弯，只要把循环倒过来，把 += 改成 -= 就可以了。代码如下：

#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;
}
## 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 具体实现

```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}
  1. 对于任意整数 n\geqslant0,k\geqslant0 有
\omega_{N}^{k}=\omega_{2N}^{2k}
  1. 对于任意整数 n\geqslant0 有
\omega_{N}^{\frac{N}{2}}=-1
  1. 对于任意整数 n\geqslant0,i\geqslant0,j\geqslant0 有
\omega_{N}^{i+j}=\omega_{N}^{i}\times\omega_{N}^{j}
  1. 对于任意整数 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 的非递归写法，代码如下:

#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 ~ \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}
  1. 对于任意整数 n\geqslant0,k\geqslant0 有
g_{N}^{k}=g_{2N}^{2k}
  1. 对于任意整数 n\geqslant0 有
g_{N}^{\frac{N}{2}}=-1
  1. 对于任意整数 n\geqslant0,i\geqslant0,j\geqslant0 有
g_{N}^{i+j}=g_{N}^{i}\times g_{N}^{j}
  1. 对于任意整数 n\geqslant0,i\geqslant0,j\geqslant0 有
g_{N}^{ij}=(g_{N}^{i})^j

5.6.2 能选的质数

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

当 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 换个角度看 \bigoplus 卷积

我们再回忆一下4.2节的内容，我们在做 \bigoplus 卷积时，其实可以做 lg_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只有初一，如今重整时可能会有一些没发现的错误，如有欢迎指正。