【真·通俗易懂】FFT 学习笔记
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∈\mathbb Z)\),那么这就是一个单位根
由复数相乘,辐角相加可知这玩意的 \(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<n/2\),代入 \(\omega_n^k\)
- \(k<n/2\),代入 \(\omega_n^{k+n/2}\)
比对一下两个式子,发现只是正负号的区别
这里如果我们知道 \(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 <cstdio>
#include <cmath>
#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<len;k++)sav[k]=f[k];
for (int k=0;k<len/2;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<len/2;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<len;k++)f[k]=sav[k];
}
int main()
{
scanf("%d",&n);
for (int i=0;i<n;i++)scanf("%lf",&f[i].x);
//一开始都是实数,虚部为0
for(m=1;m<n;m<<=1);
//把长度补到2的幂,不必担心高次项的系数,因为默认为0
dft(f,m);
for(int i=0;i<m;++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)\)
可得:
那么结论就是:
通过这个柿子可以把点值还原
证明
分讨:
-
\(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 \]
证毕
这个东西本质是单位根反演()
或者可以理解成范德蒙德矩阵求逆……
蒟蒻狂悲
相当于把 \(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 <cstdio>
#include <cmath>
#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<len;k++)sav[k]=f[k];
for (int k=0;k<len/2;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<len/2;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<len;k++)f[k]=sav[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<n;++i)f[i]=f[i]*p[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<bits/stdc++.h>
#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<<" = "<<x<<endl;
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<p[i]) swap(f[i],f[p[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<n;j+=i)
{
CP buf(1,0);
for(int k=j;k<j+len;k=-~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<<ans[i];
}
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<bits/stdc++.h>
#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<<" = "<<x<<endl;
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<p[i]) swap(f[i],f[p[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<n;j+=i)
{
CP buf(1,0);
for(int k=j;k<j+len;k=-~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<bits/stdc++.h>
#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<<" = "<<x<<endl;
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<mod)?x+y:x+y-mod;}
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<p[i]) swap(f[i],f[p[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<n;j+=i)
{
int buf(1);
for(int k=j;k<j+len;k=-~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<<ans[i];
}
else fd(i,0,m) cout<<mul(f[i],invn)<<' ';
}
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;
}
本文来自博客园,作者:whrwlx,转载请注明原文链接:https://www.cnblogs.com/whrwlx/p/18459333
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· 实操Deepseek接入个人知识库
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· 【.NET】调用本地 Deepseek 模型
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库