【真·通俗易懂】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\)

\[\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<n/2\),代入 \(\omega_n^{k+n/2}\)

\[\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 <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)\)

可得:

\[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 <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;
}
posted @   whrwlx  阅读(296)  评论(2编辑  收藏  举报
相关博文:
阅读排行:
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· 实操Deepseek接入个人知识库
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· 【.NET】调用本地 Deepseek 模型
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
点击右上角即可分享
微信分享提示