From FFT to NTT

ACM 数论
ACM 数论
发布日期:   2022-08-22
更新日期:   2022-08-31
文章字数:   5.6k
阅读时长:   26 分
阅读次数:  

1、前言

FFT全称（Fast Fourier Transformation）快速傅里叶变换

NTT全称（Fast Number Theoretical Transformation）快速数论变换

1.1、引入问题

如何计算多项式$f(x)$与$g(x)$的乘积？

例如：计算$f(x)=1+x^2+x^3$与$g(x)=1+x^5+x^6+x^7$的乘积

通常而言，我们思考暴力的算法：

显而易见，这个算法的时间复杂度是$O(n^2)$

那么有快速的计算方法吗？

当然是有的，于1965年由J.W.库利和T.W.图基提出FFT算法，可以在$O(nlogn)$的时间复杂度下计算出答案。

1.2、定义一些记号：

# $f(x)$表示一个多项式

#对于一个$n$次多项式$f(x)$的第$m$项，记为$f[m]$;换言之，$f(x)=\sum_{i=0}^{n}f[i]x^i$

1.3、多项式乘法理解

例如：现有两个$n$次多项式$f(x)$与$g(x)$,需要计算这两个多项式的乘积.

不妨设$h(x)=f(x)\times g(x)$

那么通过1.1中的例子，我们可以发现，其实$h[k]=\sum_{i=0}^{k}f[i]g[k-i]$

换言之，我们可以理解为：要想乘出$x^k$,需要$f(x)$的$x^i$项和$g(x)$的$x^{k-i}$项相乘，然后将他们的相乘的贡献累加便是$h[k]$

1.4、卷积

将形如:$h[k]=\sum_{i⊕j}f[i]g[j]$的式子称为卷积，其中⊕在数学中称为假设符号（假设什么运算）

例如，在多项式乘法中，为加法卷积：$h[k]=\sum_{i+j}f[i]g[j]$

因此，对于卷积可以利用FFT构造多项式进行计算！

2、FFT思想

2.1、多项式的点值表达式

我们可以从

得出推论：在平面直角坐标系中，$n+1$个点就能唯一确定一个$n$次多项式。

因此，对于一个$n$次多项式$f(x)$而言，我们可以用$n+1$个点来描述一个多项式；换言之，两个多项式相乘，等价于两个多项式的对应点值相乘。(点值的乘法对应着整个多项式的乘法,也就是浓缩了整个多项式)

将能唯一表达一个$n$次多项式$f(x)$的$n+1$个点,称为$f(x)$的点值表达式。

由于$n$次多项式$f(x)$与$n$次多项式$g(x)$,他们的乘积一定是一个$2n+1$次的多项式，所以我们只需要对$f(x)$和$g(x)$找到$2n+1$个对应的点，然后将他们对应点的纵坐标相乘，便可以确定多项式$f(x)\times g(x)$的点值表达式；

最后我们将多项式$f(x)\times g(x)$的点值表达式转化为多项式的一般表达式（系数表达式）便得到了我们想要的多项式乘积了。

而这也是FFT算法的流程

2.2、DFT&IDFT

DFT：把多项式的系数表达式转化为点值表达式的变换过程

IDFT：把多项式的点值表达式转化为系数表达式的变换过程

下面这张图便对比了FFT算法与暴力算法的区别

3、FFT的实现

3.1、DFT

​ 由于我们想要多项式的点值表达式，所以我们需要利用一些方法获得多项式的上的点；但是细想一下会发现，

​ 这样一来，我们通过随便找$n+1$个点的方式，会让我们的时间复杂度和暴力没什么区别！

​ 那么有快速找到$n+1$个点的方法么？答案当然是有的！

​ 我们可以这样考虑问题：

3.1.1、Good idea 1：函数的奇偶性

​ 首先，我们不考虑一般多项式，只考虑一个简单的多项式$f(x)=x^2$

那么，如果我们要从$f(x)$中选择8个点，我们如何选择？

这里我们需要注意$f(x)$的一个优美的性质————$f(x)$是偶函数，所以我们计算出一个点$(x,f(x))$后，便可以立马知道点$(-x,f(x))$也一定是$y(x)$上的点。

那么，如果$f(x)=x^3$,应该如何处理喃？

显而易见，我们会发现$f(x)=x^3$是一个奇函数，它与偶函数的区别在于，奇函数上的点是关于原点对称的，换言之，如果我们计算出点$(x,f(x))$,那就可以知道点$(-x,-f(x))$也是多项式上的点！

既然这样，对于多项式$f(x)=1+2x^2+3x^3+4x^4+5x^5$

那么我可以将$f(x)$中按幂数分离，因此可得$f(x)=(1+2x^2+4x^4)+(3x^3+5x^5)$

这样，对于$f_o(x^2),f_e(x^2)$都是关于$x^2$的多项式，我们便可以利用函数的奇偶性快速找到6个点，然后利用点值表达式的计算方式计算出该函数的点值表达式。

然后我们可以推广到一般多项式$f(x)=a_0+a_1x+a_2x^2+…..a_nx^{n}$

而对于$f_e(x)与f_o(x)$,我们可以重复$f(x)$的奇偶分离的方式；这样我们可以把对于$n$次多项式转化为求$\frac{n}{2}$次多项式（类似分治的思想），所以这个思路下，求多项式点值表达式的时间复杂度是$O(nlogn)$

因此，只要我们最后找到的点为合法的点，那么算法就是合法的！

但是这里出现了一个问题：对于$f(x)$中的$[\pm x_1,\pm x_2…\pm x\frac{n}{2}]$而言，每一个$\pm x_i$都是$\pm$配对的；但是对于$f_o(x^2)与f_e(x^2)$中的$[x_1^2,x_2^2..x_n^2]$而言,没有一个$x_i^2$是配对的;所以在实数域下，是找不到满足条件的点的！

不过解决问题的的方法总比问题多，所以我们考虑把数域扩展到复数域！

3.1.2、Good idea 2：复数域下的单位根

3.1.2.1解决疑惑：为何复数域下，满足3.1.1中的分治找点

一个简单例子：求多项式$f(x)=1+2x+3x^3$的点值表达式

3.1.2.2引入单位根

单位根：在复平面上，以单位圆点为起点，单位圆的 $n$ 等分点为终点，在单位圆上可以得到 $n$ 个复数，设幅角为正且最小的复数为$n$次单位根，记作$\omega_n$

#推导$\omega_n$

#推导单位根的一些性质:

(1)$w_n^{a}\times w_n^{b}=w_n^{i+j}$

(2)$w_{rn}^{rk}=w_n^{k}$

(3)$wn^{2k}=w{\frac{n}{2}}^{k}$

(4)如果$n$是偶数, $w_n^{k+\frac{n}{2}}=-w_n^{k}$

3.1.2.3代入单位根

现已知$f(x)=f_e(x^2)+xf_o(x^2)$

(1)代入$x=w_n^{k},k<\frac{n}{2}$

(2)代入$x=w_n^{k+\frac{n}{2}},k<\frac{n}{2}$

既然这样，那么我们再3.1.1中$O(nlogn)$求多项式的思路是可行的！！！！

3.1.3：优化分治（蝴蝶变换）

我们用下面一张图来表示3.1.1中的分治过程：

显然,对于上述过程是可以使用递归的写法,但是递归的写法常数过大，

但是我们需要观察每个数最终的位置,最终发现对于幂数$i$而言，其实$i$最终在序列中的位置其实是$i$在二进制下翻转（蝴蝶变换）得到的数

因此我们可以把每个数都放在其最终位置上，不断往上迭代即可！

#附一个蝴蝶变换的小结论：

蝴蝶变换过后的序列，任取三个数都不构成等差数列。

3.2、IDFT

现在我们已经知道$k$次多项式$f(x)$的点值表达式：${(w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),(w_n^2,f(w_n^2)),….,(w_n^n,f(w_n^n))} (k\leq n,n=2^m-1)$

那么可一把点值表达式改写成与之对应的矩阵的形式：

那么在这个矩阵乘法中：

最后我们注意到对于$w_n^{n-j}$，其实它为$w_n^j$的共轭复数；因此在IDFT中，只需将DFT中的乘以单位根变为乘以其共轭复数即可，但是分治结束后，需要除以$n$

其实上述需要求逆的矩阵是一个范德蒙德矩阵，也可以利用范德蒙德矩阵求逆的过程

好了！耗时三天终于敲完FFT的公式推导了！！！下面的代码实现就比较容易了！

4、FFT代码实现

直接上模版：

double PI = acos(-1.0);
//复数类
struct Complex
{
    double r, i;
    Complex() {}
    Complex(double r, double i) : r(r), i(i) {}
    inline void real(const double &x) { r = x; }
    inline double real() { return r; }
    inline Complex operator+(const Complex &rhs) const
    {
        return Complex(r + rhs.r, i + rhs.i);
    }
    inline Complex operator-(const Complex &rhs) const
    {
        return Complex(r - rhs.r, i - rhs.i);
    }
    inline Complex operator*(const Complex &rhs) const
    {
        return Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
    }
    inline Complex operator*(const double &rhs) const
    {
        return Complex(r * rhs, i * rhs);
    }
    inline void operator/=(const double &x)
    {
        r /= x, i /= x;
    }
    inline void operator*=(const Complex &rhs)
    {
        *this = Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
    }
    inline void operator+=(const Complex &rhs)
    {
        r += rhs.r, i += rhs.i;
    }
    inline Complex conj()
    {
        return Complex(r, -i);
    }
};
// i的二进制翻转数组
int bitrev[maxn];
// FFT预处理
void fft_prepare(int n)
{
    int bits = 0;
    //由于FFT需要项数为2的整数次方倍，所以多项式（形式）次数为第一个大于n的二的正整数次方，bits表示多少位
    while ((1 << bits) < n)
    {
        bits++;
    }
    for (int i = 0; i < n; i++)
    {
        //递推式：二进制翻转
        // 0 1 2 3 4 5 6 7
        // 0|4|2|6|1|5|3|7
        bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (bits - 1));
    }
}

//快速傅里叶变换系数表达式转换为点值表达式
void fft(Complex *a, int n, int type = 1)
{
    //求出要迭代的序列
    for (int i = 0; i < n; i++)
    {
        if (i < bitrev[i])
        {
            swap(a[i], a[bitrev[i]]);
        }
    }
    //待合并区间的中点
    for (int mid = 1; mid < n; mid <<= 1)
    {
        //单位根
        Complex Wn{cos(PI / mid), type * sin(PI / mid)};
        // R是区间的右端点，j表示已经到哪个位置
        for (int R = mid << 1, j = 0; j < n; j += R)
        {
            Complex w{1, 0}; //幂
            //枚举左半部分
            for (int k = 0; k < mid; k++, w *= Wn)
            {
                //蝴蝶效应
                Complex x = a[j + k], y = w * a[j + mid + k];
                a[j + k] = x + y;
                a[j + k + mid] = x - y;
            }
        }
    }
    // type为-1则进行逆变换，即IDFT
    if (type == -1)
    {
        for (int i = 0; i < n; i++)
        {
            a[i].r /= n;
        }
    }
}

代码有注释，就不细说了！

5、FFT的瓶颈

虽然FFT的在时间复杂度上非常优秀，但是计算机在计算时仍然会遇到一些问题！

例如：由于复数类中使用的是double类型，而且在计算单位根存在除法，除此之外计算机计算下$\pi$始终是一个逼近值，因此必然也必然会存在误差！所以在特殊情况下计算机计算会出现一定的误差！

但是有没有解决办法喃？

可以说在一定情况下是可以避免误差的！这个时候就要引出原根

6、To NTT

原根：

通过定义我们可以知道，原根是一种定义在模意义下的定义，因此对于FFT遇到的的精度问题，我们也只能保证在模意义下的答案正确

我们通过推导可以发现原根具有和单位根几乎一样的性质，因此我们可以通过原根推出基本相同的结论

这样一来我们就不必考虑复数域了！！

也可以理解为NTT其实就是把FFT中的单位根换成了原根

6.1 NTT代码实现

const int g = 3, gi = 332748118, mod = 998244353;
// 998244353的一个原根为3且998244353-1=2^23*119，3在模998244353意义下的逆元为332748118
int qp(int a, int n)
{
    int ans = 1;
    while (n)
    {
        if (n & 1)
        {
            ans = ans * a % mod;
        }
        a = a * a % mod;
        n >>= 1;
    }
    return ans;
}
// i的二进制翻转数组
int bitrev[maxn];
// NTT预处理
void ntt_prepare(int n)
{
    int bits = 0;
    //由于NTT需要项数为2的整数次方倍，所以多项式（形式）次数为第一个大于n的二的正整数次方，bits表示多少位
    while ((1 << bits) < n)
    {
        bits++;
    }
    for (int i = 0; i < n; i++)
    {
        //递推式：二进制翻转
        // 0 1 2 3 4 5 6 7
        // 0|4|2|6|1|5|3|7
        bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (bits - 1));
    }
}
// NTT，type=1时系数表示法转点值表示法，否则点值表示法转系数表示法
void NTT(int *a, int n, int type = 1)
{
    for (int i = 0; i < n; i++)
    {
        if (i < bitrev[i])
        {
            swap(a[i], a[bitrev[i]]);
        }
    }
    //在这之前NTT与FFT无异
    for (int mid = 1; mid < n; mid <<= 1)
    {
        //单位原根g_n
        int gn = qp(type == 1 ? g : gi, (mod - 1) / (mid << 1));
        //枚举具体区间，j也就是区间右端点
        for (int R = mid << 1, j = 0; j < n; j += R)
        {
            int g0 = 1;
            //合并，记得要取模
            for (int k = 0; k < mid; k++, g0 = gn * g0 % mod)
            {
                int x = a[j + k], y = g0 * a[mid + j + k] % mod;
                a[j + k] = (x + y) % mod;
                a[j + k + mid] = (x - y + mod) % mod;
            }
        }
    }
}

可以发现，NTT和FFT代码大体逻辑是相同的！！！只是从单位根换成了原根

6.2 tips

在抽象代数中，原根就是循环群的生成元。这个概念只在模 $m$缩剩余系关于乘法形成的群中有“原根”这个名字，在一般的循环群中都称作“生成元”。

并非每个模$m$缩剩余系关于乘法形成的群都是循环群，存在原根就表明它同构于循环群，如果不存在原根就表明不同构。

文章作者: fatzard
文章链接: http://wjknowledge.top/2022/08/22/from-fft-to-ntt/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 fatzard !
ACM 数论
赏
评论
  目录
本站总访问量次