zoukankan      html  css  js  c++  java
  • FFT和NTT学习笔记_基础

    FFT和NTT学习笔记


    算法导论

    参考(贺)

    http://picks.logdown.com/posts/177631-fast-fourier-transform

    https://blog.csdn.net/qq_38944163/article/details/81835205

    https://www.cnblogs.com/RabbitHu/p/FFT.html



    概述

    目的

    (O(nlg_n))的时间复杂度计算多项式乘法

    多项式的表达

    • 系数表达: ({a_0, a_1, ..., a_{n-1}})
    • 点值表达: ({(x_0,y_0), (x_1,y_1), ..., (x_{n-1},y_{n-1})})
      • 插值多项式唯一性: (x)各不相同得到唯一的多项式

    策略

    graph TD A[A,B的系数表达] -- 普通乘法n^2 --> B[C的系数表达] A[A,B的系数表达] -- 扩展成2n并求值nlgn --> A2[A,B在2n个单位复数根的点值] A2[A,B在2n个单位复数根的点值] -- 点值乘法n --> B2[C在2n个单位复数根的点值] B2[C在2n个单位复数根的点值] -- 插值nlgn --> B[C的系数表达]

    关于单位复数根

    复平面

    复数平面(complex plane)是用水平的实轴与垂直的虚轴建立起来的复数的几何表示。
    一个复数的实部用沿着 x-轴 的位移表示,虚部用沿着 y-轴 的位移表示。

    复数乘法: 模长相乘, 幅角相加(幅角:向量与x轴正向夹角)

    单位复数根

    感性理解:

    将复平面上的单位圆从 x 轴起分成 n 等分, 从 x 轴起标号 0..n-1

    那么由 “模长相乘, 幅角相加” 可得, 标号为 1 的那个复数的 (k) 次方, 就是标号为 (k)

    所以将这些数记为 (omega_n^0, omega_n^1, dots, omega_n^{n-1})

    理性证明:

    • (n)次单位复数根 是满足 (omega^n=1) 的复数 (omega), 恰好有 (n)

    (e^{i heta}=cos( heta)+isin( heta))

    (e^{2pi i}=cos(2pi)+isin(2pi)=1)

    • 所以对于 (k=0,1,...n-1) , 这些根是 (e^{2pi i k/n})
    • 有复数乘法可得这 (n) 个单位复数根均匀分布在以复平面原点为圆心的单位半径的圆周上
    • (omega_n=e^{2pi i/n}) 称为 (n)次单位根,其他单位根都是 (omega_n) 的幂次,并有 (omega_n^{j+k}=omega_n^{(j+k) mod n})

    单位复数根基本性质

    (消去引理)

    • 对任何整数 (以及ngeq0,kgeq0,以及dgeq0) , (omega_{dn}^{dk}=omega_n^k)
    • 感性理解: 原来分成 (n) 等分, 取第 (k) 个, 和分成 (dn) 等分, 取第 (dk) 个当然一样
    • 证明:根据指数形式定义式

    (折半引理)

    感性理解:

    (omega_n^{k + n/2} = -omega_n^k)
    复平面单位圆上关于原点对称

    理性证明:

    • 如果 (n>0) 为偶数,那么 (n)(n) 次单位复数根的平方的集合就是 (n/2)(n/2) 次单位复数根的集合
    • 证明((omega_n^{k+n/2})^2=omega_n^{2k+n}=omega_n^{2k}omega_n^n=(omega_n^k)^2=omega_{n/2}^k)
      • 也可以根据 (omega_n^{n/2}=-1) 得到 (omega_n^{k+n/2}=-omega_n^k)
    • 由此可见他可以递归让子问题的规模缩小一半

    (求和引理)

    [sum_{j=0}^{n-1}(omega_n^k)^j = egin{cases} n, & k = 0 \ 0, & k eq 0 \ end{cases} ]

    感性理解:

    (n)(2) 的次幂
    (k=0) 显然
    (k eq 0)
    首先一个偶数 (n) 的等分点的和是 (0) (对称点), 那么对于 (n=2^k) 显然适用
    假如 (k)(n) 的因数, 那么 (k) 也是 (2) 的因数, 它会取一个等分的循环取若干次, 和是 (0)
    否则互质, 那么它会取遍这 (n) 个点, 和也是 (0)

    • 证明:等比数列求和 (sum_{j=0}^{n-1}(omega_n^k)^j=frac{1*(1-omega_n^{kn})}{1-omega_n^k}=0)
    • 因为要求 (k) 不能被 (n) 整除,保证分母不为 (0)

    概述

    系数转点值: 求出 (n) 个单位根的点值 ((y_0, y_1, dots, y_{n-1}))

    点值转系数: 将这 (n) 个点值作为一个新的多项式 (B(x)) 的系数, 用单位根的倒数求一次点值 ((z_0, z_1, dots, z_{n-1}))
    展开可得

    [egin{align*} z_k &= sum_{i = 0}^{n - 1} y_i(omega_n^{-k})^i \ &= sum_{i = 0}^{n - 1}(sum_{j = 0}^{n - 1} a_j(omega_n^i)^j)(omega_n^{-k})^i \ &= sum_{j = 0}^{n - 1}a_j(sum_{i = 0}^{n - 1}(omega_n^{j - k})^i) end{align*} ]

    根据之前的求和引理可得: (z_k = n cdot a_k)
    于是系数 (a_k = z_k / n)

    也就是说, 这两个过程都可以通过求点值来完成,
    所以要完成两个多项式的相乘, 只需要先求一遍他们的点值, 点值相乘, 在转成系数就行了

    DFT

    点值向量(y=(y_0, y_1, ..., y_{n-1}))就是系数向量(a=(a_0, a_1, ..., a_{n-1}))离散傅里叶变换(DFT), 也记为(y=DFT_n(a))

    FFT

    快速傅里叶变换(FFT),利用单位复数根的特殊性质,在(O(nlg_n))的时间内计算出(DFT_n(a))

    分治策略,采用(A(x))中偶数下标和奇数下标的系数,分别定义两个次数界为(n/2)的多项式

    • (A_0(x) = a0 + a2*x + ... + a_{n-2}*x^{n/2-1})
    • (A_1(x) = a1 + a3*x + ... + a_{n-1}*x^{n/2-1})

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

    带入单位根可以发现

    [egin{align*} A(omega_n^k) &= A_0(omega_n^{2k}) + omega_n^kA_1(omega_n^{2k}) \ &= A_0(omega_{frac{n}{2}}^{k}) + omega_n^kA_1(omega_{frac{n}{2}}^{k}) end{align*} ]

    [egin{align*} A(omega_n^{k + frac{n}{2}}) &= A_0(omega_n^{2k + n}) + omega_n^{k + frac{n}{2}}A_1(omega_n^{2k + n}) \ &= A_0(omega_{frac{n}{2}}^{k} imes omega_n^n) + omega_n^{k + frac{n}{2}} A_1(omega_{frac{n}{2}}^{k} imes omega_n^n) \ &= A_0(omega_{frac{n}{2}}^{k}) - omega_n^kA_1(omega_{frac{n}{2}}^{k}) end{align*} ]

    这两个关于远点对称的单位根的点值只是相差一个符号, 所以计算 (A_0(omega_{frac{n}{2}}^{k}), omega_n^kA_1(omega_{frac{n}{2}}^{k})) 就可得到两个单位根的取值

    因为应用了(omega_n^k)的正负数形式,所以称其为旋转因子

    可得到以下伪代码,时间复杂度(O(nlg_n))

    RECURSIVE_FFT(a) // 计算DFTn(a)
        n = a.length
        if n == 1
            return a // a * (x^0)
        a[0] = (a(0), a(2), ..., a(n-2))
        a[1] = (a(1), a(3), ..., a(n-1))
        y[0] = RECURSIVE_FFT(a[0])
        y[1] = RECURSIVE_FFT(a[1])
        for k = 0 to n/2-1
            y[k] = y[0][k] + omega[n][k] * y[1][k]
            y[k+n/2] = y[0][k] - omega[n][k] * y[1][k]
        return y;
    

    IDFT

    通过点值得到系数

    IFFT

    在概述中说明了, 只需要将点值作为系数的到新的多项式, 再带入单位根的倒数, 得到这个多项式的点值

    然后点值 (/n) 即可得到系数


    高效FFT实现

    画出系数的递归树可得叶子的标号的规律:

    蝴蝶操作中

    graph LR id1[y0k] --> A id2[+] --> ret1[y0k+omega*y1k] id3[y1k] --> B id4[-] --> ret2[y0k-omega*y1k] A --> id2[+] A --> id4[-] B --> id2[+] B --> id4[-] ret1[y0k + omega*y1k] --> ret3[y k] ret2[y0k - omega*y1k] --> ret4[y k+n/2]

    可以观察到(y_k^{[1]})的位置其实就是(y_{k+n/2})的位置,相当于两个位置上的(y)值进行蝴蝶操作并只对这两个位置进行修改

    观察下图递归FFT输入的向量树

    graph TD id1[a0,a1,a2,a3,a4,a5,a6,a7] --- id2[a0,a2,a4,a6] id1[a0,a1,a2,a3,a4,a5,a6,a7] --- id3[a1,a3,a5,a7] id2[a0,a2,a4,a6] --- id4[a0,a4] id2[a0,a2,a4,a6] --- id5[a2,a6] id3[a1,a3,a5,a7] --- id6[a1,a5] id3[a1,a3,a5,a7] --- id7[a3,a7] id4[a0,a4] --- id8[a0,000] id4[a0,a4] --- id9[a4,100] id5[a2,a6] --- id10[a2,010] id5[a2,a6] --- id11[a6,110] id6[a1,a5] --- id12[a1,001] id6[a1,a5] --- id13[a5,101] id7[a3,a7] --- id14[a3,011] id7[a3,a7] --- id15[a7,111]

    其叶子有这样的规律:

    (rev(x))(x)的二进制表示的倒串所形成的(lg_n)位的数,则(rev(x))位上为(a_x)

    因为每次的叶子是一样的( (n) 相同), 所以可以预处理, 并且可以 (O(n)) 递推 rev[]

    注意 预处理叶子位置不是扫到 (len / 2)

    不妨从下往上模拟递归的过程

    方便起见(省去二维的数组), 将每一层中各多项式的各点值标号, 将递归的写成这样

    {
        FFT(a, n / 2);
        FFT(a + n / 2, n / 2);
        for (int i = 0; i < n / 2; ++ i) // 枚举 omega_n^i
        {
            buf[i] = a[i] + a[i + m];
            buf[i + m] = a[i] - omega_{n / 2}^{-i} * a[i + m];
        }
    }
    

    转成非递归, 现将原系数变成叶子的状态

    for (int l = 2; l <= len; l <<= 1)
        {
    	int ll = l >> 1;
    	for (int i = 0; i < len; i += l)
    	{
    	    for (int j = 0; j < l / 2; ++ j)
    	    {
    		CPD x = ome[len / l * j];
    		if (opt) x = inv[len / l * j];
    		buf[i + j] = a[i + j] + x * a[i + j + ll];
    		buf[i + j + ll] = a[i + j] - x * a[i + j + ll];
    	    }
    	}
    	for (int i = 0; i < len; ++ i) a[i] = buf[i];
        }
    

    注意这里的单位根可以利用引理对应到 (n) 次单位根上

    卡常提示:

    • 然后注意不要预处理 ome[], inv[], 循环的时候直接乘上第一个单位根即可
    • 手写复数, 并且不要用构造函数(常数小了一半)
    • 2 * PI / l 可以改为 PI / mid
    // P1919 A*B(FFT)
    //#pragma GCC optimize(2)
    #include <iostream>
    #include <cstdio>
    #include <cmath>
    #include <cstring>
    #include <algorithm>
    
    using namespace std;
    
    struct Complex
    {
        double x, y;
        Complex operator + (const Complex &b) const { return (Complex){x + b.x, y + b.y}; }
        Complex operator - (const Complex &b) const { return (Complex){x - b.x, y - b.y}; }
        Complex operator * (const Complex &b) const { return (Complex){x * b.x - y * b.y, x * b.y + y * b.x}; }
    } ;
    
    const int MAXN = (1 << 21) + 3;
    char ch[MAXN];
    
    int lena, lenb, n, dgt;
    Complex a[MAXN], b[MAXN];
    void read(Complex * a, int & len)
    {
        scanf("%s", ch + 1); len = strlen(ch + 1);
        for (int i = 0; i < len; ++ i) a[i].x = ch[len - i] - '0';
    }
    
    int rev[MAXN];
    void init(int n, int dgt)
    {
        for (int i = 0; i < n; ++ i) 
    	rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (dgt - 1));
    }
    
    const double PI = acos(-1.0);
    void FFT(Complex * a, int len, int opt)
    {
        for (int i = 0; i < len; ++ i) 
    	if (i < rev[i]) swap(a[i], a[rev[i]]);
        for (int mid = 1; mid < len; mid <<= 1)
        {
    	Complex omen = (Complex){cos(PI / mid), opt * sin(PI / mid)} ;
    	for (int i = 0; i < len; i += (mid << 1))
    	{
    	    Complex ome = (Complex){1, 0} ;
    	    for (int j = 0; j < mid; ++ j, ome = ome * omen)
    	    {
    		Complex t = ome * a[i + j + mid];
    		a[i + j + mid] = a[i + j] - t;
    		a[i + j] = a[i + j] + t;
    	    }
    	}
        }
    }
    
    int ret[MAXN];
    
    /*
    20191212
    0724~0757~0839
    a * b FFT
     */
    
    int main()
    {
        read(a, lena); read(b, lenb);
        n = 1; dgt = 0;
        while (n < lena + lenb) n <<= 1, ++ dgt;
        init(n, dgt);
        FFT(a, n, 1); 
        FFT(b, n, 1);
        for (int i = 0; i < n; ++ i) a[i] = a[i] * b[i];
        FFT(a, n, -1);
        for (int i = 0; i < n; ++ i) 
        {
    	ret[i] += (int)(a[i].x / n + 0.5);
    	ret[i + 1] += ret[i] / 10;
    	ret[i] %= 10;
        }
        while (!ret[n - 1] && n > 0) -- n;
        for (int i = n - 1; i >= 0; -- i) printf("%d", ret[i]);
    //    printf("#%d
    ", cnt);
        return 0;
    }
    /*
    76543210
    76543210
    
    6543
    21
     */
    

    总结

    FFT 利用了单位根的以下性质

    from http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-17

    首先来看 FFT 中能在 O(nlogn) 时间内变换用到了单位根 ω 的什么性质 ...

    自己解释一遍:

    1. 单位根 (omega_n^0, dots, omega_n^{n-1}) 是单位圆的等分点(不同), 所以有了 "消去引理"
    2. (omega_n^{frac{n}{2}+k}=-omega_n^k, omega_n^2=omega_{frac{n}{2}})
      • 根据奇偶系数拆分的表达式, 这一层的 (m) 个多项式都要带入 (n) 个值, 可以通过下一层 (2m) 个多项式带入 (n/2) 个值求得, 由此可知 (log)
    3. (sum_{j=0}^{n-1}(omega_n^k)^j = egin{cases} n, & k = 0 \ 0, & k eq 0 \ end{cases}) 使得点值转系数的逆变换可以同样方式完成

    其中第一点是最基本的, 他推出了后面的性质

    NTT

    在取一些模数的意义下, 原根具有上述性质

    区别在于原根的手动计算不同

    void NTT(LL * a, int len, int sgn)
    {
        for (int i = 0; i < len; ++ i) 
    	if (i < rev[i]) swap(a[i], a[rev[i]]);
        for (int mid = 1; mid < len; mid <<= 1)
        {
    	for (int i = 0; i < len; i += (mid << 1))
    	{
    	    for (int j = 0; j < mid; ++ j)
    	    {
    		LL t = (sgn == 1 ? g[len / mid / 2 * j] : ig[len / mid / 2 * j]) * a[i + j + mid] % MOD;
    		a[i + j + mid] = (a[i + j] + MOD - t) % MOD;
    		a[i + j] = (a[i + j] + t) % MOD;
    	    }
    	}
        }
    }
    
    invn = inv(n); 
    g[0] = ig[0] = 1; 
    g[1] = fpow(3, (MOD - 1) / n); ig[1] = inv(g[1]);
    for (int i = 2; i < n; ++ i) 
        g[i] = g[i - 1] * g[1] % MOD, ig[i] = ig[i - 1] * ig[1] % MOD;
    
  • 相关阅读:
    wget(转)
    852. Peak Index in a Mountain Array
    617. Merge Two Binary Trees
    814. Binary Tree Pruning
    657. Judge Route Circle
    861. Score After Flipping Matrix
    832. Flipping an Image
    461. Hamming Distance
    654. Maximum Binary Tree
    804. Unique Morse Code Words
  • 原文地址:https://www.cnblogs.com/Kuonji/p/10354302.html
Copyright © 2011-2022 走看看