2023.2.9【数学】快速傅里叶变换

2023.2.9【模板】快速傅里叶变换（FFT）

好多天没写博客了qwq

题目描述

给定一个 n 次多项式 F(x)，和一个 m 次多项式 G(x)。

请求出 F(x) 和 G(x) 的卷积。

朴素（正常）思路

枚举计算的每一位，交叉相乘加起来计算答案，时间复杂度O($n^2$)

原地爆炸

这个时候就需要用到NOIP不考的FFT了

前置芝士

有关复数详见lxw (tqqqqqqqqqqqqqqqqqql%%%%%%%%%%%%%)巨佬的博客

复数 - Ricky2007 - 博客园 (cnblogs.com)

*泰勒展开

泰勒中值定理：若$f(x)$在$x_0$处有$n$阶导数，那么存在$x_0$邻域中的$x$,有

$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{f^{″}(x_0)}{2}{(x-x_0)}^2+...+\frac{f^{(n)}(x_0)}{n!}{(x-x_0)}^n+o({(x-x_0)}^n)$

泰勒展开经常用于将一个难以计算的函数逼近（或等效）于一个多项式，为了消掉$(x-x_0)$

，我们经常将函数在0处展开，即$x_0=0$

*欧拉公式

$e^{ix}=isinx+cosx$

证明：将三个柿子泰勒展开可以得到：

$e^x=1+x+\frac{1}{2!}{x}^2+\frac{1}{3!}{x}^3+...+\frac{1}{n!}{x}^n$

$sinx=x-\frac{1}{3!}{x}^3+\frac{1}{5!}{x}^5-\frac{1}{7!}{x}^7+...$

$cosx=1-\frac{1}{2!}{x}^2+\frac{1}{4!}{x}^4-\frac{1}{6!}{x}^6+...$

我们（神奇地）发现：$sinx+cosx=1+x-\frac{1}{2!}{x}^2-\frac{1}{3!}{x}^3+\frac{1}{4!}{x}^4+\frac{1}{5!}{x}^5...$

这难道不和$e^x$长得很像吗？

为了将$e^x$中一定项转为负号，我们引入$i=\sqrt{-}1$ :

$e^{ix}=1+ix-\frac{1}{2!}x-\frac{1}{3!}i{x}^3+....$

$isin$ $x+cosx=1+ix-\frac{1}{2!}x-\frac{1}{3!}i{x}^3...$

$e^{ix}=isin$ $x+cosx$

特殊地，代入$x=\pi$得：$e^{i\pi }=isin\pi +cos\pi =-1$

这不就是那个被很多人当作壁纸的公式吗，我们来思考它的实际作用

在一个笛卡尔坐标系中，一个从源点指向$(x,y)$的向量，可以使用与$x$轴夹角$\theta$和长度$\rho$来表达（极坐标），即$(\theta ,\rho$)

我们将Y坐标指定为复数域，X坐标指定为实数域，发现x + yi可以表示一个点

($\rho isin\theta +\rho cos\theta$)

这样以来，再使用欧拉公式转化，就得到了：$\rho e^{i\theta }$

一个单项就可以表达一个坐标，省去许多计算麻烦

*离散傅里叶级数

对于序列$<c_0,c_1,c_2....c_{n-1}>$ ，定义其$k$次离散傅里叶级数（DF）为：

$h({\omega }^k)=c_0+c_1{\omega }^k+c_2{\omega }^{2k}+...+c_{n-1}{\omega }^{(n-1)k}$

其中$\omega =e^{\frac{2i\pi }{n}}$

可以看作在极坐标内，一个单位角度是$\frac{2\pi }{n}$，每次走$k$单位的角度，一圈是$n$单位角度

（B站神犇up@3Blue1Brown的示意图）

*离散傅里叶变换 DFT

对于序列$<c_0,c_1,c_2,....c_{n-1}>$，构造一个新的序列：

$<h({\omega }^0),h({\omega }^1),h({\omega }^2),...,h({\omega }^{n-1})>$

这个过程叫做离散傅里叶变换，反之，我们称对于一个构造好的序列，求原序列的过程叫做的逆傅里叶变换（IFT）

*定理：一个序列$<c_0,c_1,...,c_{n-1}>$的傅里叶变换后的新序列$<h({\omega }^0),h({\omega }^1),...,h({\omega }^{n-1})>$其$(-k)$次DF值恰为原序列$c_k$的$n$倍

证明：设$h({\omega }^k)={\mathrm{\Sigma }}_{x=0}^{n-1}{\omega }^{kx}{c}_x,g({\omega }^{-k})={\mathrm{\Sigma }}_{y=0}^{n-1}{\omega }^{-ky}h({\omega }^y)$

代入$h({\omega }^y)$得：$g({\omega }^{-k})={\mathrm{\Sigma }}_{y=0}^{n-1}{\omega }^{-ky}{\mathrm{\Sigma }}_{x=0}^{n-1}{\omega }^{xy}{c}_x$

提出求和符号：$g({\omega }^{-k})={\mathrm{\Sigma }}_{y=0}^{n-1}{\mathrm{\Sigma }}_{x=0}^{n-1}{\omega }^{(x-k)y}{c}_x$

$={\mathrm{\Sigma }}_{x=0}^{n-1}{c}_x{\mathrm{\Sigma }}_{y=0}^{n-1}{\omega }^{(x-k)y}$

Key1. $x=k$时：

${\mathrm{\Sigma }}_{y=0}^{n-1}{\omega }^{(x-k)y}={\mathrm{\Sigma }}_{y=0}^{n-1}1=n$

Key2.$x\ne k$时：

${\mathrm{\Sigma }}_{y=0}^{n-1}[{\omega }^{(x-k)}{]}^y=[{\omega }^{(x-k)}{]}^0+[{\omega }^{(x-k)}{]}^1+...+[{\omega }^{(x-k)}{]}^{n-1}=\frac{1-[{\omega }^{(x-k)}{]}^n}{1-{\omega }^{(x-k)}}$

为什么$x^0+x^1+x^2+...+x^{n-1}=\frac{1-x^n}{1-x}$？

设其为$f(x)$，将上式乘$x$,可以得到

$x^1+x^2+...+x^n=f(x)-1+x^n=xf(x)$

$(x-1)f(x)=x^n-1,f(x)=\frac{x^n-1}{x-1}$

$\omega =e^{\frac{2i\pi }{n}}\to [{\omega }^{(x-k)}{]}^n=({\omega }^n{)}^{x-k}=1^{x-k}=1$

${\mathrm{\Sigma }}_{y=0}^{n-1}[{\omega }^{x-k}{]}^y=\frac{1-1}{1-{\omega }^{x-k}}=0$，不计入答案

$g({\omega }^{-k})={\mathrm{\Sigma }}_{x=0}^{n-1}{c}_x\ast [(x=k)?n:0]=n{c}_k$

$QED$

通过*****,我们就可以构造出一种计算多项式乘法的新方法：

1.将F(x)和G(x)进行傅里叶变换

2.将两者点乘($x\in [0,n-1]$上的值分别相乘)

3.将所得序列进行$(-k)$次方逆傅里叶变换，结果除以$n$

$O(n^2)$？

怎样加快计算过程？

考虑多项式$h(x)=c_0+c_{1}x+...+c_{n-1}{x}^{n-1}$

提出奇数$odd(x)=c_1+c_{3}x+c_5{x}^2+...+c_{n-1}{x}^{\frac{n}{2}-1}$

偶数$even(x)=c_0+c_{2}x+c_4{x}^2+...+c_{n-2}{x}^{\frac{n}{2}-1}$

可推得合并式$h(x)=even(x^2)+x\ast odd(x^2)$

代入：$h({\omega }^k)=even({\omega }^{2k})+{\omega }^k\ast odd({\omega }^{2k})$

可以发现：$h({\omega }^{k+\frac{n}{2}})=even({\omega }^{2k+n})+{\omega }^{k+\frac{n}{2}}\ast odd({\omega }^{2k+n})$

$=even({\omega }^{2k}\ast {\omega }^n)+{\omega }^{k+\frac{n}{2}}\ast odd({\omega }^{2k}\ast {\omega }^n)$

$\omega =e^{\frac{2i\pi }{n}}\to {\omega }^n=1,{\omega }^{\frac{n}{2}}=-1$

$h({\omega }^{k+\frac{n}{2}})=even({\omega }^{2k})-{\omega }^k\ast odd({\omega }^{2k})$

$h({\omega }^k)=even({\omega }^{2k})+{\omega }^k\ast odd({\omega }^{2k})$

其中通过${\omega }^{2k}$，我们将要计算的$h({\omega }^k)$和$h({\omega }^{k+\frac{n}{2}})$共计n项，转化成$even$和$odd$共计$n$项，每次将序列拆成奇项和偶项两部分，共拆$lo{g}_{2}n$层，共计时间复杂度$O(n{\log }_{2}n)$

$————FFT$快速傅里叶变换

合并过程图解：例：$n=8$

小技巧：怎样快速将序列拆成奇项和偶项？

考虑到我们要将偶项也拆成其中的奇和偶、再拆为奇和偶...我们可以发现，每次拆项都是按照数字二进制的第K位排序,最后的顺序就是将序列按照最低位为第一关键字，次低位为第二关键字...来排序（或者说是拆分数组），一个数的排名就是它二进制拆分的倒序，于是我们在程序开始预处理每个数字的二进制翻转：

$re{v}_i=[(re{v}_{i>>1})>>1]|(i\mathrm{\&}1)<<d$

(其中$d$为最大数的二进制总位数)

在$FFT$开始前预处理一遍数组，若$i<rev[i]$，则$swap(x[i],x[rev[i]])$即可，这样可以做到不重不漏地将每个数翻转一次，然后直接按照下标分治即可

Code

#include<bits/stdc++.h>
using namespace std;
const int N = 3e6 + 5;
const double PI = acos(-1.0);
int n,m,rev[N],tt = 1,tw = 0;
struct Complex{
    double r,c;
}a[N],b[N];
Complex operator +(Complex x,Complex y)
{
    Complex z;
    z.r = x.r + y.r;
    z.c = x.c + y.c;
    return z;
}
Complex operator -(Complex x,Complex y)
{
    Complex z;
    z.r = x.r - y.r;
    z.c = x.c - y.c;
    return z;
}
Complex operator *(Complex x,Complex y)
{
    Complex z;
    z.r = x.r * y.r - x.c * y.c;
    z.c = x.c * y.r + x.r * y.c;
    return z;
}
inline int read()
{
    int s = 0,w = 1;
    char k = getchar();
    while(k > '9' || k < '0')
    {
        if(k == '-') w = -w;
        k = getchar();
    }
    while(k <= '9' && k >= '0')
    {
        s = s * 10 + k - '0';
        k = getchar();
    }
    return s * w;
}
inline void FFT(Complex *x,int len,int type)
{
    for(int i = 0;i < tt;i++) if(i < rev[i]) swap(x[i],x[rev[i]]);
    for(int mid = 1;mid < tt;mid <<= 1)
    {
        Complex omega;omega.r = cos(PI / mid);omega.c = type * sin(PI / mid);
        for(int j = 0,R = mid << 1;j < tt;j += R)
        {
            Complex now;now.c = 0;now.r = 1;
            for(int k = j;k < j + mid;k++,now = now * omega)
            {
                Complex X = x[k],Y = x[k + mid];
                x[k] = X + now * Y;
                x[k + mid] = X - now * Y;
            }
        }
    }
}
int main()
{
    n = read();m = read();
    for(int i = 0;i <= n;i++) a[i].r = read();
    for(int i = 0;i <= m;i++) b[i].r = read();
    while(tt <= n + m) tt <<= 1,tw++;
    rev[0] = 0;
    for(int i = 1;i < tt;i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (tw - 1));
    FFT(a,n,1);FFT(b,m,1);
    for(int i = 0;i <= tt;i++) a[i] = a[i] * b[i];
    FFT(a,n,-1);
    for(int i = 0;i <= n + m;i++)
        printf("%d ",(int)(a[i].r / tt + 0.5));
    return 0;
}

买一送一环节

我们知道，FFT中暴力计算$\omega$的值，虽然使用了$double$，但是精度仍然会丢失很多，所以我们就要对其进行魔改 $————NTT$ 快速数论变换

考虑到在$FFT$中，每一层的单位乘积是${\omega }^k$，考虑用原根替换这个东西

考虑到答案中的每个数不会很大，我们将答案模一个素数，就可以利用它的原根

(在$int$范围内，我们一般取998244353，其最小原根是3，原根逆元为332748118，并且

$998244353=119\ast 2^{23}+1$，足够我们完成项数在$1e6$以内的运算

于是我们每一步将${\omega }_k$换作$g^{\frac{(p-1)}{k}}$即可($k$表示当前枚举的部分长度)，作为单位乘积每次乘上。

（这就是为什么需要$P-1$有很多个$2$因子）

Code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 3e6 + 5,MOD = 998244353,g = 3,gi = 332748118;
inline ll ksm(ll base,ll pts)
{
    ll ret = 1;
    for(;pts > 0;pts >>= 1,base = base * base % MOD)
        if(pts & 1)
            ret = ret * base % MOD;
    return ret;
}
int rev[N];
ll a[N],b[N],n,m,tt = 1,tw = 0;
inline void NTT(ll *x,ll type)
{
    for(int i = 0;i < tt;i++) if(i < rev[i]) swap(x[i],x[rev[i]]);
    for(int mid = 1,t = 0;mid < tt;mid <<= 1,t++)
    {
        ll dom = ksm((type == 1) ? g : gi,(MOD - 1) >> (t + 1));
        for(int j = 0,R = mid << 1;j < tt;j += R)
        {
            ll w = 1;
            for(int k = j;k < j + mid;k++,w = w * dom % MOD)
            {
                ll X = x[k],Y = x[k + mid] * w % MOD;
                x[k] = (X + Y) % MOD;
                x[k + mid] = (X - Y + MOD) % MOD;
            }
        }
    }
}
int main()
{
    cin>>n>>m;
    for(int i = 0;i <= n;i++) cin>>a[i];
    for(int i = 0;i <= m;i++) cin>>b[i];
    while(tt <= n + m) tt <<= 1,tw++;
    for(int i = 0;i < tt;i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (tw - 1));
    NTT(a,1);NTT(b,1);
    for(int i = 0;i <= tt;i++) a[i] = a[i] * b[i] % MOD;
    NTT(a,-1);
    ll inv = ksm(tt,MOD - 2);
    for(int i = 0;i <= n + m;i++) cout<<a[i] * inv % MOD<<" ";
    return 0;
}