快速傅里叶变换|快速数论变换

FFT——快速傅里叶变换

什么是FFT

(FFT)全称((Fast Fourier Transformation))，中文名：快速傅里叶离散变换。

这个名字听起来很高级，实际上也很高级，它是(DFT)和(IDFT)的加速版本，用于加速多项式乘法。

  • 接下来我们把某个多项式记为“大写字母+元素”，如多项式：

(A(x)=a_nx^n+a_{n-1}x^{n-1}+…+a_2x^2+a_1x^1+a_0)

(B(x)=b_mx^m+b_{m-1}x^{m-1}+…+b_2x^2+b_1x^1+b_0)

你也可以理解为(A/B)是关于(x)的函数

现在，我们要求(A(x)*B(x))，即两个多项式的乘积(卷积)，那么我们很容易想到的就是(O(N^2))的做法：把每一项依次分别相乘再相加。

[sum^{n+m-1}_{i=0}(sum^{i}_{j=0}a_i*b_{i-j})*x^i ]

这样的时间复杂度明显不尽人意，那么(FFT)就是把这种算式计算优化到(NlogN)的算法。

前置知识

  • 复数

我们定义虚数 (i)， 且 (i^2 = -1)，即 (i = sqrt {-1})，那么复数就形如 (a+bi)，(i) 是一个虚数单位，(a) 叫做复数的实部， (b) 叫做复数的虚部，那么 (a+bi) 在复平面（高中会讲）可以类似于平面直角坐标系，表示为点 ((a, b))。

复数的运算法则

  1. 加法：实部和虚部分别相加。

  2. 减法：实部和虚部分别相减。

  3. 乘法：像一次多项式一样相乘。

  4. 除法：分子分母同乘分母的共轭。

共轭：两个复数实部相同，虚部相反。

上面除法分子分母同乘分母的共轭，把分母有理化成实数。

    struct complex { //复数 
        long double a, b;
        complex() { a = b = 0;}
        //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    };

这边大家直接记一个重要结论：

复数相乘，辐角相加，模长相乘。

单位根

对于复数(w),若整数 (n) 满足 (w^n=1)，则称 (w) 是 (n) 次单位根。

扔个单位圆：

为什么扔个单位圆呢？因为单位圆上任意一点到原点的距离，即模长都为(1)。

那么，对于一个复数，只有模长为(1)的复数才有可能成为(n)次单位根!

现在，我们在复平面上把复数锁定在了这个单位圆上。

设复数(w)模长为(q)，其幅角为(mπ)，(m)为([0,2])间实数，那么，其(n)次方幅角为(nmπ)。

则有(nmπ=kπ,k∈N)，即(nm=k)。

  1. 若(k=0),因为(n>0)，所以方程有且只有(m=0)一解。

  2. 若(k>0)，则(k=1,2,3,…)，易知(m)有(n-1)个不重解

综上，复数(w)的(n)次单位根有(n)个。

并且，复数的(n)次单位根(n)等分单位圆。

单位根的书写：

单位根用符号(ω)表示，从1开始，沿单位圆逆时针把单位根从(0)开始标号。

然后单位根还有个神奇的东西(ω^{k}_{n}=ω_{n}^{kmod n})

单位根的性质

  1. 对于任意(n)，(ω_n^0=1)。

  2. (ω_n^k=(ω_n^1)^k)。

  3. ((ω_n^m)^k=(ω_n^k)^m)

  4. (ω_n^i*ω_n^j=ω_{n}^{i+j})。

  5. (ω_{2n}^{2m}=ω_n^m)。

  6. 若 (2 mid n)，(ω_n^{m+frac{n}{2}}=-ω_n^m)

多项式的表达方法

  1. 系数表达法：这个如同我们上面写的多项式即为系数表达法。

  2. 点值表达法

定义多项式(F(x))：

我们知道，(1)个n元一次方程最少需要(n)个式子求解，而我们把系数表达(F(x))抽象为一个函数(F(x))，那么这个函数至少需要(n+1)个点确定下来，那么我们任取(n+1)个不同的数构成集合(U={k_1,k_2,…,k_{n+1}})，对(F(x))分别求值得到一个新的集合，那么将两个集合合并成点集，有：

[T={(k_1,F(k_1)),k_2,F(k_2),…,(k_{n+1},F(n+1))} ]

这便是(F(x))的点值表达式，并且，点值表达法下的乘法运算更加简单：

多项式(P,Q)分别取点((x,y_1))和((x,y_2))，则(P*Q)可得到点((x,y_1*y_2))。令(S=P*Q)，则(C(x_i)=y_{1_i}*y_{2_i})。

可见，点值表达法下，多项式乘法是(O(N))的。

  • （补充）乘法本质：卷积

根据我们上面暴力计算(A(x)*B(x))，不妨设(C=A*B)，那么：

(C[k]=sum^{k}_{i=0}A[i]B[k-i])

(A[i],B[k-1])分别为(A)中(x^i)的系数和(B)中(x^{k-i})的系数，那么我们知道：(x^i*x^{k-i}=x^k)，故所有的(x^i*x^{k-i})的系数对(x)的次项系数均有贡献，需要累加。

那么就可以写成：(C[k]=sumlimits_{i+j=k}A[i]B[j])

那么形如(C[k]=sumlimits_{i⊕j=k}A[i]B[j])((⊕)为某种运算)的式子称为卷积。

那么，多项式的乘法就是对于加法的卷积。

我们仍然可以对多项式相乘的模板打个暴力：

    #include<bits/stdc++.h>
    using namespace std;

    const int SIZE = 1e6 + 10;

    long long A[SIZE], B[SIZE], C[SIZE];
    int n, m;

    int main() {
        scanf("%d %d", &n, &m);
        n++, m++;
        for (int i = 0; i < n; i++) scanf("%lld", A + i);
        for (int i = 0; i < m; i++) scanf("%lld", B + i);
        for (int k = 0; k < n+m-1; k++)
            for (int i = 0; i <= k; i++)
                C[k] += A[i] * B[k - i];
        for (int i = 0; i< n + m - 1; i++) 
            printf("%lld ", C[i]);
        return 0;
    }

模板题亲测 44 分。

算法理论

上面的主要部分还是写了一个(O(N^2))的暴力，那我们要寻找一个更快的办法：

就是它！(FFT)算法

前面我们已经说过：(FFT)算法是一个(O(NlogN))的算法，但常数较大(相信你们会卡常)，并且，(FFT)是(DFT)和(IDFT)的快速算法：

思想：

首先你要会点值表达式，当然我们已经讲过。

然后我们来看(A(x)=x^2+3x-2)（红色）与(B(x)=-x^2+6x+4)(绿色)的取点得到点的点值表达式:

在(A(x))上取三个点（紫色）((-4,2)(-2,-4)(1,2))；在(B(x))上取三个点((0,4)(2,12)(6,4))，然后转为系数表达(O(N^2))相乘然后再带入(x)得到(C)的点值表达式？这时绝对不行的，我们上面说过，点值表达式相乘可以在(O(n))内解决，那怎么办呢？

我们可以取(3)个(x)相同的点：

这样就行了？那是不可能的。我们知道，两个二次项式相乘得到的是一个四次项式，而(C(x)=y_1*y_2)，只能得到三个点，这肯定是不行的。

这很明显我们需要(2n+1)个点，不够怎么办呢？再添两个不就好了？
下面添加了(I,J)(横坐标相同)，(G,H)(横坐标相同)，这样就行了。

那么，我们只需要进行(2n+1)次乘法即可，省略常数（都说了常数巨大了），复杂度(O(N))

神仙问题：(tql) (but) 这有用吗

当然有用啊：

算法流程：

[ ext{系数表达式}implies ext{点值表达式} implies ext{系数表达式} ]

这不就有个用了吗？

但这个思路仅仅是(DFT+IDFT)，而(FFT)是用来加速(DFT)和(IDFT)的。

如果让我们求点值表达式，那我们肯定会带整数进去，因为这样好计算。但是傅里叶，（不知道是不是脑抽了），把复数带进了多项式中，然后神奇的事情就发生了……

我们把多项式如下拆开：

(F(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+…+a_{n-1}x^{n-1}))

（保证(n=2^k,k∈N^*))

然后设两个有(n/2)次项的多项式(FL(x))和(FR(x))

[FL(x)=a_0+a_2x+…+a_{n-2}x^{frac{n}{2}-1} ]

[FR(x)=a_1+a_3x+…+a_{n-1}x^{frac{n}{2}-1} ]

[F(x)=FL(x^2)+xFR(x^2) ]

设(k<frac{n}{2})，把(ω_n^k)代入(F(x))：(F(ω_n^k)=FL(ω_n^{2k})+ω_n^kFR(ω_n^{2k}))

接着，就有(F(ω_n^k)=FL(ω_{n/2}^{k})+ω_n^kFR(ω_{n/2}^{k}))

那么，如果我们知道(FL(x),FR(x))分别在(ω^0_{n/2},ω^1_{n/2},…,ω^{n/2-1}_{n/2})处的点值表示，那么我们就可以(O(N))内求出(F(x))在(ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n})处的点值表示了。

接下来我们的目标就是求(F(x))在(ω^{n/2}_{n},ω^{n/2+1}_{n},…,ω^{n-1}_{n})处的点值表示。

然后……

继续设(k<n/2)，然后把(ω^{k+n/2}_{n})代入(F(x))中。

运用一系列单位根的性质，读者可以自行证明出：

[F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2}) ]

然后我们就可以求出(F(x))在(ω^0_{n},ω^1_{n},…,ω^{n-1}_{n})处的点值表示，接着，我们就实现了(DFT)（听着好累的样子……）

但是我们还不知道(FL(x),FR(x))在(ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n})处的点值表示。难道要暴力代入吗？这样反而会使时间复杂度上升。

但是我们想想：(FL(x))和(FR(x))是由原来问题(F(x))转化来的，怎么转化的？看到这里，你应该想到了分治。我们知道，分治的时间复杂度是(logN)，好像完成了一半了？

然后，因为(FL(x),FR(x))是(F(x))的子问题，那么我们可以对它们继续进行分解成一个个小的子问题，然后再合并，类似于归并排序。

(FFT)分治的合并过程根据：

1. (F(ω^{k}_{n})=FL(ω^{k}_{n/2})+ω^{k}_{n}FR(ω^{k}_{n/2}))

2. (F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2}))

从理论到实现

上文的所有证明都基于(2=2^k,k∈N^*)，但是我们做题的时候并没有这一点，所以我们可以补一些系数为(0)的项，这对计算结果没有影响。

1. 预处理单位根

我们知道一个性质(ω_n^k=(ω_n^1)^k)

那么，我们只要知道(ω_n^1)就可以求出(ω_n^0,ω_n^1,ω_n^2,…,ω_n^{n-1})。

怎么求(ω_n^1)？

解三角形不就好了？已知(|OA|=1),幅角是(frac{1}{n})圆周，那么：

(由于(c++)下三角函数用弧度制表示，以下全部使用弧度制。)

设幅角(α=frac{2π}{n}),，则(ω_n^1(cosig(frac{2π}{n} ig),sinig(frac{2π}{n} ig)))。

这边要牢记，(x)是实轴，(y)是虚轴：

    struct complex { //复数 
        long double a, b;
        complex() { a = b = 0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    }w[SIZE],temp, buf;

    void unit_roots(int n) {
        const long double pi = 3.14159265358979;
        //complex ;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        w[0].a = 1, w[0].b = 0;
        buf.a = 1, buf.b = 0;
        for(int i = 1; i < n ; i++) {
            buf *= temp;
            w[i] = buf; 
        }
        for(int i = 0; i < n; i++)
            printf("w_%d %.3Lf %.3Lf
", i, w[i].a, w[i].b);	
    }

然后我们就可以实现(DFT)了

因为还没讲(IDFT)，先写(DFT):

递归(DFT)实现（复数板子上面有了下面就不写了）

    int n, m;
    Complex f[SIZE];
    void dft(Complex *f, int len) { //当前子问题长度
        if(len == 0) return;
        Complex fl[len + 10], fr[len + 10];
        for(int i = 0 ; i < len; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        dft(fl, len >> 1);
        dft(fr, len >> 1);
        len *= 2;
        Complex temp, buf;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        buf.a = 1;
        for(int i = 0 ; i < len / 2; i++) {
            f[i] = fl[i] + buf * fr[i];
            f[i + len / 2] = fl[i] - buf * fr[i];
            buf *= temp;
        }
    }
    int main() {

        scanf("%d", &n);
        for(int i = 0; i < n; i++) scanf("%Lf", &f[i].a);
        for(m = 1; m < n; m <<= 1);//补到2的幂次方
        dft(f, m >> 1);
        for(int i = 0; i < m; i++) printf("%.3Lf ", f[i].a);
        return 0;
    }

事实上虚数可以用(小心TLE)：

            complex <double> f;

然后我们得到了一堆点值表达……

这并没有什么用

于是，我们现在来学习(IDFT)——(DFT)的逆运算。

刚刚的多项式为(F(x))，系数是(A)数列，点值表示成(M)，(M[k]=sum limits_{i=0}^{n-1}(ω_n^{k})^i*F[i])

然后就有(F[k]=frac{1}{n}sum limits_{i=0}^{n-1}(ω_n^{-k})^i*M[i])

结论记住就好了，可以自己尝试证明。

我们可以先求出和式，然后再除以(n)。

接下来我们求出(ω_n^0,ω_n^{-1},…ω_n^{-n+1})

并且仍有：(ω_n^{-j}=(ω_n^{-1})^j)

所以我们还是求出(ω_0^{-1})即可，并且，(ω_0^{-1})与(ω_0^{1})关于(x)轴对称。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const long double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        long double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
        void conj() { imag = -imag; } 
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }	
    int n, m, p;		
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%Lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%Lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
	return 0;
}

接着TLE了……(丝毫没有卡常的欲望)

因为递归常数巨大，并且很容易爆空间。(但实际上是因为写了(long double)炸了)

下面是AC的递归代码，以后不要写long double 了

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 

        void conj() { imag = -imag;}
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }	
    int n, m, p;		
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    }

实际上跑的很快的递归评测记录。

那有没有更优化的办法呢？这是肯定的。

蝴蝶变换

先看一张图

这张图是什么？这不是我们分治的顺序吗？

我在一开始和结束都打上了二进制，你有没有什么发现？

是的，每一个数的二进制都反过来了。

然后我们的任务就转换成求二进制反序：

    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l)//l是二进制位数
    //rev[i]是i翻转得到的

结论牢记即可。

还记得(DP)吗？（貌似没什么关系）

我们用循环实现(FFT)，也就像(DP)一样定义状态、阶段决策：

  • 第一层循环枚举变换区间大小（阶段），第二层循环枚举开头（决策），第三层处理区间。
    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex operator += (const complex &x) {
            *this = *this + x;
            return *this;
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
    }a[SIZE], b[SIZE];
    int n, m, p, L = -1;
    int rev[SIZE];

    void fft(complex *f, int inv) {
        for(int i = 0 ; i < p ; i++) 
            if(i < rev[i]) 
                swap(f[i], f[rev[i]]);
        for(int len = 2; len <= p; len <<= 1) {
            int length = len >> 1;
            complex temp;
            temp.real = cos(pi / length);
            temp.imag = sin(pi / length);
            if(inv) temp.imag = -temp.imag;
            for(int l = 0; l < p; l += len) {
                complex buf; buf.real = 1;
                for(int k = l; k < l + length; k++) {
                    complex t = buf * f[k + length];
                    f[k + length] = f[k] - t;
                    f[k] += t;
                    buf *= temp;
                }
            }
        }
    }

    int main() {
        scanf("%d %d", &n, &m);
        for(int i = 0; i <= n; i++)
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++)
            scanf("%lf", &b[i].real);
        for(p = 1; p < n + m + 1; p <<= 1, L++);
        //l是位数，因为会多算一位，一开始初值直接给-1
        for(int i = 1; i < p; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
        fft(a, 0), fft(b, 0);
        for(int i = 0; i < p; i++)
            a[i] *= b[i];
        fft(a, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    }

评测记录

NTT——快速数论变换

  • 首先，学习(NTT)算法你要学会(FFT)（理解也可）（刚刚我们才讲完你不会忘了吧……）（我认为我写的很详细了呢(QAQ)）。

  • 接下来，你要确保你会原根这个东西（起码要知道），不过我在后面还会再提，具体可以先看我的这一篇博客：同余相关，在比较下面的部分（由于还待完善，所以只是写了一点）。

  • 然后就是希望大家看了这一篇博客对(NTT)有所了解，可能我讲的不好，但是起码你要知道(NTT)是什么。

前置知识

（怎么又是前置知识）

(FFT)

上面已经说了先理解(FFT)，并且我也给了写的比较详细的博客，这里不再展开。

但其实(FFT)并不是重点，而是里面用于优化的很重要的一个东西：(color{red}{ ext{单位根}})。

  • 原根

但是如果 (FFT) 要求取模呢？

[color{red}{What? ext{你不会？}} ]

好像因为 (FFT) 用的单位根不太好取模（复数），所以我们可能需要找一个东西来替代下单位根，那么就要求这个东西和单位根具有类似的性质。然后科学家们就开始寻找了……（也不知道是谁找到的）接下来他们找到了原根。

首先我们要了解一下啥是原根。（如果你看了上面给出的博客你应该知道了）

原根的定义

其实在那一篇博客里的原根定义很容易让人摸不着头脑……为啥呢？

[color{green}{δ_m(a)=φ(m) ext{？} what is this?} ]

因为它很不好求……那我们不妨换一种说法：

设 (p) 为质数，则根据 (color{orange}{ ext{费马小定理}})，有：

[ ext{对于任意 }p mid a,a^{p-1}equiv1 pmod{p} ]

[color{red}{ ext{那么称 $g$ 为模 $p$ 的原根，当且仅当}{g^0,g^1,…,g^{p-2}} ext{对} p ext{的模数互不相同。}} ]

这样可能会更好理解些？至于为什么是上面那个绿绿的，我也不知道啊……

  • 补充：本原单位根

其实定义很简单，对于一些比较特殊的(n)次单位根，它的(0,1,…,n-1)次幂恰好生成了所有的(n)次单位根，我们把这样的单位根称为本原单位根。

在复数域 (C) 上，(omega_n=expfrac{2pi i}{n}=cosfrac{2pi}{n}+isinfrac{2pi}{n}) 是一个本原单位根。

（以上不是重点，下面这句话才是重点）

在有限域（即模质数）上，本原单位根和数论中的原根有关！

这也是我补充本原单位根的原因了……因为后面我们还会提到它。

并且我们可以证明的是，模质数的原根总是存在的（不知道为啥你就问度娘吧）。

并且原根的性质和单位根非常相似哦~

(NTT) 里原根的性质

换句话说，在模 (p) 意义下，(g) 可以被看做一个 (p-1) 次本原单位根。

反正感性理解就好，结论不就是用来记的嘛……

设 (nmid p-1) ,我们不妨令 (omega_n=g^{frac{p-1}{n}})（这个其实很像单位根的求法），然后你就可以得到一大串东西了：（好像也不多）

[omega_n^n=(g^{frac{p-1}{n}})^n=g^{p-1}equiv1pmod p ]

并且因为单位根的定义，你会知道 (omega_n^0,omega_n^0,…,omega_n^{n-1}) 互不相同……这不就完了嘛 $QAQ $

因此此时的 (omega_n) 在模 (p) 意义下具有 (n) 次本原单位根的性质，所以我们利用这个东西来定义 (DFT) 即可，这被称为 (color{red}{NTT})。

原根的求法

求模质数原根的办法： 将 (p-1) 分解质因数，即 (p-1=p_1^{a_1}p_2^{a2}…p_{n}^{a_n}) 为标准的(p-1)分解式，那么若 (g^{frac{p-1}{p_i}} e 1pmod p) ，则(g)是模 (p) 意义下的一个原根。

求模合数原根的办法： 其实你只需要把上面的 (p-1) 换成 (φ(p)) 即可。

证明：（假设法）
  先对于第一种模质数的情况：
  假设存在一个 (t<φ(p)=p-1) 使得 (g^tequiv1pmod p) 且 (forall i∈[1,n],g^{frac{p-1}{p_i}} e 1 pmod p)。由裴蜀定理可知：一定存在一组 ((a,b)) 使得：(a*t=(p-1)*b+gcd(t,p-1)) 成立。由欧拉定理/费马小定理：(g^{p-1}equiv 1 pmod p)。所以 (g^{a*t}equiv g^{(p-1)*b+gcd(t,p-1)}equiv1pmod p)。因为 (t<p-1)，所以又有 (gcd(t,p-1)<p-1)。又因为 (gcd(t,p-1)mid p-1)，于是 (exists i ∈[1,n],gcd(t,p-1)mid{g^{frac{p-1}{p_i}}})。那么就可以得到 (g^{frac{p-1}{p_i}}equiv g^{gcd(t, p-1)}equiv 1 pmod p)。
 因此假设不成立。所以原命题成立。
(Q.E.D.)

那么上述结论自然可以推广到模合数而把 (p-1) 改为 (φ(p)) 了。

以上的证明是针对于这一句话的：(δ_m(a)=φ(m))，也就是上述求法成立的条件。我们从另一个性质出发，即 (g^0,g^1,…,g^{p-1}) 模(p)的余数各不相同，那么可以得知这些余数的集合为：({1,2,…p-1})（不分顺序）。我们知道取模是满足乘法性质的。也就是说，如果(exists i∈[1,n])使得(g^{frac{p-1}{p^i}}=1)，我们知道唯有 (g^0 mod p=1 mod p = 1) ，那么很明显 (frac{p-1}{p^i}) 是不会为 (0) 的，那么一旦出现第一个数字的重复，也就是 (1)，那么很明显不满足原根的上述性质。

那么我就很好奇，为什么重复的(1)一定会出现在某个(g^{frac{p-1}{p^i}})当中呢？上面其实也已经证明，如果存在某个 (t<φ(p)) 使得 (g^tequiv1pmod p) 成立，那么必定会存在一个 (g^frac{p-1}{p^i}mod p=1) 。并且我认为这样可能更好理解。

反正证明这种东西感性理解一下就好了……

然后原根你也会求了，所以用原根 (g) 去代替 (FFT) 中的 (omega_n) (即单位根),就可以完成 (FFT) 了。

NTT 模数

首先根据 (FFT) 我们知道 (g^{frac{p-1}{n}})中的 (n) 是 (2) 的幂次，所以在做 (FFT) 的时候我们最好构造形如 (p=a*2^k+1)的质数，这样就可以满足上面的需求了。这样的质数最好取大一点（精度），所以这样的质数可以取：

(p=a*2^k+1) $ $ (a) (k) (g)
(3) (1) (1) (2)
(5) (1) (2) (2)
(17) (1) (4) (3)
(97) (3) (5) (5)
(193) (3) (6) (5)
(257) (1) (8) (3)
(7681) (15) (9) (17)
(12289) (3) (12) (11)
(40961) (5) (13) (3)
(65537) (1) (16) (3)
(786433) (3) (18) (10)
(5767169) (11) (19) (3)
(7340033) (7) (20) (3)
(23068673) (11) (21) (3)
(104857601) (25) (22) (3)
(167772161) (5) (25) (3)
(469762049) (7) (26) (3)
(color{red}{998244353}) (119) (23) (3)
(color{red}{1004535809}) (479) (21) (3)
(1998585857) (953) (21) (3)
(2013265921) (15) (27) (31)
(2281701377) (17) (27) (3)
(3221225473) (3) (30) (5)
(75161927681) (35) (31) (3)
(77309411329) (9) (33) (7)
(206158430209) (3) (36) (22)
(2061584302081) (15) (37) (7)
(2748779069441) (5) (39) (3)
(6597069766657) (3) (41) (5)
(39582418599937) (9) (42) (5)
(79164837199873) (9) (43) (5)
(263882790666241) (15) (44) (7)
(1231453023109121) (35) (45) (3)
(1337006139375617) (19) (46) (3)
(3799912185593857) (27) (47) (5)
(4222124650659841) (15) (48) (10)
(7881299347898369) (7) (50) (6)
(31525197391593473) (7) (52) (3)
(180143985094819841) (5) (55) (6)
(1945555039024054273) (27) (56) (5)
(4179340454199820289) (29) (57) (3)

理论到实现

第一件事，我们先来回忆一下 (FFT)，而 (FFT) 分为两个步骤：(DFT) 与 (IDFT)。

然后我们知道单位根的定义是对于复数域 (C)，有 (z^n=1) 的复数就是一个 (n)次单位根，(n) 次单位根有 (n) 个，为：(e^{frac{2pi k}{n}i}, k = {0,1,2…,n-1})（(i)是虚数单位）。（你用三角表示也是可以的，但是不好理解）

然后这里就用 (g^{frac{p-1}{n}}) 代替 上面那个就可以，因为你会发现他们具有相同的性质。

唯一具有卷积性质的变换就是 (DFT)，按照上面的式子，(DFT) 的变换式就是：

[X_k=sum_{n=0}^{N-1}x_ne^{-frac{2pi}{N}nki} k=0,1,2…,N-1 ]

(DFT) 反演（也就是逆变换）公式有：

[x_n=frac{1}{N}sum_{k=0}^{N-1}X_ke^{frac{2pi}{N}nki} n=0,1,2…N-1 ]

那么由于 (g^{frac{p-1}{n}}) 和 (e^{frac{2pi k}{n}i}, k = {0,1,2…,n-1}) 等价，所以：

[e^{-frac{2pi}{N}nki}equiv g^{frac{p-1}{N}}pmod p ]

接着我们就得到了快速数论变换的公式：

[X_k=sum_{n=0}^{N-1}x_ng^{nk}pmod p k=0,1,2…,N-1 ]

反演：

[x_n=frac{1}{N}sum_{k=0}^{N-1}X_kg^{-nk}pmod p n=0,1,2…N-1 ]

这样我们就成功把复数域转到了整数域，然后剩下的东西在(mod p)意义下考虑即可。上面也已经讲过素数的构造方法。

  • 补充：原根为何可以代替单位根？

考虑单位根的性质（也就是我们为什么使用单位根）：

  1. (n) 个单位根互不相同。那么(g^0,g^1,…,g^{p-2}) 在模 (p) 意义下也不相同，成立。

  2. (z^n=1)（(z∈C)）。根据费马小定理：(g^{p-1}equiv 1pmod p)，成立。

  3. 单位根对称分布。其实这对于原根也成立。

证明：
  (∵(g^{frac{p-1}{2}})^2equiv g^{p-1}equiv 1 pmod p)
  又(∵g^i e g^j(i e j ext{且}0≤i,j≤p-1))
  (∴g^{frac{p-1}{2}}= e 1, ext{则}g^{frac{p-1}{2}}=-1)
(Q.E.D.)

好了，我们真的可以使用 (NTT) 来完成 (FFT) 所能完成和所不能完成的事了。

可以尝试过一下模板题：(P3803)【模板】多项式乘法((FFT))

虽然题目没有要求模数，但是仍然可以取一个比较大的模数从而完成。

这里取 (p=998244353,g=3)。

我们结合代码来理解 (NTT) :

#include<bits/stdc++.h>
using namespace std;
const int N = (1 << 21) + 100, P = 998244353, G = 3;

inline int read() {
	static int x, f;
	static char c;
	while(!isdigit(c = getchar()) && c != '-');
	x = (f = c == '-') ? 0 : c - '0';
	while(isdigit(c = getchar())) 
		x = (x << 3) + (x << 1) + (c ^ 48);
	return f ? -x + 1 : x;
}

typedef long long ll;
inline ll power(ll a, ll k) {
	ll base = 1;
	for(; k; k >>= 1) {
		if(k & 1) base = (base * a) % P;
		a = (a * a) % P;
	}	
	return base;
}

int rev[1 << 21], limit = 1;
ll a[N], b[N];
inline void NTT(ll *f, int inv) {
	for(int i = 0; i < limit; i++)
		if(i < rev[i]) swap(f[i], f[rev[i]]);
	
	for(int len = 1; len < limit; len <<= 1) {
		ll temp = power(G, (P - 1) / (len << 1));
		if(inv) temp = power(temp, P - 2);
		//反NTT时，除法全部变为乘逆元 
		for(int l = 0; l < limit; l += len * 2) {
			ll w = 1;
			for(int k = 0; k < len; k++) {
				int x = f[l + k], y = w * f[l + k + len] % P;
				f[l + k] = (x + y) % P;
				f[l + k + len] = (x - y + P) % P;
				w = (w * temp) % P;
			}
		}
	} 
}

int main() {
	int n = read(), m = read(), L = -1;
	for(int i = 0; i <= n; i++) a[i] = (read() + P) % P;
	for(int i = 0; i <= m; i++) b[i] = (read() + P) % P;
	for( ; limit <= n + m; limit <<= 1, L++); 
	for(int i = 0; i < limit; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
	NTT(a, 0), NTT(b, 0);
	for(int i = 0; i < limit; i++)
		a[i] = (a[i] * b[i]) % P;
	NTT(a, 1);
	ll inv = power(limit, P - 2);
	for(int i = 0; i <= n + m; i++)
		printf("%lld ", (a[i] * inv) % P);
	
	return 0;
}

然后真的比(FFT)快好多呀……模板题评测记录：

( ext{递归}FFT)：点我

( ext{蝴蝶变换}FFT)：点我

(NTT)：点我

Upd20200918:文章当中还有一些错误，近期本人正在矫正。

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/Ning-H/p/12072626.html