快速傅里叶变换(FFT)学习笔记(其二)(NTT)

再探快速傅里叶变换(FFT)学习笔记(其二)(NTT)

目录

• 再探快速傅里叶变换(FFT)学习笔记(其二)(NTT)
  • 写在前面
    • 一些约定
  • 前置知识
    • 同余类和剩余系
    • 欧拉定理
    • 阶
    • 原根
      • 求原根
  • NTT
    • NTT的定义
    • 从单位根到原根
      • 常用NTT模数表
    • NTT的实现

写在前面

为了不使篇幅过长，预计将把学习笔记分为四部分:

1. DFT,IDFT,FFT的定义,实现与证明:快速傅里叶变换(FFT)学习笔记(其一)
2. NTT的实现与证明:快速傅里叶变换(FFT)学习笔记(其二)
3. 任意模数NTT与FFT的优化技巧
4. 多项式相关操作

一些约定

1. ([p(x)]=egin{cases}1,p(x)为真 \ 0,p(x)为假 end{cases})
2. 本文中序列的下标从0开始
3. 若(s)是一个序列,(|s|)表示(s)的长度
4. 若大写字母如(F(x))表示一个多项式,那么对应的小写字母如(f)表示多项式的每一项系数，即(F(x)=sum_{i=0}^{n-1} f_ix^i)
5. 由于本文主要讨论数论知识，本文中的未知数如未经说明，均为正整数。
6. 在处理同余式的时候，可能会出现(equiv)和(=)混用的情况,此时的(=)均代表(equiv)

前置知识

同余类和剩余系

定义1.1:若整数(a)和整数(b)除以正整数(m)的余数相等，则称(a,b)模(m)同余，记作(a equiv b (mod m))

定义1.2:把全体整数按照模(m)的余数分为若干个互不相交的集合,每一个集合称为模(m)的同余类(剩余类)(residue(congruence) class).余数为(a)的同余类记作(ar{a})

对于(forall a in [0,m-1]),(ar{a}={a+km }(k in mathbb{Z}))

同余类中的每个元素都可以拿来代表该同余类，称为该同余类的代表数

定义1.3:模(m)的同余类一共有(m)个,分别为(ar{0},ar{1},dots ar{m-1}).在它的每一个同余类中取一个元素作为代表，所有这些代表元素组成的集合称为(m)的完全剩余系(complete residue system)

定义1.4 ([1,m])中与(m)互质的数代表的同余类组成(m)的简化(既约)剩余系(reduced residue system)

例如,10的简化剩余系为({ar{1},ar{3},ar{7},ar{9}})

定理1.1 (m)的简化剩余系关于模(m)乘法封闭

证明:若(forall a,b in [1,m])与(m)互质,那么(ab)也不可能与(m)含有相同的质因子。根据辗转相除法我们知道(gcd(ab mod m,m)=gcd(ab,m)=1),也就是说(ab mod m)也属于(m)的简化剩余系

欧拉定理

我们知道(varphi(m))表示([1,m])中与(m)互质的b数的个数。那么易得(m)的简化剩余系的大小为(varphi(m))

定理2.1 若(gcd(a,n)=1),则(a^{varphi(n)} equiv 1 (mod n))

证明:设(n)的简化剩余系为({ar{a_1},ar{a_2},dots ar{a_{varphi(n)}} })

(forall a_i,a_j,a_i eq a_j (mod n))

又因为(gcd(a,n)=1), 故(a(a_i-a_j) eq 0 (mod n)),即(aa_i eq aa_j(mod n))

又因为简化剩余系关于模(n)乘法封闭,故(ar{aa_i})也在(n)的简化剩余系中。因此({ar{a_1},ar{a_2},dots ar{a_{varphi(n)}} })和({ar{aa_1},ar{aa_2},dots ar{aa_{varphi(n)}} })都可以表示(n)的简化剩余系(可以理解为对应顺序不同).所以有

(a^{varphi(n)}a_1a_2dots a_{varphi(n)} equiv a_1a_2dots a_{varphi(n)} (mod n))

所以(a^{varphi(n)} equiv 1 (mod n))

定理2.2,若(p)为质数且(gcd(a,p)=1),(a^{p-1} equiv 1 (mod p))

由定理2.1显然

定理2.1被称为欧拉定理,而定理2.2被称为费马小定理

阶

根据欧拉定理我们知道,存在正整数(d leq m-1),使得(a^d equiv 1 (mod m)).(显然(varphi(m) leq m-1),(m)为素数时取等).但是,最小的(d)是多少呢?

定义3.1:设(m>1),且(gcd(a,m)=1),那么使得(a^r equiv 1 (mod m))的最小正整数(r)称为(a)对模(m)的阶(指数)(order),记为(delta_m(a))

定理3.1: 对于正整数(d>1),(a^d equiv 1 (mod m))的充要条件是(delta_m(a)|d).

证明:

由于(m>1,gcd(a,m)=1),易得(forall j geq 0,a^j mod m e 0)((a)与(m)没有公共质因数,那么(a^j)也一定和(m)没有公共的质因数).

不妨设(a^j=q_jm+r_j(0<r_j<m))

这样,取(j=0,1,2 dots m-1),其中(m)个余数(r_0,r_1,dots r_{m-1})最多取(m-1)个值,根据抽屉原理,必有两个相等,记为(r_i,r_k).不妨设(0 leq k<i<m),那么(m(q_i-q_k)=a^i-a^k=a^k(a^{i-k}-1)). 所以(m|a^k(a^{i-k}-1))

又因为(gcd(a^k,m)=1),所以把互质的因数去掉也不影响整除性,所以(m|(a^{i-k}-1))，即(a^{i-k} equiv 1 (mod m)).取(d=i-k)就证明了该定理.

定理3.2:设(k)是非负整数,则有

[delta_m(a^k)=frac{delta_m(a)}{gcd(delta_m(a),k)} ]

设(t=gcd(delta_m(a),k)),则(k=q_1t,delta_m(a)=q_2t(q_1,q_2 in mathbb{N}))

令(delta_m(a^k)=s),则(a^{ks} equiv 1(mod m)).根据定理3.1,(q_2t|ks).

所以(q_2|ks),又因为(gcd(q_2,k)=1),所以(q_2|s).又根据原根的最小性知,(s)是第一个满足(q_2|s)的数，因此(q_2=s)。

也就是说(delta_m(a^k)=q_2=frac{delta_m(a)}{t}=frac{delta_m(a)}{gcd(delta_m(a),k)})

原根

定义4.1:设(m geq 1,gcd(a,m)=1),若(delta_m(a)=varphi(m)),则称(a)为(m)的原根

定理4.1:若(m)有原根(g),模(m)的简化剩余系可以表示为({ar{g^0},ar{g^1},dots ar{g^{varphi(m)-1}} })

证明:

引理4.1.1: (g^0,g^1 dots g^{varphi(m)-1})两两模(m)不同余

反证法,假设存在(0 leq i<j leq varphi(m)-1),使得(g^i equiv g^j (mod m)),那么(g^{j-i} equiv 1(mod m)),即(j-i< varphi(m)),这与原根的最小性矛盾

引理4.1.2:(g^0,g^1 dots g^{varphi(m)-1})均与(m)互质

由于(gcd(g,m)=1),所以(g^k mod m eq 0)

结合两个引理和简化剩余系的定义，我们证明了定理4.1

定理4.1的推论:若(p)是素数,(g^0,g^1 dots g^{p-2})恰好构成(1,2,dots p-1)的一个排列

证明:(varphi(p)=p-1),显然得证

定理4.2 (m)的原根有(varphi(varphi(m)))个

证明:根据定理4.2,(g^0,g^1 dots g^{varphi(m)-1})遍历(m)的简化剩余系,也就是说我们只需要考虑(g)的幂(g^k).

若(a^k)是(m)的原根,(delta_m(a^k)=varphi(m)=delta_m(a)).

根据定理3.2,(delta_m(a^k)=frac{delta_m(a)}{gcd(delta_m(a),k)})

那么(gcd(delta_m(a),k)=1),即(gcd(varphi(m),k)=1).满足条件的(k)的个数就是(varphi(varphi(m)))

求原根

定理4.3:对于任意正整数(m geq 2),(varphi(m))的所有的不同的素因数为(q_1,q_2 dots q_s),那么(g)是模(p)的原根的充要条件是:

[forall 1 leq j leq s,g^{frac{varphi(m)}{q_j}} eq 1 (mod m) ]

证明:

根据欧拉定理,(g^{varphi(m)}=1(mod m))

根据定理3.1,(delta_m(g)|varphi(m)).如果(varphi(m))的每个非自身的因数(d)都不满足条件(a^d=1(mod m))，那么根据原根的最小性,(g)就是原根。实际上我们不需要枚举每个因数,只需要枚举那些包含质因数最多的就可以了,因为如果有比(frac{varphi(m)}{q_j})小的数满足条件,同样根据定理3.1,(frac{varphi(m)}{q_j})也满足条件。也就是说,所有(frac{varphi(m)}{q_j})的因数集合的并集包含了除(varphi(m))外(varphi(m))的所有因数。

定理4.3的推论设(p)是奇素数,(p-1)的所有的不同的素因数为(q_1,q_2 dots q_s),那么(g)是模(p)的原根的充要条件是:

[forall 1 leq j leq s,g^{frac{p-1}{q_j}} eq 1 (mod p) ]

求原根没有什么太好的算法，只能暴力枚举,然后根据定理4.3判断，但好在原根一般都很小,很快能找到答案。

下标给出了一些常用大质数的原根，容易发现原根都很小

p	原根g
3	2
5	2
17	3
97	5
193	5
257	3
40961	3
65537	3
104857601	3
167772161	3
469762049	3
998244353	3
1004535809	3

NTT

NTT的定义

快速数论变换(Number Theory Transform,NTT)可以用于在模意义下计算两个整数序列的卷积，避免了FFT的精度问题。

从单位根到原根

我们发现,FFT用到了复数单位根的五大性质

1. (omega_n^0=omega_n^n=1)
2. (forall i eq j,0<i,j<n,omega_n^i eq omega_n^j)
3. (omega_{2n}^{2k}=omega_n^k)
4. (omega_n^{k+frac{n}{2}}=-omega_n^k)
5. (forall v in mathbb{N},frac{1}{n} sum_{k=0}^{n-1} omega_n^{v k}=[v mod n=0])

想到我们之前提到的原根,如果把单位根换成原根,结果会如何呢?

一般的NTT的模数只能是形如(p=a imes 2^b+1)的质数.设(p)的原根为(g),我们定义(g_n^k=(g^{frac{p-1}{n}})^{k} (mod p))

下面我们证明(g)满足单位根的5条性质,注意所有计算在模(p)下进行:

定理5.1:(g_n^0=g_n^n=1 (mod p))

证明:显然(g_n^0=1),(g_n^n=g^{p-1}=g^{varphi(p)}=1 (mod p)),最后一步用到了原根的定义

定理5.2 (forall i eq j(mod (p-1)),g_n^i eq g_n^j (mod p))

证明:根据定理4.1,(g^0,g^1 dots g^{p-2})遍历(p)的剩余系,证毕。一定要注意,(i,j)是在模(p-1)的意义下

定理5.3 (g_{2n}^{2k}=g_n^k)

证明:$$g_{2 n}^{2 k} equivleft(g^{frac{p-1}{2 n}} ight)^{2 k} equivleft(g^{frac{p-1}{n}} ight)^{k} equiv g_{n}^{k} quad(mod p)$$

定理5.4:(g_n^{k+frac{n}{2}}=-g_n^k (mod p))

由于(p)是质数,且(g_n^n=1 (mod p)),因此(g_n^{n/2}= ±1)

由定理(5.2)得(g_n^{n/2} eq g_n^1=1)

所以(g_n^{n/2}=-1 (mod p)),(g_n^{k+frac{n}{2}}=-g_n^k)

定理5.5 $$forall v in mathbb{N},frac{1}{n} sum_{k=0}^{n-1} omega_n^{v k}=[v mod n=0]$$

当(v mod n eq 0)时,根据等比数列求和公式

[egin{aligned} sum_{k=0}^{n-1}left(g_{n}^{v} ight)^{k} &=frac{left(g_{n}^{v} ight)^{n}-1}{g_{n}^{v}-1} \ &=frac{left(g_{n}^{n} ight)^{v}-1}{g_{n}^{v}-1} \ &=frac{1-1}{g_{n}^{v}-1} \ &=0 end{aligned}]

当(v mod n=0)时,(g_n^v=g^{frac{(p-1)v}{n}}),那么((p-1)|frac{(p-1)v}{n}),根据定理3.1,(g_n^v=1),故(g_n^{vk}=1), (sum_{i=0}^{n-1}left(g_{n}^{v} ight)^{i}=sum_{i=0}^{n-1} 1=n)

证毕。

至此，我们发现可以用(g_n)替代单位根. 注意到(frac{p-1}{n})必须是整数,考虑到FFT的过程是不断二分,那么(p-1)必须有足够的质因子2,这就是为什么一般的NTT模数都是(p=a imes 2^b+1)的形式

常用NTT模数表

(p=a imes 2^b+1)	(a)	(b)	原根(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

节选自该链接

NTT的实现

void NTT(ll *x,int n,int type){
    for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]); 
    for(int len=1;len<n;len*=2){
        int sz=len*2;
        ll gn1=fast_pow((type==1?G:invG),(mod-1)/sz);
        for(int l=0;l<n;l+=sz){
            int r=l+len-1;
            ll gnk=1;
            for(int i=l;i<=r;i++){
                ll tmp=x[i+len];
                x[i+len]=(x[i]-gnk*tmp%mod+mod)%mod;
                x[i]=(x[i]+gnk*tmp%mod)%mod;
                gnk=gnk*gn1%mod; 
            }
        }
    } 
    if(type==-1){
        ll invn=inv(n);
        for(int i=0;i<n;i++) x[i]=x[i]*invn%mod; 
    }
} 

可以与FFT代码对比