FTT 快速傅里叶变换

高效实现 DFT 离散傅里叶变换的一种手段，可以将两个 $n$ 次多项式相乘的时间复杂度优化到 $O(n\log n)$ 的算法

引理：给定 $n$ 个点可以确定 $n-1$ 次函数曲线的系数和系数存在映射关系

设两个多项式分别为 $F(x),G(x)$

  • 通过这两个多项式的系数经过系数乘法可以计算出 $F(x)\cdot G(x)$ 的系数，基于单位圆划分出 $n$ 个点，通过计算得到的这 $n$ 个点值，我们就能唯一确定函数 $F(x)\cdot G(x)$
  • 基于单位圆划分出 $n$ 个点，计算这两个多项式的这 $n$ 个点值，再利用点值乘法确定函数 $F(x)\cdot G(x)$

单位圆

image-20250807220538390

复平面的单位元圆上的点都可以被表示为：

$$ z=\cos \theta +i\sin \theta\ (0\le \theta < 2\pi) $$

把单元圆 $n$ 等分（$n=2^{b\in N^+}$），可以得到方程 $z^n=1$ 的 $n$ 个复数根 $\omega_n$ ，单位根为：

$$ \omega_n^k=\cos \frac{2\pi k}{n}+i\sin\frac{2\pi k}{n} $$

满足性质有：$\omega _n^k\cdot \omega_n^m=\omega_n^{k+m},(\omega _n^k)^m=\omega _n^{km}$ ，周期性：$\omega _n^{k+n}=\omega _n^k$ ，对称性：$\omega _n^{k+\frac{n}{2}}=-\omega _n^k$ ，折半性：$\omega _n^{2k}=\omega _{\frac{n}{2}}^k$

正变换

通过多项式 $A(x)$ 的系数求点值，设系数为 $a_0,a_1,\cdots,a_{n-1}$ ，则 $A(x)$ 被表示为：

$$ A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1} x^{n-1} $$

将奇次项和偶次项分离后有：

$$ A(x)=(a_0+a_2x^2+\cdots +a_{n-2}x^{n-2})+(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2})x $$

分别设 $A_0(x),A_1(x)$ 为：

$$ A_0(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1} $$

$$ A_1(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1} $$

$$ A(x)=A_0(x^2)+A_1(x^2)x $$

将 $$x=\omega _n^k\ (k<\frac{n}{2})$$ 代入得：

$$ A(\omega_n^k)=A_0(\omega_n^{2k})+A_1(\omega_n^{2k})\omega_n^k=A_0(\omega_\frac{n}{2}^{k})+A_1(\omega_\frac{n}{2}^{k})\omega_n^k $$

同理代入 $$x=\omega _n^{k+\frac{n}{2}}\ (k<\frac{n}{2})$$

$$ A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_n^{2k+n})+A_1(\omega_n^{2k+n})\omega_n^{k+\frac{n}{2}}=A_0(\omega_\frac{n}{2}^{k})-A_1(\omega_\frac{n}{2}^{k})\omega_n^{k} $$

在求出 $A_0(\omega _n^k),A_1(\omega _n^k)$ 后我们就能同时得到 $A(\omega _n^k),A(\omega _n^{k+\frac{n}{2}})$ ，所以对 $A_0,A_1$ 分别进行递归求解

特别的，基于分治 FFT 处理的多项式长度只能是 $2^{b\in N^+}$ ，所以对于不是二整数次幂的形式，都需要补齐再计算

逆变换

通过多项式 $A(x)$ 的点值求系数，设代入 $\omega _n^0,\omega _n^1,\cdots,\omega _n^{n-1}$ 得到的点值为 $y_0,y_1,\cdots,y_{n-1}\ (y_i=\sum\limits_{j=0}^{n-1} a_j(\omega_n^i)^j)$

构造多项式 $B(x)$：

$$ B(x)=y_0+y_1x+\cdots+y_{n-1}x^{n-1} $$

将单位根的倒数 $\omega _n^0,\omega _n^{-1},\cdots,\omega _n^{-(n-1)}$

$$ \omega _n^k=\cos \theta +i\sin \theta\\ \frac{1}{\omega _n^k}=\frac{1}{\cos \theta +i\sin \theta}=\cos \theta -i\sin \theta=\cos (-\theta) +i\sin (-\theta) $$

代入 $B(x)$ 得到新点值为 $z_0,z_1,\cdots,z_{n-1}$，其中：

$$ z_k=\sum\limits_{i=0}^{n-1} y_i(\omega_n^{-k})^i=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1} a_j(\omega_n^i)^j(\omega_n^{-k})^i=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i $$

$$ \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i= \left\{ \begin{array}{l} \sum\limits_{i=0}^{n-1}1=n & j=k\\ (\omega_n^{j-k})^n=(\omega_n^n)^{j-k}=0 & j\ne k \end{array} \right. $$

所以 $z_k=na_k$ ，通过 $\omega_n^{-k}$ 确定的点值逆变换后可以确定原函数的系数

基于递归实现的 FFT

void FFT(vector<complex<double>> &A,int n,int op){
    if (n==1) return;
    vector<complex<double>> A0(n/2),A1(n/2);
    for (int i=0;i<n/2;i++)
        A0[i]=A[i*2],A1[i]=A[i*2+1];
    FFT(A0,n/2,op),FFT(A1,n/2,op);
    complex<double> w1({cos(2*PI/n),sin(2*PI/n)*op});
    complex<double> wk({1,0});
    for (int i=0;i<n/2;i++){
        A[i]=A0[i]+A1[i]*wk;
        A[i+n/2]=A0[i]-A1[i]*wk;
        wk=wk*w1;
    }
}

op 参数表示当前执行的是正变换还是逆变换，区别在于对 $\omega_n^k$ 的符号控制

基于迭代实现的 FFT

image-20250808142932444

原序列分治到底层后变换的序列存在一个规律：

  • 每个数用二进制表示，将二进制翻转对称一下，就是最终那个位置的下标，例如 $a_0=001_2,a_4=100_2$ ，我们称这个变换为位逆序置换

根据定义的方法，我们可以在 $O(n\log n)$ 的时间复杂度内求出每个数变换后的结果

事实上，还可以通过 $O(n)$ 的递推时间复杂度实现，设 $n=2^k$ ，$k$ 表示二进制数的长度

设 $R_x$ 为长度为 $k$ 的二进制数 $x$ 翻转后的结果，存在这样的递推关系：

$$ R_x=\frac{R_{\frac{x}{2}}}{2} +(x\ \mathrm{mod}\ 2)\times \frac{n}{2} $$

从小到大求 $R_x$ ，$R_{\frac{x}{2}}$ 的值已知。把 $x$ 除以二（右移一位）然后翻转再除以二，就得到了 $x$ 除开二进制个位之外其他位翻转的结果，再单独考虑二进制个位的影响，如果翻转前个位是 $1$ ，那么翻转后的最高位也为 $1$

void change(vector<complex<double>> &A,int n){
    vector<int> R(n,0);
    for (int i=0;i<n;i++)
        R[i]=(R[i>>1]>>1)+((i&1)?(n>>1):0);
    for (int i=0;i<n;i++)
        if (i<R[i]) swap(A[i],A[R[i]]);
}

利用位逆序变换得到后序列，然后自底向上合并

void FFT(vector<complex<double>> &A,int n,int op){
    change(A,n);
    for (int m=2;m<=n;m<<=1){//枚举分治块的长度
        complex<double> w1({cos(2*PI/m),sin(2*PI/m)*op});
        for (int i=0;i<n;i+=m){//枚举块的个数
            complex<double> wk({1,0});
            for (int j=0;j<m/2;j++){//
                complex<double> x=A[i+j],y=A[i+j+m/2]*wk;
                A[i+j]=x+y;
                A[i+j+m/2]=x-y;
                wk=wk*w1;
            }
        }
    }
}

优化高精度乘法

void solve(){
    string sa,sb;
    cin>>sa>>sb;
    int n=sa.size(),m=sb.size();
    for (int i=n-1;i>=0;i--) a[n-i-1]=sa[i]-'0';
    for (int i=m-1;i>=0;i--) b[m-i-1]=sb[i]-'0';
    int nm=n+m;
    for (n=1;n<=nm;n<<=1);
    FFT(a,n,1),FFT(b,n,1);
    for (int i=0;i<n;i++) ab[i]=a[i]*b[i];
    FFT(ab,n,-1);
    vector<int> ans;
    for (int i=0,t=0;i<n;i++){
        t+=(ab[i].real()/n+0.5);
        ans.push_back(t%10);
        t/=10;
    }
    while (ans.size()>1&&ans.back()==0) ans.pop_back();
    for (int i=ans.size()-1;i>=0;i--) cout<<ans[i];
}

将十进制数看作一个 $x=10$ 的函数 $A(x)$ ，每个系数 $a_{i}$ 表示第 $i+1$ 位上的数

同样的，将读入的两个数 $a,b$ 转化为数组存储后利用 FFT 计算出点值，然后通过点值乘法计算 $a\times b$ 映射的点值

再通过 FFT 逆变换计算出系数，最后按位取模得到拆分下来的十进制的每一位