快速数论变换(fast number-theoretic transform,FNTT)

‎2024‎年‎1‎月‎15‎日 12:55:59 记

大数乘法与多项式乘法(Polynomial Multiplication)

两个大数如$12345$和$67890$，用多项式表达：

在进行多项式乘法的时候，得到$h(x)$的系数便可以变相表示为两个大整数之积

转为列表

逆序一下，符合高位在左的书写习惯

进位

快速傅里叶变换(Fast Fourier Transform,FFT)

回到开始的多项式，估算原积不超过$10$位，不足$2^k$位数则用$0$补齐（方便后续递归分治）

正常的乘法，需要对$f(x)$和$g(x)$的每一项都两两相乘，复杂度$O(n^2)$很高

为了描述这些多项式方程，除了采用系数表示法($n$个系数)，还可以使用描点表示法($n$个点)，例如方程$y=1x+0$可以用点$(0,0),(1,1)$唯一确定，方程$y=1x^2+2x+3$可以用点$(0,3),(1,6),(2,11)$唯一确定，即具有$n$个系数的$(n-1)$次多项式方程可以由$n$个不同的点唯一确定

那么对于$7$次多项式方程

任取$7$个

得到$7$个函数值

为了求这$7$个函数值，可以用矩阵方法表示(用列向量主要是因为横向量太长了Typora写不下)

如若使用编程实现该功能，复杂度$O(n^2)$依然没有优势

重新回到$f(x)$，wiki上对奇偶函数的代数结构有这样的描述：

偶函数的任何线性组合皆为偶函数，且偶函数会形成一个实数上的向量空间。相似地，奇函数的任何线性组合皆为奇函数，且奇函数亦会形成一个实数上的向量空间，实际上，”所有“实值函数之向量空间为偶函数和奇函数之子空间的直和。换句话说，每个定义域关于原点对称的函数都可以被唯一地写成一个偶函数和一个奇函数的相加：

类似地，将原式$f(x)$奇偶性分解

记

则有

如此一来计算$f(x_i)$的值只需计算$q_e(x_i),q_o(x_i)$的值即可，观察到$q(x)$与$f(x)$形似，也可以使用同样的逻辑进一步拆分分治，使用递归的逻辑连接起来。

为了更好的利用奇偶性质，$x_i$的值不任取，如在本例$n=8$中选取

这样一来

综上写出伪代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
def fft(a):
	#a = [a_0,a_1,a_2,...,a_{n-1}]
	#return [f(w^0),f(w^1),f(w^2),...,f(w^{n-1})]
	if (n==1) return [a_0]
	pe = [a_0,a_2,...,a_{n-2}]
	po = [a_1,a_3,...,a_{n-1}]
	q_o = fft(p_e)
    q_e = fft(p_e)
    w = cos(2pi/n) + i sin(2pi/n)
    for j in [0,1,2,...,n/2-1]:
    	f(w^j) = q_e[j] + w^j q_o[j]
    	f(w^{j+n/2}) = q_e[j] - w^j q_o[j]
    return [f(w^0),f(w^1),f(w^2),...,f(w^{n-1})]

以上，实现了从系数表示到点列表示的过程，即主要是为了实现这个计算

p.s:突然发现这好像就是我当时参加的TACA的原题啊喂

不妨把这个大矩阵记作$W$

想要实现从点列表示转成系数表示，就需要求得这个$W$矩阵的逆矩阵$W^{-1}$

离散傅里叶变换矩阵(Discrete Fourier Transform matrix,DFTmtx)

离散傅立叶变换矩阵是将离散傅立叶变换以矩阵乘法来表达的一种表示式

DFT矩阵定义

或

注意：此处$\omega$在指数上是否存在负号区别于之前的$\omega$，同时此处矩阵$W$有正规化因数$\frac{1}{\sqrt{N}}$

其逆矩阵是

下验证$WT=I$

当且仅当$i=j$时有

当$i \ne j$时有

至此$WT=I$验证完毕

快速傅里叶逆变换(Inverse Fast Fourier Transform,IFFT)

根据上面的式子，我们得到了

那么回到最先的式子

可以有

即

由于

这意味着先前得出的FFT逻辑完全可以把所有的$\omega$换成$\tau$便可完成IFFT的过程，再一次实现函数复用

现在剩余一个问题，$\omega$是一个复数，在$n$不大的情况下计算误差不大，但是当$n$很大，计算误差和三角函数计算时间都是问题，于是要引入有限域

有限域(finite field)

一点儿没学明白，查资料也写不出来多少$QAQ$

大致运算：

如在$\mathrm{mod} \space 5 $下，

根据此类性质，可以写出C函数如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
typedef long long int lli;
lli ZmodAdd(lli a, lli b){
    a += b;
    return (a<MOD_P) ? a : (a-MOD_P);
}

lli ZmodSub(lli a, lli b){
    a -= b;
    return (a>=0) ? a : (a+MOD_P);
}

lli ZmodMul(lli a, lli b){
    return (a * b) % MOD_P;
}

lli ZmodInv(lli a) {
    //根据费马小定理，乘法逆元唯一存在
    //利用扩展欧几里得算法求乘法逆元
    lli x,y;
    gcdext(a,MOD_P,&x,&y);
    //x mod p 即为所求逆元，x的范围是正负p/2
    return (x>=0) ? x : (x+MOD_P);
}

lli ZmodDiv(lli a, lli b) {
    return (a * ZmodInv(b)) % MOD_P;
}

lli gcdext(lli a,lli b,lli *x,lli *y){  //扩展欧几里得算法，网上抄的
    if(b==0){
        *x=1;
        *y=0;
        return a;
    }
    lli ret=gcdext(b,a%b,x,y);
    lli t=*x;
    *x=*y;
    *y=t-a/b*(*y);
    return ret;
}

快速数论变换(fast number-theoretic transform,FNTT)

先前FFT中选取$\omega = e^{\frac{2i\pi}{n}}$的性质主要有三：

  1. 所有计算除了系数$a_i$只需要计算$\omega ^{k}$即可，便于计算
  2. $\omega ^{k+\frac{n}{2}}=-\omega ^{k},\omega^{n}=1$
  3. $\omega ^{k},k\in Z_{n}$各不相同

而在一个有限域内，也同样有数可以满足这两条特性，如$p=84906529,\omega =213016,n=16$，$\omega^k \pmod{p} $的各值如下

这样的$\omega$也同样符合先前的性质，因此不妨把FFT中的所有运算都转为模$p$数域下的运算，$\omega$换成整数的$\omega$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int FNTT(int *pa,int w,List a,int *f){ //FNTT快速数论变换（基于FFT快速傅里叶变换）
    if (a.n == 1){
        f[0] = pa[a.begin];
        return 0;
    }
    int nOver2 = a.n>>1;
    int stepBy2 = a.step<<1;
    List pe = createList(a.begin, stepBy2, nOver2);
    List po = createList(a.begin+a.step, stepBy2, nOver2);
    int qe[nOver2];
    int qo[nOver2];
    FNTT(pa, ZmodMul(w,w), pe, qe);
    FNTT(pa, ZmodMul(w,w), po, qo);
    int t = 1;
    int qo_t;
    for (int j=0; j<nOver2; j++){
        qe_t = ZmodMul(qo[j],t);
        f[j] = ZmodAdd(qe[j] , qo_t);
        f[j+nOver2] = ZmodSub(qe[j] , qo_t);
        t = ZmodMul(w,t);
    }
    return 0;
}

那么如何选择$p,\omega$？

对于质数$p = q2^m+1$，原根$g$满足$g^{q2^k} \equiv 1\pmod{p}$，将$g_n \equiv g^q \pmod{p}$看作$\omega$的等价，则其满足相似的性质

p.s:数论知识好多，不会证明;w;

综上，写成完整C语言代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#include <stdio.h>
#define MOD_P 84906529
typedef long long int lli;

typedef struct {
    int begin;
    int step;
    int n;
}List;
List createList(int begin,int step,int n) {
    List newList;
    newList.begin = begin;
    newList.step = step;
    newList.n = n;
    return newList;
}

//Zmod(p)数域计算（c语言好像没有重载运算符，所以表达式看起来会比较麻烦~）
lli ZmodAdd(lli a, lli b);
lli ZmodSub(lli a, lli b);
lli ZmodMul(lli a, lli b);
lli ZmodInv(lli a);
lli ZmodDiv(lli a, lli b);
lli gcdext(lli a,lli b,lli *x,lli *y);  //扩展欧几里得算法

int FNTT(int *pa,int w,List a,int *f){ //FNTT快速数论变换（基于FFT快速傅里叶变换）
    if (a.n == 1){
        f[0] = pa[a.begin];
        return 0;
    }
    int nOver2 = a.n>>1;
    int stepBy2 = a.step<<1;
    List pe = createList(a.begin, stepBy2, nOver2);
    List po = createList(a.begin+a.step, stepBy2, nOver2);
    int qe[nOver2];
    int qo[nOver2];
    FNTT(pa, ZmodMul(w,w), pe, qe);
    FNTT(pa, ZmodMul(w,w), po, qo);
    int t = 1;
    int qo_t;
    for (int j=0; j<nOver2; j++){
        qe_t = ZmodMul(qo[j],t);
        f[j] = ZmodAdd(qe[j] , qo_t);
        f[j+nOver2] = ZmodSub(qe[j] , qo_t);
        t = ZmodMul(w,t);
    }
    return 0;
}


int main() {
    int w = 213016;
    int t = ZmodInv(w);
    int f[16]     = {5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0};
    int g[16]     = {1,2,3,4,5,0,0,0,0,0,0,0,0,0,0,0};
    int ntt_f[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    int ntt_g[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    int h[16]     = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    int r[16]     = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    List start = createList(0,1,16);
    FNTT(f,w,start,ntt_f);
    FNTT(g,w,start,ntt_g);
    for (int i=0; i<16; i++){
        h[i] = ZmodMul(ntt_f[i],ntt_g[i]);
    }
    FNTT(h,t,start,r);
    for (int i=0; i<16; i++){
        r[i] = ZmodDiv(r[i],16);
    }
    for (int i=0; i<16; i++){
        printf("%d ",r[i]);
    }
    return 0;
}

lli ZmodAdd(lli a, lli b){
    a += b;
    return (a<MOD_P) ? a : (a-MOD_P);
}

lli ZmodSub(lli a, lli b){
    a -= b;
    return (a>=0) ? a : (a+MOD_P);
}

lli ZmodMul(lli a, lli b){
    return (a * b) % MOD_P;
}

lli ZmodInv(lli a) {
    //根据费马小定理，乘法逆元唯一存在
    //利用扩展欧几里得算法求乘法逆元
    lli x,y;
    gcdext(a,MOD_P,&x,&y);
    //x mod p 即为所求逆元，x的范围是正负p/2
    return (x>=0) ? x : (x+MOD_P);
}

lli ZmodDiv(lli a, lli b) {
    return (a * ZmodInv(b)) % MOD_P;
}

lli gcdext(lli a,lli b,lli *x,lli *y){  //扩展欧几里得算法，网上抄的
    if(b==0){
        *x=1;
        *y=0;
        return a;
    }
    lli ret=gcdext(b,a%b,x,y);
    lli t=*x;
    *x=*y;
    *y=t-a/b*(*y);
    return ret;
}

参考学习

超硬核FFT快速傅里叶变换讲解，高效进行高精度乘法运算！ - Bilibili

最讨厌做算法题的时候出现浮点数了，有什么方法可以避免？斐波那契数列和FFT均可适用！ - Bilibili

ChatGPT

快速数论变换 - OI wiki

Finite field - Wikipedia

离散傅里叶变换矩阵 - Wikipedia

奇函数与偶函数 - Wikipedia

__END__