【参考博客】
0.简述
FFT全称：Fast Fourier Transformation
中文名： 快速傅里叶离散变换
作用：以 \(O(n \log n)\) 的复杂度计算多项式乘法
多项式乘法：
例：\((x^2+x+1)(2x^2+x-1)\)
习惯性的 可以 \(O(n^2)\) 求出 原式=\(2x^4+3x^3+2x^2-1\)
为了简洁，下面引入记号：
记 \(F(x)\) 表示一个多项式，可简称为 “多项式 \(F\) ”
若有一个 \(n\) 次多项式 \(F\)
\(F(x)=a_1+a_2x^1+...a_nx^n\)
其中 \(a_i\) 为参数，也叫系数
通常把 \(F(x)\) 的 \(k\) 次项系数称为 \(F[k]\)
即 \(F(x)=\sum_{i=0}^n F[i]x^i\)
卷积
如果要求 \(F(x)\) 与 \(G(x)\) 的乘积，如何写？
设 \(H(x)=F(x) * G(x)\)
那么 \(H[k]=\sum_{i=0}^k F[i] G[k-i]\)
等价于 \(H[k]=\sum_{i+j=k} F[i] G[j]\)
形如 \(H[k]=\sum_{i \otimes j=k} F[i] G[j]\) 的柿子称为卷积（这里点名批评一下卷·狗素耗·鸡）
那么你可能观察到了，多项式乘法就是加法卷积
1.DFT 与 IDFT
多项式的点值表达：
很简单，由算数基本定理可以得知：
如果知道 \(n+1\) 个点 \((x_0,y_0),(x_1,y_1)...(x_n,y_n)\)，即可确定一个 \(n\) 次多项式 \(F\)
那么，如果将数列 \(X=\{x_0,x_1,...,x_n\}\) 代入到 \(F(x)\) 中得到 \((x_0,y_0),(x_1,y_1)...(x_n,y_n)\)，代入到 \(G(x)\) 中得到 \((x_0,y_0'),(x_1,y_1')...(x_n,y_n')\)
那么如果 \(H(x)=F(x)*G(x)\)，即可得到将数列 \(X\) 带入到 \(H(x)\) 的点为 \((x_0,y_0 \times y_0'),(x_1,y_1 \times y_1')...(x_n,y_n \times y_n')\)
我们可以发现，只要知道了 \(F(x)\) 与 \(G(x)\) 的点值表达，即可知道 \(H(x)\) 的点值表达，进而可以确定唯一的多项式 \(H(x)\)
但是我们发现 \(H(x)\) 是 \(2n\) 次多项式，而 \(F(x)\) 和 \(G(x)\) 是 \(n\) 次多项式，这之中少了 \(n\) 个点，怎么办呢？
再找 \(n\) 个不就行了，这时只需要做 \(2n+1\) 次乘法即可。
那么我们可以把系数表达转换为点值表 -> 点值表达相乘 -> 把点值表达转换为系数表达。
这就是 FFT 的算法流程。
“把系数表达转换为点值表达” 的算法叫做 DFT
“把点值表达转换为系数表达” 的算法叫做 IDFT（DFT 的逆运算）
2.单位根
DFT 代入什么由你自己决定,只要点值个数足够
事实证明,找一些奇奇怪怪的东西代入进去是个好想法（）
上古之时,有一位达捞横空出世,他就是傅里叶！
他把单位根 \(\omega^0_{n+1}...\omega^n_{n+1}\) 代入了多项式
然后 \(\omega^0_{n+1}...\omega^n_{n+1}\) 有一些神奇性质
然后分治，就得到了 \(O(n \log n)\) 的 FFT
复数
详见高中数学必修二（好像是）
这里重点说一下 复数相乘：
模长相同（勾股定理即可）
辐角相加（这个得证相似很麻烦）
接下来说说 单位根：
\(n\) 次单位根（\(n\) 为正整数）是 \(n\) 次幂为 \(1\) 的复数。
即为 \(x^n=1\) 的复数解
圆心在原点且半径为 \(1\) 的圆叫 单位圆：
在复平面上，单位圆上的点模长都为 \(1\)
由于 \(n\) 次单位根是 \(n\) 次幂为 \(1\) 的复数
考虑一个复数 \(x\)：
若 \(|x|>1\) 则 \(|x^n|=|x|^n>1\)
若 \(|x|=1\) 则 \(|x^n|=|x|^n=1\)
若 \(|x|<1\) 则 \(|x^n|=|x|^n<1\)
所以只有模长等于 \(1\) 的复数才 有可能 成为 \(n\) 次单位根
容易找到 : 幅角为 \(0,\frac{1}{n},...\frac{n-1}{n}\) 圆周 的复数,都是单位根,共 \(n\) 个
也就是说,假如一个复数,其幅角为 \(\frac{a}{n} (0≤a 由复数相乘，辐角相加可知这玩意的 \(n\) 次方幅角是 \(a\) 倍圆周 所以 \(n\) 次单位根 \(n\) 等分单位圆 我们从 \(1\) 开始（\(1\) 一定是单位根）,沿着单位圆 逆时针 把单位根标上号 注： 虽然我们只承认 \(\omega^0_n,\omega^1_n...\omega^{n-1}_n\) 但也有 \(\omega^k_n,k\in(-\infty,0)\cup[n,\infty)\) 的情况 而且有 \(\omega^k_n=\omega^{k\%n}_n\) 单位根性质： \(\omega^k_n=(\omega^1_n)^k\) 由辐角相加可以证明 \(\omega^j_n\times \omega^k_n=\omega^{j+k}_n\) 也是可以从辐角相加的角度理解 \(\omega^{kp}_{np}=\omega^k_n\) 也也是可以从辐角相加的角度理解 若 \(n\) 是偶数，\(\omega^{k+n/2}_n=-\omega^k_n\) 乘以 \(\omega^{n/2}_n\) 相当于旋转 \(180\) 度，也就是关于原点对称，也就是取相反数 3.DFT 加速 我们来讲讲 FFT 的分治方式 我们按然后按指数奇偶把多项式分成两部分 \(F(x)=(F[0]+F[2]x^2+...+F[n-2]x^{n-2})+(F[1]+F[3]x^2+...+F[n-1]x^{n-1})\) 这里保证 \(n\) 是 \(2\) 的整幂，不会出现分得不均匀的情况，如果不满足，可以在高位上补 \(0\) 又设两个有 \(n/2\) 项的多项式 \(FL(x)\) 和 \(FR(x)\) \(FL(x)=F[0]+F[2]x+...+F[n-2]x^{n/2-1}\) \(FL(x)=F[1]+F[3]x+...+F[n-1]x^{n/2-1}\) 则可以得出 \(F(x)=FL(x^2)+xFR(x^2)\) 我们把 \(\omega_n^k\) 代入 \(F(x)\) \(k \[\begin{aligned} F(\omega_n^k)&=FL((\omega_n^k)^2)+\omega_n^k FR((\omega_n^k)^2)\\ &=FL(\omega_{n/2}^k)+\omega_n^k FR(\omega_{n/2}^k) \end{aligned} \] \(k \[\begin{aligned} F(\omega_n^k)&=FL((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2} FR((\omega_n^{k+n/2})^2)\\ &=FL(\omega_{n}^{2k+n})+\omega_n^{k+n/2}FR(\omega_{n}^{2k+n})\\ &=FL(\omega_{n}^{2k})+\omega_n^{k+n/2}FR(\omega_{n}^{2k})\\ &=FL(\omega_{n}^{2k})-\omega_n^{k}FR(\omega_{n}^{2k})\\ &=FL(\omega_{n/2}^{k})-\omega_n^{k}FR(\omega_{n/2}^{k}) \end{aligned} \] 比对一下两个式子，发现只是正负号的区别 这里如果我们知道 \(FL(x)\) 和 \(FR(x)\) 在 \(\omega_{n/2}^{0},\omega_{n/2}^{1},...\omega_{n/2}^{n/2-1}\) 的点值表达，就能 \(O(n)\) 求出 \(F(x)\) 在 \(\omega_{n}^{0},\omega_{n}^{1},..\omega_{n}^{-1}\) 的点值表达 问题在于我们不知道 \(FL(x)\) 和 \(FR(x)\) 的点值表达，怎么办？ 暴力代入？ 并不是，上面的操作，可以看做一个分治过程 这个可以一直 分治 下去,直到多项式只剩下一个项为止 4.DFT 代码实现 还有一个小细节 上文有一句话：“保证 \(n\) 是 \(2\) 的整幂，不会出现分得不均匀的情况” 实际上，\(n\) 不一定是 \(2\) 的正整数次幂 我们可以补项，在最高次添加一些 系数 为 \(0\) 的项，不会影响计算结果 讲完了这些我们可以开始写 DFT 了 1.复数结构体 根据复数的四则运算重载 CP是 \(complex\) 的简称（直接用会 CE，别用 STL,会被卡常） struct complex struct CP { double a,b;//a+bi CP (double A=0,double B=0){a=A,b=B;} CP operator+(CP const &B) const{return CP(a+B.a,b+B.b);} CP operator-(CP const &B) const{return CP(a-B.a,b-B.b);} CP operator*(CP const &B) const{return CP(a*B.a-b*B.b,a*B.b+b*B.a);} CP operator/(CP const &B) const { double t=B.a*B.a+B.b*B.b; return CP((a*B.a+b*B.b)/t,(b*B.a-a*B.b)/t); } CP operator+=(CP const &B){return (*this)=(*this)+B;} CP operator-=(CP const &B){return (*this)=(*this)-B;} CP operator*=(CP const &B){return (*this)=(*this)*B;} CP operator/=(CP const &B){return (*this)=(*this)/B;} }; 2.预处理单位根 直接欧拉公式 可以考虑弧度制（详见高中数学必修二） 可得 \(\omega_{n}^{1}\) 坐标为 \({\Big(} \cos(\frac{2 \pi}{n}),\sin(\frac{2 \pi}{n}) {\Big)}\) 只要求出 \(\omega_n^1\) 并自乘 \(k\) 次可得到 \(\omega_m^k\) 3.简单实现 DFT 这里贴一个大佬的 （其实是我不想写了） DFT #include #include #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} //除法没用 }f[Maxn<<1],sav[Maxn<<1]; void dft(CP *f,int len) { if (len==1)return ;//边界 //指针的使用比较巧妙 CP *fl=f,*fr=f+len/2; for (int k=0;k for (int k=0;k {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];} dft(fl,len/2); dft(fr,len/2);//处理子问题 //由于每次使用的单位根次数不同(len次单位根),所以要重新求。 CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0); for (int k=0;k //这里buf = (len次单位根的第k个) sav[k]=fl[k]+buf*fr[k];//(1) sav[k+len/2]=fl[k]-buf*fr[k];//(2) //这两条语句具体见Tips的式子 buf=buf*tG;//得到下一个单位根。 }for (int k=0;k } int main() { scanf("%d",&n); for (int i=0;i //一开始都是实数,虚部为0 for(m=1;m //把长度补到2的幂,不必担心高次项的系数,因为默认为0 dft(f,m); for(int i=0;i printf("(%.4f,%.4f)\n",f[i].x,f[i].y); return 0; } 5.IDFT理论与FFT初步实现 首先说一个 结论 : 把 DFT 中的 \(\omega_n^1\) 换成 \(\omega_n^-1\)，做完之后除以 \(n\) 即可 从这里可以拓展出 单位根反演，但是我不会…… 还是多项式 \(F(x)\)，设我们变换之后，得到的点值数列为 \(G\) 即 \(G={\text {DFT}}(F)\) 可得： \[G_k=\sum_{i=0}^{n-1}(\omega_n^k)^iF[i] \] 那么结论就是： \[n\times F[k]=\sum_{i=0}^{n-1}(\omega_n^{-k})^iG_i \] 通过这个柿子可以把点值还原 证明 \[\begin{aligned} 右边&=\sum_{i=0}^{n-1}(\omega_n^{-k})^iG_i\\ &代入可得\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\omega_n^{ij}\omega_n^{-ik}F[j]\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\omega_n^{i(j-k)}F[j]\\ \end{aligned} \] 分讨： \(j=k\) 贡献为： \[\sum_{i=0}^{n-1}\omega^0_nF[k]=nF[k] \] \(j \not =k\) 设 \(p=j-k\) 贡献为： \[\sum_{i=0}^{n-1}\omega^{ip}_nF[k+p]=\omega_n^p \sum_{i=0}^{n-1}\omega^{i}_nF[k+p]=\omega_n^p(\frac{\omega_n^{0}-1}{\omega_n^p-1})F[k+p]=0 \] 证毕 这个东西本质是单位根反演（） 或者可以理解成范德蒙德矩阵求逆…… 蒟蒻狂悲 \[n\times F[k]=\sum_{i=0}^{n-1}(\omega_n^{-k})^iG_i \] 相当于把 \(G\) 数列当做系数,再代入一遍求值 不同的是，这次代入的是 \(\omega^0_n,\omega^{-1}_n...\omega^{-n+1}_n\) 如何做？ 求出 \(\omega^{-1}_n\) 即可：\({\Big (} \cos(\frac{2 \pi}{n}),-\sin(\frac{2 \pi}{n}) {\Big)}\) 至此我们已经写出了第一个版本的 FFT 代码很好写，和 DFT 不同的地方很少 前文我们说过，IDFT 和 DFT 只有一个符号的区别 那么我们何不减少一下代码量呢：（大佬的） FFT #include #include #define Maxn 1350000 using namespace std; const double Pi=acos(-1); int n,m; struct CP { CP (double xx=0,double yy=0){x=xx,y=yy;} double x,y; CP operator + (CP const &B) const {return CP(x+B.x,y+B.y);} CP operator - (CP const &B) const {return CP(x-B.x,y-B.y);} CP operator * (CP const &B) const {return CP(x*B.x-y*B.y,x*B.y+y*B.x);} //除法没用 }f[Maxn<<1],p[Maxn<<1],sav[Maxn<<1]; //flag=1 -> DFT flag=0 -> IDFT void fft(CP *f,int len,bool flag) { if (len==1)return ; CP *fl=f,*fr=f+len/2; for (int k=0;k for (int k=0;k {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];} fft(fl,len/2,flag); fft(fr,len/2,flag); CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0); if (!flag)tG.y*=-1; for (int k=0;k sav[k]=fl[k]+buf*fr[k]; sav[k+len/2]=fl[k]-buf*fr[k]; buf=buf*tG; }for (int k=0;k } int main() { scanf("%d%d",&n,&m); for (int i=0;i<=n;i++)scanf("%lf",&f[i].x); for (int i=0;i<=m;i++)scanf("%lf",&p[i].x); for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂 fft(f,n,1);fft(p,n,1);//DFT for(int i=0;i fft(f,n,0); for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49)); return 0; } 6. FFT的精细实现 “分开” 的优化 具体的来说： 原来的递归版（数组下标，先偶后奇,从 0 开始）： 第 1 层 0 1 2 3 4 5 6 7 第 2 层 0 2 4 6|1 3 5 7 第 3 层 0 4|2 6|1 5|3 7 第 4 层 0|4|2|6|1|5|3|7 我们要求出第 \(4\) 层的数组状况,然后才能往上合并， 这个东西貌似叫做“蝴蝶变换”？ 具体怎么求呢？ 注意到最后的序列是原序列的 二进制反转 如何得到二进制翻转后的数列呢？ 可以用递推 \(O(n)\) 搞定 待写 “合并” 的优化 观察合并的代码片段 待写 注意到指针fl[k] <--> f[k] . . . fr[k] <--> f[k+len/2] 那么我们完全可以用如下代码来代替： 待写 这就规避了所有的数组 \(\text{Copy}\) 可以写出如下代码： 待写 三角函数优化 三角函数の求值很慢，我们却使用了 \(O(n\log n)\) 次三角函数求值，可以考虑预处理来减小常数 不过一种更有效的方式是：迭代实现 Code #include #define int long long #define ll long long #define fd(i,a,b) for(int i=(a);i<=(b);i=-~i) #define bd(i,a,b) for(int i=(a);i>=(b);i=~-i) #define db(x) cout<<"DEBUG "<<#x<<" = "< using namespace std; const int N=1e6+3e5+509,M=1e6+509,mod=998244353; const double Pi=acos(-1); struct CP { double a,b;//a+bi CP (double A=0,double B=0){a=A,b=B;} CP operator+(CP const &B) const{return CP(a+B.a,b+B.b);} CP operator-(CP const &B) const{return CP(a-B.a,b-B.b);} CP operator*(CP const &B) const{return CP(a*B.a-b*B.b,a*B.b+b*B.a);} CP operator/(CP const &B) const { double t=B.a*B.a+B.b*B.b; return CP((a*B.a+b*B.b)/t,(b*B.a-a*B.b)/t); } CP operator+=(CP const &B){return (*this)=(*this)+B;} CP operator-=(CP const &B){return (*this)=(*this)-B;} CP operator*=(CP const &B){return (*this)=(*this)*B;} CP operator/=(CP const &B){return (*this)=(*this)/B;} }f[N<<1],g[N<<1]; int p[N<<1]; int n,m,ans[N<<1]; void FFT(CP *f,bool flag) { fd(i,0,n-1) if(i for(int i=2;i<=n;i<<=1) { int len=i>>1; CP tg(cos(2*Pi/i),sin(2*Pi/i)); if(!flag) tg.b*=-1; for(int j=0;j { CP buf(1,0); for(int k=j;k { CP x=buf*f[len+k]; f[len+k]=f[k]-x; f[k]+=x,buf*=tg; } } } } inline void Print(CP *f,bool flag=0)//flag=1 十进位 =0 不进位 { if(flag) { fd(i,0,n-1) { ans[i]+=(int)(f[i].a/n+0.49); ans[i+1]+=ans[i]/10,ans[i]%=10; } while(n>-1&&!ans[n]) --n; if(n==-1) cout<<"0\n"; else bd(i,n,0) cout< } else { fd(i,0,m) cout<<(int)(f[i].a/n+0.49)<<' '; } } char s1[N],s2[N]; inline void Init(bool flag=0) { if(flag) { scanf("%s%s",s1,s2); n=strlen(s1)-1;m=strlen(s2)-1; fd(i,0,n) f[i].a=(double)(s1[n-i]-'0'); fd(i,0,m) g[i].a=(double)(s2[m-i]-'0'); } else { cin>>n>>m; fd(i,0,n) cin>>f[i].a; fd(i,0,m) cin>>g[i].a; } } signed main() { // #define FJ #ifdef FJ freopen(".in","r",stdin); freopen(".out","w",stdout); #else // freopen("A.in","r",stdin); // freopen("A.out","w",stdout); #endif // ios::sync_with_stdio(0); // cin.tie(0); cout.tie(0); Init(1); for(m+=n,n=1;n<=m;n<<=1); fd(i,0,n-1) p[i]=(p[i>>1]>>1)|((i&1)?n>>1:0); FFT(f,1);FFT(g,1); fd(i,0,n) f[i]*=g[i]; FFT(f,0); Print(f,1); return 0; } 这个版本是最经典的 建议背诵全文 “三次变两次” 优化 根据 \((a+bi)(c+di)=ac-bd+(ad+bc) \times i\) 假设我们要求 \(F(x) \times G(x)\) 则有 \(H(x)=F(x)+i \times G(x)\) 则 \(H(x)^2=F(X)^2-G(x)^2+2i \times F(x)G(x)\) 发现 \(H(x)^2\) 的虚部为 \(2i \times F(x)G(x)\) 也就是说求出 \(H(x)^2\) 后把虚部除 \(2\) 即可 Code #include #define int unsigned long long #define ll long long #define fd(i,a,b) for(int i=(a);i<=(b);i=-~i) #define bd(i,a,b) for(int i=(a);i>=(b);i=~-i) #define db(x) cout<<"DEBUG "<<#x<<" = "< using namespace std; const int N=1e6+3e5+509,M=1e6+509,mod=998244353; const double Pi=acos(-1); struct CP { double a,b;//a+bi CP (double A=0,double B=0){a=A,b=B;} CP operator+(CP const &B) const{return CP(a+B.a,b+B.b);} CP operator-(CP const &B) const{return CP(a-B.a,b-B.b);} CP operator*(CP const &B) const{return CP(a*B.a-b*B.b,a*B.b+b*B.a);} CP operator/(CP const &B) const { double t=B.a*B.a+B.b*B.b; return CP((a*B.a+b*B.b)/t,(b*B.a-a*B.b)/t); } CP operator+=(CP const &B){return (*this)=(*this)+B;} CP operator-=(CP const &B){return (*this)=(*this)-B;} CP operator*=(CP const &B){return (*this)=(*this)*B;} CP operator/=(CP const &B){return (*this)=(*this)/B;} }f[N<<1]; int p[N<<1]; int n,m; void FFT(CP *f,bool flag) { fd(i,0,n-1) if(i for(int i=2;i<=n;i<<=1) { int len=i>>1; CP tg(cos(2*Pi/i),sin(2*Pi/i)); if(!flag) tg.b*=-1; for(int j=0;j { CP buf(1,0); for(int k=j;k { CP x=buf*f[len+k]; f[len+k]=f[k]-x; f[k]+=x,buf*=tg; } } } } void Print(CP *f,int len) { fd(i,0,m) cout<<(int)(f[i].b/n/2+0.49)<<' '; } signed main() { // #define FJ #ifdef FJ freopen(".in","r",stdin); freopen(".out","w",stdout); #else // freopen("A.in","r",stdin); // freopen("A.out","w",stdout); #endif ios::sync_with_stdio(0); cin.tie(0); cout.tie(0); cin>>n>>m; fd(i,0,n) cin>>f[i].a; fd(i,0,m) cin>>f[i].b; for(m+=n,n=1;n<=m;n<<=1); fd(i,0,n-1) p[i]=(p[i>>1]>>1)|((i&1)?n>>1:0); FFT(f,1); fd(i,0,n-1) f[i]*=f[i]; FFT(f,0); Print(f,m); return 0; } 7.后话 再见啦~~ 欲知后事如何，好像没有下文分解了…… 先放个 NTT 板子 NTT #include #define int long long #define ll long long #define fd(i,a,b) for(int i=(a);i<=(b);i=-~i) #define bd(i,a,b) for(int i=(a);i>=(b);i=~-i) #define db(x) cout<<"DEBUG "<<#x<<" = "< using namespace std; const int N=1e6+5e5+509,M=1e6+509,mod=998244353,G=3; int f[N<<1],g[N<<1]; int p[N<<1]; int n,m,ans[N<<1]; inline int qpow(int x,int y=mod-2) { x%=mod;int re=1; while(y) { if(y&1) (re*=x)%=mod; (x*=x)%=mod,y>>=1; } return re; } const int invG=qpow(G); inline int add(int x,int y){return (x+y inline int del(int x,int y){return (x-y<0)?x-y+mod:x-y;} inline int mul(int x,int y){return (ll)x*y%mod;} void NTT(int *f,bool flag) { fd(i,0,n-1) if(i for(int i=2;i<=n;i<<=1) { int len=i>>1; int tg=qpow(flag?G:invG,(mod-1)/i); for(int j=0;j { int buf(1); for(int k=j;k { int x=mul(buf,f[len+k]); f[len+k]=del(f[k],x), f[k]=add(f[k],x), buf=mul(buf,tg); } } } } inline void Print(int *f,bool flag=0)//flag=1 十进位 =0 不进位 { int invn=qpow(n); if(flag) { fd(i,0,n-1) { ans[i]+=mul(f[i],invn); ans[i+1]+=ans[i]/10; ans[i]%=10; } while(n>-1&&!ans[n]) --n; if(n==-1) cout<<"0\n"; else bd(i,n,0) cout< } else fd(i,0,m) cout< } char s1[N],s2[N]; inline void Init(bool flag=0) { if(flag) { scanf("%s%s",s1,s2); n=strlen(s1)-1;m=strlen(s2)-1; fd(i,0,n) f[i]=(int)(s1[n-i]-'0'); fd(i,0,m) g[i]=(int)(s2[m-i]-'0'); } else { cin>>n>>m; fd(i,0,n) cin>>f[i]; fd(i,0,m) cin>>g[i]; } } signed main() { // #define FJ #ifdef FJ freopen(".in","r",stdin); freopen(".out","w",stdout); #else // freopen("A.in","r",stdin); // freopen("A.out","w",stdout); #endif // ios::sync_with_stdio(0); // cin.tie(0); cout.tie(0); Init(0); for(m+=n,n=1;n<=m;n<<=1); fd(i,0,n-1) p[i]=(p[i>>1]>>1)|((i&1)?n>>1:0); NTT(f,1);NTT(g,1); fd(i,0,n-1) f[i]=mul(f[i],g[i]); NTT(f,0); Print(f,0); return 0; }