FFT——快速傅里叶变换
什么是FFT
F F T F F T 全称(F a s t F o u r i e r T r a n s f o r m a t i o n F a s t F o u r i e r T r a n s f o r m a t i o n ),中文名:快速傅里叶离散变换。
这个名字听起来很高级,实际上也很高级,它是D F T D F T 和I D F T I D F T 的加速版本,用于加速多项式乘法。
接下来我们把某个多项式记为“大写字母+元素”,如多项式:
A ( x ) = a n x n + a n − 1 x n − 1 + … + a 2 x 2 + a 1 x 1 + a 0 A ( x ) = a n x n + a n − 1 x n − 1 + … + a 2 x 2 + a 1 x 1 + a 0
B ( x ) = b m x m + b m − 1 x m − 1 + … + b 2 x 2 + b 1 x 1 + b 0 B ( x ) = b m x m + b m − 1 x m − 1 + … + b 2 x 2 + b 1 x 1 + b 0
你也可以理解为A / B A / B 是关于x x 的函数
现在,我们要求A ( x ) ∗ B ( x ) A ( x ) ∗ B ( x ) ,即两个多项式的乘积(卷积),那么我们很容易想到的就是O ( N 2 ) O ( N 2 ) 的做法:把每一项依次分别相乘再相加。
n + m − 1 ∑ i = 0 ( i ∑ j = 0 a i ∗ b i − j ) ∗ x i ∑ i = 0 n + m − 1 ( ∑ j = 0 i a i ∗ b i − j ) ∗ x i
这样的时间复杂度明显不尽人意,那么F F T F F T 就是把这种算式计算优化到N l o g N N l o g N 的算法。
前置知识
我们定义虚数 i i , 且 i 2 = − 1 i 2 = − 1 ,即 i = √ − 1 i = − 1 ,那么复数就形如 a + b i a + b i ,i i 是一个虚数单位,a a 叫做复数的实部, b b 叫做复数的虚部,那么 a + b i a + b i 在复平面(高中会讲)可以类似于平面直角坐标系,表示为点 ( a , b ) ( a , b ) 。
复数的运算法则
加法:实部和虚部分别相加。
减法:实部和虚部分别相减。
乘法:像一次多项式一样相乘。
除法:分子分母同乘分母的共轭。
共轭 :两个复数实部相同,虚部相反。
上面除法分子分母同乘分母的共轭,把分母有理化成实数。
这边大家直接记一个重要结论:
复数相乘,辐角相加,模长相乘。
单位根
对于复数w w ,若整数 n n 满足 w n = 1 w n = 1 ,则称 w w 是 n n 次单位根。
扔个单位圆:
为什么扔个单位圆呢?因为单位圆上任意一点到原点的距离,即模长都为1 1 。
那么,对于一个复数,只有模长为1 1 的复数才有可能成为n n 次单位根!
现在,我们在复平面上把复数锁定在了这个单位圆上。
设复数w w 模长为q q ,其幅角为m π m π ,m m 为[ 0 , 2 ] [ 0 , 2 ] 间实数,那么,其n n 次方幅角为n m π n m π 。
则有n m π = k π , k ∈ N n m π = k π , k ∈ N ,即n m = k n m = k 。
若k = 0 k = 0 ,因为n > 0 n > 0 ,所以方程有且只有m = 0 m = 0 一解。
若k > 0 k > 0 ,则k = 1 , 2 , 3 , … k = 1 , 2 , 3 , … ,易知m m 有n − 1 n − 1 个不重解
综上,复数w w 的n n 次单位根有n n 个。
并且,复数的n n 次单位根n n 等分单位圆。
单位根的书写:
单位根用符号ω ω 表示,从1开始,沿单位圆逆时针把单位根从0 0 开始标号。
然后单位根还有个神奇的东西ω k n = ω k mod n n ω n k = ω n k mod n
单位根的性质
对于任意n n ,ω 0 n = 1 ω n 0 = 1 。
ω k n = ( ω 1 n ) k ω n k = ( ω n 1 ) k 。
( ω m n ) k = ( ω k n ) m ( ω n m ) k = ( ω n k ) m
ω i n ∗ ω j n = ω i + j n ω n i ∗ ω n j = ω n i + j 。
ω 2 m 2 n = ω m n ω 2 n 2 m = ω n m 。
若 2 ∣ n 2 ∣ n ,ω m + n 2 n = − ω m n ω n m + n 2 = − ω n m
多项式的表达方法
系数表达法:这个如同我们上面写的多项式即为系数表达法。
点值表达法
定义多项式F ( x ) F ( x ) :
我们知道,1 1 个n元一次方程最少需要n n 个式子求解,而我们把系数表达F ( x ) F ( x ) 抽象为一个函数F ( x ) F ( x ) ,那么这个函数至少需要n + 1 n + 1 个点确定下来,那么我们任取n + 1 n + 1 个不同的数构成集合U = k 1 , k 2 , … , k n + 1 U = k 1 , k 2 , … , k n + 1 ,对F ( x ) F ( x ) 分别求值得到一个新的集合,那么将两个集合合并成点集,有:
T = ( k 1 , F ( k 1 ) ) , k 2 , F ( k 2 ) , … , ( k n + 1 , F ( n + 1 ) ) T = ( k 1 , F ( k 1 ) ) , k 2 , F ( k 2 ) , … , ( k n + 1 , F ( n + 1 ) )
这便是F ( x ) F ( x ) 的点值表达式,并且,点值表达法下的乘法运算更加简单:
多项式P , Q P , Q 分别取点( x , y 1 ) ( x , y 1 ) 和( x , y 2 ) ( x , y 2 ) ,则P ∗ Q P ∗ Q 可得到点( x , y 1 ∗ y 2 ) ( x , y 1 ∗ y 2 ) 。令S = P ∗ Q S = P ∗ Q ,则C ( x i ) = y 1 i ∗ y 2 i C ( x i ) = y 1 i ∗ y 2 i 。
可见,点值表达法下,多项式乘法是O ( N ) O ( N ) 的。
根据我们上面暴力计算A ( x ) ∗ B ( x ) A ( x ) ∗ B ( x ) ,不妨设C = A ∗ B C = A ∗ B ,那么:
C [ k ] = ∑ k i = 0 A [ i ] B [ k − i ] C [ k ] = ∑ i = 0 k A [ i ] B [ k − i ]
A [ i ] , B [ k − 1 ] A [ i ] , B [ k − 1 ] 分别为A A 中x i x i 的系数和B B 中x k − i x k − i 的系数,那么我们知道:x i ∗ x k − i = x k x i ∗ x k − i = x k ,故所有的x i ∗ x k − i x i ∗ x k − i 的系数对x x 的次项系数均有贡献,需要累加。
那么就可以写成:C [ k ] = ∑ i + j = k A [ i ] B [ j ] C [ k ] = ∑ i + j = k A [ i ] B [ j ]
那么形如C [ k ] = ∑ i ⊕ j = k A [ i ] B [ j ] C [ k ] = ∑ i ⊕ j = k A [ i ] B [ j ] (⊕ ⊕ 为某种运算)的式子称为卷积。
那么,多项式的乘法就是对于加法的卷积。
我们仍然可以对多项式相乘的模板打个暴力:
模板题亲测 44 分。
算法理论
上面的主要部分还是写了一个O ( N 2 ) O ( N 2 ) 的暴力,那我们要寻找一个更快的办法:
就是它!F F T F F T 算法
前面我们已经说过:F F T F F T 算法是一个O ( N l o g N ) O ( N l o g N ) 的算法,但常数较大(相信你们会卡常),并且,F F T F F T 是D F T D F T 和I D F T I D F T 的快速算法:
思想:
首先你要会点值表达式,当然我们已经讲过。
然后我们来看A ( x ) = x 2 + 3 x − 2 A ( x ) = x 2 + 3 x − 2 (红色)与B ( x ) = − x 2 + 6 x + 4 B ( x ) = − x 2 + 6 x + 4 (绿色)的取点得到点的点值表达式:
在A ( x ) A ( x ) 上取三个点(紫色)( − 4 , 2 ) ( − 2 , − 4 ) ( 1 , 2 ) ( − 4 , 2 ) ( − 2 , − 4 ) ( 1 , 2 ) ;在B ( x ) B ( x ) 上取三个点( 0 , 4 ) ( 2 , 12 ) ( 6 , 4 ) ( 0 , 4 ) ( 2 , 12 ) ( 6 , 4 ) ,然后转为系数表达O ( N 2 ) O ( N 2 ) 相乘然后再带入x x 得到C C 的点值表达式?这时绝对不行的,我们上面说过,点值表达式相乘可以在O ( n ) O ( n ) 内解决,那怎么办呢?
我们可以取3 3 个x x 相同的点:
这样就行了?那是不可能的。我们知道,两个二次项式相乘得到的是一个四次项式,而C ( x ) = y 1 ∗ y 2 C ( x ) = y 1 ∗ y 2 ,只能得到三个点,这肯定是不行的。
这很明显我们需要2 n + 1 2 n + 1 个点,不够怎么办呢?再添两个不就好了?
下面添加了I , J I , J (横坐标相同),G , H G , H (横坐标相同),这样就行了。
那么,我们只需要进行2 n + 1 2 n + 1 次乘法即可,省略常数(都说了常数巨大了),复杂度O ( N ) O ( N )
神仙问题:t q l t q l b u t b u t 这有用吗
当然有用啊:
算法流程:
系 数 表 达 式 ⟹ 点 值 表 达 式 ⟹ 系 数 表 达 式 系数表达式 ⟹ 点值表达式 ⟹ 系数表达式
这不就有个用了吗?
但这个思路仅仅是D F T + I D F T D F T + I D F T ,而F F T F F T 是用来加速D F T D F T 和I D F T I D F T 的。
如果让我们求点值表达式,那我们肯定会带整数进去,因为这样好计算。但是傅里叶 ,(不知道是不是脑抽了 ),把复数带进了多项式中,然后神奇的事情就发生了……
我们把多项式如下拆开:
F ( x ) = ( a 0 + a 2 x 2 + … + a n − 2 x n − 2 ) + ( a 1 x 1 + a 3 x 3 + … + a n − 1 x n − 1 ) F ( x ) = ( a 0 + a 2 x 2 + … + a n − 2 x n − 2 ) + ( a 1 x 1 + a 3 x 3 + … + a n − 1 x n − 1 )
(保证n = 2 k , k ∈ N ∗ n = 2 k , k ∈ N ∗ )
然后设两个有n / 2 n / 2 次项的多项式F L ( x ) F L ( x ) 和F R ( x ) F R ( x )
F L ( x ) = a 0 + a 2 x + … + a n − 2 x n 2 − 1 F L ( x ) = a 0 + a 2 x + … + a n − 2 x n 2 − 1
F R ( x ) = a 1 + a 3 x + … + a n − 1 x n 2 − 1 F R ( x ) = a 1 + a 3 x + … + a n − 1 x n 2 − 1
F ( x ) = F L ( x 2 ) + x F R ( x 2 ) F ( x ) = F L ( x 2 ) + x F R ( x 2 )
设k < n 2 k < n 2 ,把ω k n ω n k 代入F ( x ) F ( x ) :F ( ω k n ) = F L ( ω 2 k n ) + ω k n F R ( ω 2 k n ) F ( ω n k ) = F L ( ω n 2 k ) + ω n k F R ( ω n 2 k )
接着,就有F ( ω k n ) = F L ( ω k n / 2 ) + ω k n F R ( ω k n / 2 ) F ( ω n k ) = F L ( ω n / 2 k ) + ω n k F R ( ω n / 2 k )
那么,如果我们知道F L ( x ) , F R ( x ) F L ( x ) , F R ( x ) 分别在ω 0 n / 2 , ω 1 n / 2 , … , ω n / 2 − 1 n / 2 ω n / 2 0 , ω n / 2 1 , … , ω n / 2 n / 2 − 1 处的点值表示,那么我们就可以O ( N ) O ( N ) 内求出F ( x ) F ( x ) 在ω 0 n , ω 1 n , … , ω n / 2 − 1 n ω n 0 , ω n 1 , … , ω n n / 2 − 1 处的点值表示了。
接下来我们的目标就是求F ( x ) F ( x ) 在ω n / 2 n , ω n / 2 + 1 n , … , ω n − 1 n ω n n / 2 , ω n n / 2 + 1 , … , ω n n − 1 处的点值表示。
然后……
继续设k < n / 2 k < n / 2 ,然后把ω k + n / 2 n ω n k + n / 2 代入F ( x ) F ( x ) 中。
运用一系列单位根的性质,读者可以自行证明出:
F ( ω k + n / 2 n ) = F L ( ω k n / 2 ) − ω k n F R ( ω k n / 2 ) F ( ω n k + n / 2 ) = F L ( ω n / 2 k ) − ω n k F R ( ω n / 2 k )
然后我们就可以求出F ( x ) F ( x ) 在ω 0 n , ω 1 n , … , ω n − 1 n ω n 0 , ω n 1 , … , ω n n − 1 处的点值表示,接着,我们就实现了D F T D F T (听着好累的样子…… )
但是我们还不知道F L ( x ) , F R ( x ) F L ( x ) , F R ( x ) 在ω 0 n , ω 1 n , … , ω n / 2 − 1 n ω n 0 , ω n 1 , … , ω n n / 2 − 1 处的点值表示。难道要暴力代入吗?这样反而会使时间复杂度上升。
但是我们想想:F L ( x ) F L ( x ) 和F R ( x ) F R ( x ) 是由原来问题F ( x ) F ( x ) 转化来的,怎么转化的?看到这里,你应该想到了分治 。我们知道,分治的时间复杂度是l o g N l o g N ,好像完成了一半了?
然后,因为F L ( x ) , F R ( x ) F L ( x ) , F R ( x ) 是F ( x ) F ( x ) 的子问题,那么我们可以对它们继续进行分解成一个个小的子问题,然后再合并,类似于归并排序。
F F T F F T 分治的合并过程根据:
1. F ( ω k n ) = F L ( ω k n / 2 ) + ω k n F R ( ω k n / 2 ) F ( ω n k ) = F L ( ω n / 2 k ) + ω n k F R ( ω n / 2 k )
2. F ( ω k + n / 2 n ) = F L ( ω k n / 2 ) − ω k n F R ( ω k n / 2 ) F ( ω n k + n / 2 ) = F L ( ω n / 2 k ) − ω n k F R ( ω n / 2 k )
从理论到实现
上文的所有证明都基于2 = 2 k , k ∈ N ∗ 2 = 2 k , k ∈ N ∗ ,但是我们做题的时候并没有这一点,所以我们可以补一些系数为0 0 的项,这对计算结果没有影响。
1. 预处理单位根
我们知道一个性质 ω k n = ( ω 1 n ) k ω n k = ( ω n 1 ) k
那么,我们只要知道ω 1 n ω n 1 就可以求出ω 0 n , ω 1 n , ω 2 n , … , ω n − 1 n ω n 0 , ω n 1 , ω n 2 , … , ω n n − 1 。
怎么求ω 1 n ω n 1 ?
解三角形不就好了?已知| O A | = 1 | O A | = 1 ,幅角是1 n 1 n 圆周,那么:
(由于c + + c + + 下三角函数用弧度制表示,以下全部使用弧度制。)
设幅角α = 2 π n α = 2 π n ,,则ω 1 n ( c o s ( 2 π n ) , s i n ( 2 π n ) ) ω n 1 ( c o s ( 2 π n ) , s i n ( 2 π n ) ) 。
这边要牢记,x x 是实轴,y y 是虚轴:
然后我们就可以实现D F T D F T 了
因为还没讲I D F T I D F T ,先写D F T D F T :
递归D F T D F T 实现(复数板子上面有了下面就不写了)
事实上虚数可以用(小心TLE ):
然后我们得到了一堆点值表达……
这并没有什么用
于是,我们现在来学习I D F T I D F T ——D F T D F T 的逆运算。
刚刚的多项式为F ( x ) F ( x ) ,系数是A A 数列,点值表示成M M ,M [ k ] = n − 1 ∑ i = 0 ( ω k n ) i ∗ F [ i ] M [ k ] = ∑ i = 0 n − 1 ( ω n k ) i ∗ F [ i ]
然后就有F [ k ] = 1 n n − 1 ∑ i = 0 ( ω − k n ) i ∗ M [ i ] F [ k ] = 1 n ∑ i = 0 n − 1 ( ω n − k ) i ∗ M [ i ]
结论记住就好了,可以自己尝试证明。
我们可以先求出和式,然后再除以n n 。
接下来我们求出ω 0 n , ω − 1 n , … ω − n + 1 n ω n 0 , ω n − 1 , … ω n − n + 1
并且仍有 :ω − j n = ( ω − 1 n ) j ω n − j = ( ω n − 1 ) j
所以我们还是求出ω − 1 0 ω 0 − 1 即可,并且,ω − 1 0 ω 0 − 1 与ω 1 0 ω 0 1 关于x x 轴对称。
接着TLE了……(丝毫没有卡常的欲望 )
因为递归常数巨大,并且很容易爆空间。(但实际上是因为写了l o n g d o u b l e l o n g d o u b l e 炸了)
下面是AC的递归代码,以后不要写long double 了
实际上跑的很快的递归评测记录 。
那有没有更优化的办法呢?这是肯定的。
蝴蝶变换
先看一张图
这张图是什么?这不是我们分治的顺序吗?
我在一开始和结束都打上了二进制,你有没有什么发现?
是的,每一个数的二进制都反过来了。
然后我们的任务就转换成求二进制反序:
结论牢记即可。
还记得D P D P 吗?(貌似没什么关系 )
我们用循环实现F F T F F T ,也就像D P D P 一样定义状态、阶段决策:
第一层循环枚举变换 区间大小(阶段),第二层循环枚举开头(决策),第三层处理区间。
评测记录
NTT——快速数论变换
首先,学习N T T N T T 算法你要学会F F T F F T (理解也可)(刚刚我们才讲完你不会忘了吧……)(我认为我写的很详细了呢Q A Q Q A Q )。
接下来,你要确保你会原根这个东西(起码要知道),不过我在后面还会再提,具体可以先看我的这一篇博客:同余相关 ,在比较下面的部分(由于还待完善,所以只是写了一点)。
然后就是希望大家看了这一篇博客对N T T N T T 有所了解,可能我讲的不好,但是起码你要知道N T T N T T 是什么。
前置知识
(怎么又是前置知识 )
F F T F F T
上面已经说了先理解F F T F F T ,并且我也给了写的比较详细的博客,这里不再展开。
但其实F F T F F T 并不是重点,而是里面用于优化的很重要的一个东西:单 位 根 单位根 。
但是如果 F F T F F T 要求取模呢?
W h a t ? 你 不 会 ? W h a t ? 你不会?
好像因为 F F T F F T 用的单位根不太好取模(复数),所以我们可能需要找一个东西来替代下单位根,那么就要求这个东西和单位根具有类似的性质。然后科学家们就开始寻找了……(也不知道是谁找到的)接下来他们找到了原根 。
首先我们要了解一下啥是原根。(如果你看了上面给出的博客你应该知道了)
原根的定义
其实在那一篇博客里的原根定义很容易让人摸不着头脑……为啥呢?
δ m ( a ) = φ ( m ) ? w h a t i s t h i s ? δ m ( a ) = φ ( m ) ? w h a t i s t h i s ?
因为它很不好求……那我们不妨换一种说法:
设 p p 为质数,则根据 费 马 小 定 理 费马小定理 ,有:
对 于 任 意 p ∤ a , a p − 1 ≡ 1 ( mod p ) 对于任意 p ∤ a , a p − 1 ≡ 1 ( mod p )
那 么 称 g 为 模 p 的 原 根 , 当 且 仅 当 { g 0 , g 1 , … , g p − 2 } 对 p 的 模 数 互 不 相 同 。 那么称 g 为模 p 的原根,当且仅当 { g 0 , g 1 , … , g p − 2 } 对 p 的模数互不相同。
这样可能会更好理解些?至于为什么是上面那个绿绿的,我也不知道啊 ……
其实定义很简单,对于一些比较特殊的n n 次单位根,它的0 , 1 , … , n − 1 0 , 1 , … , n − 1 次幂恰好生成了所有的n n 次单位根,我们把这样的单位根称为本原单位根。
在复数域 C C 上,ω n = exp 2 π i n = cos 2 π n + i sin 2 π n ω n = exp 2 π i n = cos 2 π n + i sin 2 π n 是一个本原单位根。
(以上不是重点,下面这句话才是重点)
在有限域(即模质数)上,本原单位根和数论中的原根有关!
这也是我补充本原单位根的原因了……因为后面我们还会提到它。
并且我们可以证明的是,模质数的原根总是存在的 (不知道为啥你就问度娘吧)。
并且原根的性质和单位根非常相似哦~
N T T N T T 里原根的性质
换句话说,在模 p p 意义下,g g 可以被看做一个 p − 1 p − 1 次本原单位根。
反正感性理解就好,结论不就是用来记的嘛……
设 n ∣ p − 1 n ∣ p − 1 ,我们不妨令 ω n = g p − 1 n ω n = g p − 1 n (这个其实很像单位根的求法),然后你就可以得到一大串东西了:(好像也不多 )
ω n n = ( g p − 1 n ) n = g p − 1 ≡ 1 ( mod p ) ω n n = ( g p − 1 n ) n = g p − 1 ≡ 1 ( mod p )
并且因为单位根的定义,你会知道 ω 0 n , ω 0 n , … , ω n − 1 n ω n 0 , ω n 0 , … , ω n n − 1 互不相同……这不就完了嘛 Q A Q Q A Q
因此此时的 ω n ω n 在模 p p 意义下具有 n n 次本原单位根的性质 ,所以我们利用这个东西来定义 D F T D F T 即可,这被称为 N T T N T T 。
原根的求法
求模质数原根的办法: 将 p − 1 p − 1 分解质因数,即 p − 1 = p a 1 1 p a 2 2 … p a n n p − 1 = p 1 a 1 p 2 a 2 … p n a n 为标准的p − 1 p − 1 分解式,那么若 g p − 1 p i ≠ 1 ( mod p ) g p − 1 p i ≠ 1 ( mod p ) ,则g g 是模 p p 意义下的一个原根。
求模合数原根的办法: 其实你只需要把上面的 p − 1 p − 1 换成 φ ( p ) φ ( p ) 即可。
证明:(假设法)
先对于第一种模质数的情况:
假设存在一个 t < φ ( p ) = p − 1 t < φ ( p ) = p − 1 使得 g t ≡ 1 ( mod p ) g t ≡ 1 ( mod p ) 且 ∀ i ∈ [ 1 , n ] , g p − 1 p i ≠ 1 ( mod p ) ∀ i ∈ [ 1 , n ] , g p − 1 p i ≠ 1 ( mod p ) 。由裴蜀定理可知:一定存在一组 ( a , b ) ( a , b ) 使得:a ∗ t = ( p − 1 ) ∗ b + g c d ( t , p − 1 ) a ∗ t = ( p − 1 ) ∗ b + g c d ( t , p − 1 ) 成立。由欧拉定理/费马小定理:g p − 1 ≡ 1 ( mod p ) g p − 1 ≡ 1 ( mod p ) 。所以 g a ∗ t ≡ g ( p − 1 ) ∗ b + g c d ( t , p − 1 ) ≡ 1 ( mod p ) g a ∗ t ≡ g ( p − 1 ) ∗ b + g c d ( t , p − 1 ) ≡ 1 ( mod p ) 。因为 t < p − 1 t < p − 1 ,所以又有 g c d ( t , p − 1 ) < p − 1 g c d ( t , p − 1 ) < p − 1 。又因为 g c d ( t , p − 1 ) ∣ p − 1 g c d ( t , p − 1 ) ∣ p − 1 ,于是 ∃ i ∈ [ 1 , n ] , g c d ( t , p − 1 ) ∣ g p − 1 p i ∃ i ∈ [ 1 , n ] , g c d ( t , p − 1 ) ∣ g p − 1 p i 。那么就可以得到 g p − 1 p i ≡ g g c d ( t , p − 1 ) ≡ 1 ( mod p ) g p − 1 p i ≡ g g c d ( t , p − 1 ) ≡ 1 ( mod p ) 。
因此假设不成立。所以原命题成立。
Q . E . D . Q . E . D .
那么上述结论自然可以推广到模合数而把 p − 1 p − 1 改为 φ ( p ) φ ( p ) 了。
以上的证明是针对于这一句话的:δ m ( a ) = φ ( m ) δ m ( a ) = φ ( m ) ,也就是上述求法成立的条件。我们从另一个性质出发,即 g 0 , g 1 , … , g p − 1 g 0 , g 1 , … , g p − 1 模p p 的余数各不相同,那么可以得知这些余数的集合为:{ 1 , 2 , … p − 1 } { 1 , 2 , … p − 1 } (不分顺序)。我们知道取模是满足乘法性质的。也就是说,如果∃ i ∈ [ 1 , n ] ∃ i ∈ [ 1 , n ] 使得g p − 1 p i = 1 g p − 1 p i = 1 ,我们知道唯有 g 0 mod p = 1 mod p = 1 g 0 mod p = 1 mod p = 1 ,那么很明显 p − 1 p i p − 1 p i 是不会为 0 0 的,那么一旦出现第一个数字的重复,也就是 1 1 ,那么很明显不满足原根的上述性质。
那么我就很好奇,为什么重复的1 1 一定会出现在某个g p − 1 p i g p − 1 p i 当中呢?上面其实也已经证明,如果存在某个 t < φ ( p ) t < φ ( p ) 使得 g t ≡ 1 ( mod p ) g t ≡ 1 ( mod p ) 成立,那么必定会存在一个 g p − 1 p i mod p = 1 g p − 1 p i mod p = 1 。并且我认为这样可能更好理解。
反正证明这种东西感性理解一下就好了……
然后原根你也会求了,所以用原根 g g 去代替 F F T F F T 中的 ω n ω n (即单位根),就可以完成 F F T F F T 了。
NTT 模数
首先根据 F F T F F T 我们知道 g p − 1 n g p − 1 n 中的 n n 是 2 2 的幂次,所以在做 F F T F F T 的时候我们最好构造形如 p = a ∗ 2 k + 1 p = a ∗ 2 k + 1 的质数,这样就可以满足上面的需求了。这样的质数最好取大一点(精度),所以这样的质数可以取:
p = a ∗ 2 k + 1 p = a ∗ 2 k + 1
a a
k k
g g
3 3
1 1
1 1
2 2
5 5
1 1
2 2
2 2
17 17
1 1
4 4
3 3
97 97
3 3
5 5
5 5
193 193
3 3
6 6
5 5
257 257
1 1
8 8
3 3
7681 7681
15 15
9 9
17 17
12289 12289
3 3
12 12
11 11
40961 40961
5 5
13 13
3 3
65537 65537
1 1
16 16
3 3
786433 786433
3 3
18 18
10 10
5767169 5767169
11 11
19 19
3 3
7340033 7340033
7 7
20 20
3 3
23068673 23068673
11 11
21 21
3 3
104857601 104857601
25 25
22 22
3 3
167772161 167772161
5 5
25 25
3 3
469762049 469762049
7 7
26 26
3 3
998244353 998244353
119 119
23 23
3 3
1004535809 1004535809
479 479
21 21
3 3
1998585857 1998585857
953 953
21 21
3 3
2013265921 2013265921
15 15
27 27
31 31
2281701377 2281701377
17 17
27 27
3 3
3221225473 3221225473
3 3
30 30
5 5
75161927681 75161927681
35 35
31 31
3 3
77309411329 77309411329
9 9
33 33
7 7
206158430209 206158430209
3 3
36 36
22 22
2061584302081 2061584302081
15 15
37 37
7 7
2748779069441 2748779069441
5 5
39 39
3 3
6597069766657 6597069766657
3 3
41 41
5 5
39582418599937 39582418599937
9 9
42 42
5 5
79164837199873 79164837199873
9 9
43 43
5 5
263882790666241 263882790666241
15 15
44 44
7 7
1231453023109121 1231453023109121
35 35
45 45
3 3
1337006139375617 1337006139375617
19 19
46 46
3 3
3799912185593857 3799912185593857
27 27
47 47
5 5
4222124650659841 4222124650659841
15 15
48 48
10 10
7881299347898369 7881299347898369
7 7
50 50
6 6
31525197391593473 31525197391593473
7 7
52 52
3 3
180143985094819841 180143985094819841
5 5
55 55
6 6
1945555039024054273 1945555039024054273
27 27
56 56
5 5
4179340454199820289 4179340454199820289
29 29
57 57
3 3
理论到实现
第一件事,我们先来回忆一下 F F T F F T ,而 F F T F F T 分为两个步骤:D F T D F T 与 I D F T I D F T 。
然后我们知道单位根的定义是对于复数域 C C ,有 z n = 1 z n = 1 的复数就是一个 n n 次单位根,n n 次单位根有 n n 个,为:e 2 π k n i , k = { 0 , 1 , 2 … , n − 1 } e 2 π k n i , k = { 0 , 1 , 2 … , n − 1 } (i i 是虚数单位)。(你用三角表示也是可以的,但是不好理解)
然后这里就用 g p − 1 n g p − 1 n 代替 上面那个就可以,因为你会发现他们具有相同的性质。
唯一具有卷积性质的变换就是 D F T D F T ,按照上面的式子,D F T D F T 的变换式就是:
X k = N − 1 ∑ n = 0 x n e − 2 π N n k i k = 0 , 1 , 2 … , N − 1 X k = ∑ n = 0 N − 1 x n e − 2 π N n k i k = 0 , 1 , 2 … , N − 1
D F T D F T 反演(也就是逆变换)公式有:
x n = 1 N N − 1 ∑ k = 0 X k e 2 π N n k i n = 0 , 1 , 2 … N − 1 x n = 1 N ∑ k = 0 N − 1 X k e 2 π N n k i n = 0 , 1 , 2 … N − 1
那么由于 g p − 1 n g p − 1 n 和 e 2 π k n i , k = { 0 , 1 , 2 … , n − 1 } e 2 π k n i , k = { 0 , 1 , 2 … , n − 1 } 等价,所以:
e − 2 π N n k i ≡ g p − 1 N ( mod p ) e − 2 π N n k i ≡ g p − 1 N ( mod p )
接着我们就得到了快速数论变换 的公式:
X k = N − 1 ∑ n = 0 x n g n k ( mod p ) k = 0 , 1 , 2 … , N − 1 X k = ∑ n = 0 N − 1 x n g n k ( mod p ) k = 0 , 1 , 2 … , N − 1
反演:
x n = 1 N N − 1 ∑ k = 0 X k g − n k ( mod p ) n = 0 , 1 , 2 … N − 1 x n = 1 N ∑ k = 0 N − 1 X k g − n k ( mod p ) n = 0 , 1 , 2 … N − 1
这样我们就成功把复数域转到了整数域,然后剩下的东西在mod p mod p 意义下考虑即可。上面也已经讲过素数的构造方法。
考虑单位根的性质(也就是我们为什么使用单位根):
n n 个单位根互不相同。那么g 0 , g 1 , … , g p − 2 g 0 , g 1 , … , g p − 2 在模 p p 意义下也不相同,成立。
z n = 1 z n = 1 (z ∈ C z ∈ C )。根据费马小定理:g p − 1 ≡ 1 ( mod p ) g p − 1 ≡ 1 ( mod p ) ,成立。
单位根对称分布。其实这对于原根也成立。
证明:
∵ ( g p − 1 2 ) 2 ≡ g p − 1 ≡ 1 ( mod p ) ∵ ( g p − 1 2 ) 2 ≡ g p − 1 ≡ 1 ( mod p )
又∵ g i ≠ g j ( i ≠ j 且 0 ≤ i , j ≤ p − 1 ) ∵ g i ≠ g j ( i ≠ j 且 0 ≤ i , j ≤ p − 1 )
∴ g p − 1 2 =≠ 1 , 则 g p − 1 2 = − 1 ∴ g p − 1 2 =≠ 1 , 则 g p − 1 2 = − 1
Q . E . D . Q . E . D .
好了,我们真的可以使用 N T T N T T 来完成 F F T F F T 所能完成和所不能完成 的事了。
可以尝试过一下模板题:P 3803 P 3803 【模板】多项式乘法(F F T F F T )
虽然题目没有要求模数,但是仍然可以取一个比较大的模数从而完成。
这里取 p = 998244353 , g = 3 p = 998244353 , g = 3 。
我们结合代码来理解 N T T N T T :
然后真的比F F T F F T 快好多呀……模板题评测记录:
递 归 F F T 递归 F F T :点我
蝴 蝶 变 换 F F T 蝴蝶变换 F F T :点我
N T T N T T :点我
Upd20200918:文章当中还有一些错误,近期本人正在矫正。
__EOF__
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 如何编写易于单元测试的代码
· 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】