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

FFT——快速傅里叶变换

什么是FFT

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

这个名字听起来很高级,实际上也很高级,它是DFTIDFT的加速版本,用于加速多项式乘法。

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

A(x)=anxn+an1xn1++a2x2+a1x1+a0

B(x)=bmxm+bm1xm1++b2x2+b1x1+b0

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

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

i=0n+m1(j=0iaibij)xi

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


前置知识

  • 复数

我们定义虚数 i, 且 i2=1,即 i=1,那么复数就形如 a+bii 是一个虚数单位,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 满足 wn=1,则称 wn 次单位根。

扔个单位圆:

为什么扔个单位圆呢?因为单位圆上任意一点到原点的距离,即模长都为1

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

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

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

则有nmπ=kπ,kN,即nm=k

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

  2. k>0,则k=1,2,3,,易知mn1个不重解

综上,复数wn次单位根有n个。

并且,复数的n次单位根n等分单位圆。

单位根的书写:

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

然后单位根还有个神奇的东西ωnk=ωnkmodn

单位根的性质

  1. 对于任意nωn0=1

  2. ωnk=(ωn1)k

  3. (ωnm)k=(ωnk)m

  4. ωniωnj=ωni+j

  5. ω2n2m=ωnm

  6. 2nωnm+n2=ωnm


多项式的表达方法

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

  2. 点值表达法

定义多项式F(x)

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

T=(k1,F(k1)),k2,F(k2),,(kn+1,F(n+1))

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

多项式P,Q分别取点(x,y1)(x,y2),则PQ可得到点(x,y1y2)。令S=PQ,则C(xi)=y1iy2i

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

  • (补充)乘法本质:卷积

根据我们上面暴力计算A(x)B(x),不妨设C=AB,那么:

C[k]=i=0kA[i]B[ki]

A[i],B[k1]分别为Axi的系数和Bxki的系数,那么我们知道:xixki=xk,故所有的xixki的系数对x的次项系数均有贡献,需要累加。

那么就可以写成:C[k]=i+j=kA[i]B[j]

那么形如C[k]=ij=kA[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(N2)的暴力,那我们要寻找一个更快的办法:

就是它!FFT算法

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


思想:

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

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

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

我们可以取3x相同的点:

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

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

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

神仙问题:tql but 这有用吗

当然有用啊:

算法流程:

系数表达式点值表达式系数表达式

这不就有个用了吗?

但这个思路仅仅是DFT+IDFTFFT是用来加速DFTIDFT的。

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

我们把多项式如下拆开:

F(x)=(a0+a2x2++an2xn2)+(a1x1+a3x3++an1xn1)

(保证n=2k,kN)

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

FL(x)=a0+a2x++an2xn21

FR(x)=a1+a3x++an1xn21

F(x)=FL(x2)+xFR(x2)

k<n2,把ωnk代入F(x)F(ωnk)=FL(ωn2k)+ωnkFR(ωn2k)

接着,就有F(ωnk)=FL(ωn/2k)+ωnkFR(ωn/2k)

那么,如果我们知道FL(x),FR(x)分别在ωn/20,ωn/21,,ωn/2n/21处的点值表示,那么我们就可以O(N)内求出F(x)ωn0,ωn1,,ωnn/21处的点值表示了。

接下来我们的目标就是求F(x)ωnn/2,ωnn/2+1,,ωnn1处的点值表示。

然后……

继续设k<n/2,然后把ωnk+n/2代入F(x)中。

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

F(ωnk+n/2)=FL(ωn/2k)ωnkFR(ωn/2k)

然后我们就可以求出F(x)ωn0,ωn1,,ωnn1处的点值表示,接着,我们就实现了DFT听着好累的样子……

但是我们还不知道FL(x),FR(x)ωn0,ωn1,,ωnn/21处的点值表示。难道要暴力代入吗?这样反而会使时间复杂度上升。

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

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

FFT分治的合并过程根据:

1. F(ωnk)=FL(ωn/2k)+ωnkFR(ωn/2k)

2. F(ωnk+n/2)=FL(ωn/2k)ωnkFR(ωn/2k)

从理论到实现

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


1. 预处理单位根

我们知道一个性质ωnk=(ωn1)k

那么,我们只要知道ωn1就可以求出ωn0,ωn1,ωn2,,ωnn1

怎么求ωn1

解三角形不就好了?已知|OA|=1,幅角是1n圆周,那么:

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

设幅角α=2πn,,则ωn1(cos(2πn),sin(2πn))

这边要牢记,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\n", 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数列,点值表示成MM[k]=i=0n1(ωnk)iF[i]

然后就有F[k]=1ni=0n1(ωnk)iM[i]

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

我们可以先求出和式,然后再除以n

接下来我们求出ωn0,ωn1,ωnn+1

并且仍有ωnj=(ωn1)j

所以我们还是求出ω01即可,并且,ω01ω01关于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并不是重点,而是里面用于优化的很重要的一个东西:单位根


  • 原根

但是如果 FFT 要求取模呢?

What?你不会?

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


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

原根的定义

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

δm(a)=φ(m)what is this?

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

p 为质数,则根据 费马小定理,有:

对于任意 pa,ap11(modp)

那么称 g 为模 p 的原根,当且仅当{g0,g1,,gp2} p 的模数互不相同。

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

  • 补充:本原单位根

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

在复数域 C 上,ωn=exp2πin=cos2πn+isin2πn 是一个本原单位根。

(以上不是重点,下面这句话才是重点)

在有限域(即模质数)上,本原单位根和数论中的原根有关!

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

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

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

NTT 里原根的性质

换句话说,在模 p 意义下,g 可以被看做一个 p1 次本原单位根。

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

np1 ,我们不妨令 ωn=gp1n(这个其实很像单位根的求法),然后你就可以得到一大串东西了:(好像也不多

ωnn=(gp1n)n=gp11(modp)

并且因为单位根的定义,你会知道 ωn0,ωn0,,ωnn1 互不相同……这不就完了嘛 QAQ

因此此时的 ωn 在模 p 意义下具有 n 次本原单位根的性质,所以我们利用这个东西来定义 DFT 即可,这被称为 NTT

原根的求法

求模质数原根的办法:p1 分解质因数,即 p1=p1a1p2a2pnan 为标准的p1分解式,那么若 gp1pi1(modp) ,则g是模 p 意义下的一个原根。

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

证明:(假设法)
  先对于第一种模质数的情况:
  假设存在一个 t<φ(p)=p1 使得 gt1(modp)i[1,n],gp1pi1(modp)。由裴蜀定理可知:一定存在一组 (a,b) 使得:at=(p1)b+gcd(t,p1) 成立。由欧拉定理/费马小定理:gp11(modp)。所以 gatg(p1)b+gcd(t,p1)1(modp)。因为 t<p1,所以又有 gcd(t,p1)<p1。又因为 gcd(t,p1)p1,于是 i[1,n],gcd(t,p1)gp1pi。那么就可以得到 gp1piggcd(t,p1)1(modp)
 因此假设不成立。所以原命题成立。
Q.E.D.

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

以上的证明是针对于这一句话的:δm(a)=φ(m),也就是上述求法成立的条件。我们从另一个性质出发,即 g0,g1,,gp1p的余数各不相同,那么可以得知这些余数的集合为:{1,2,p1}(不分顺序)。我们知道取模是满足乘法性质的。也就是说,如果i[1,n]使得gp1pi=1,我们知道唯有 g0modp=1modp=1 ,那么很明显 p1pi 是不会为 0 的,那么一旦出现第一个数字的重复,也就是 1,那么很明显不满足原根的上述性质。

那么我就很好奇,为什么重复的1一定会出现在某个gp1pi当中呢?上面其实也已经证明,如果存在某个 t<φ(p) 使得 gt1(modp) 成立,那么必定会存在一个 gp1pimodp=1 。并且我认为这样可能更好理解。

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

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

NTT 模数

首先根据 FFT 我们知道 gp1n中的 n2 的幂次,所以在做 FFT 的时候我们最好构造形如 p=a2k+1的质数,这样就可以满足上面的需求了。这样的质数最好取大一点(精度),所以这样的质数可以取:

p=a2k+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
998244353 119 23 3
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 分为两个步骤:DFTIDFT

然后我们知道单位根的定义是对于复数域 C,有 zn=1 的复数就是一个 n次单位根,n 次单位根有 n 个,为:e2πkni, k={0,1,2,n1}i是虚数单位)。(你用三角表示也是可以的,但是不好理解)

然后这里就用 gp1n 代替 上面那个就可以,因为你会发现他们具有相同的性质。

唯一具有卷积性质的变换就是 DFT,按照上面的式子,DFT 的变换式就是:

Xk=n=0N1xne2πNnki      k=0,1,2,N1

DFT 反演(也就是逆变换)公式有:

xn=1Nk=0N1Xke2πNnki      n=0,1,2N1

那么由于 gp1ne2πkni, k={0,1,2,n1} 等价,所以:

e2πNnkigp1N(modp)

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

Xk=n=0N1xngnk(modp)      k=0,1,2,N1

反演:

xn=1Nk=0N1Xkgnk(modp)      n=0,1,2N1

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

  • 补充:原根为何可以代替单位根?

考虑单位根的性质(也就是我们为什么使用单位根):

  1. n 个单位根互不相同。那么g0,g1,,gp2 在模 p 意义下也不相同,成立。

  2. zn=1zC)。根据费马小定理:gp11(modp),成立。

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

证明:
  (gp12)2gp11(modp)
  又gigj(ij0i,jp1)
  gp12=≠1,gp12=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快好多呀……模板题评测记录:

递归FFT点我

蝴蝶变换FFT点我

NTT点我


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


__EOF__

本文作者Ning-H
本文链接https://www.cnblogs.com/Ning-H/p/12072626.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本文版权归作者和博客园共有,欢迎转载,但未经作者同意必须在文章页面给出原文连接,否则保留追究法律责任的权利。
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   Ning-H  阅读(1862)  评论(0编辑  收藏  举报
编辑推荐:
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
阅读排行:
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】
点击右上角即可分享
微信分享提示