快速傅里叶变换与快速数论变换（一）

最新推荐文章于 2024-09-29 15:06:28 发布

Ivanzn

最新推荐文章于 2024-09-29 15:06:28 发布

阅读量1.8k

点赞数 1

本文链接：https://blog.csdn.net/Ivanzn/article/details/89059495

版权

转载请注明

4.3.1递归分治法(自顶向下，常数比较大)

1.前言

我计划将分成两章讲解“快速傅里叶变换与快速数论变换”。这一篇博客是我在研究和学习了3~4天的FFT和FNT的基础之上诞生的。

本片博客主要讲解FFT和FNT的最基本的知识。

FFT和FNT是用来解决数个多项式相乘的问题的。

2.前提知识

2.1秦九韶算法

2.1.1介绍

秦九韶算法（在西方叫做“霍纳法则”）能在 O(n) 的时间内计算次多项式的值

该多项式可以表示为： $a_0+a_1 x+a_2x^2+...+a_{n-1}x^{n-1}$

2.1.2算法的具体过程

$a_0+a_1 x+a_2x^2+...+a_{n-1}x^{n-1}$

$\Rightarrow a_0+(a_1+a_2x+...+a_{n-1}x^{n-2})x$

$\Rightarrow a_0+(a_1+(a_2+..+a_{n-1}x^{n-3})x)x$

$\Rightarrow a_0+(a_1+(a_2+(..a_{n-1}x..)x)x)x$

2.1.3 CODE

double f (int n,double a[],double x){//n代表a[]的大小，数组的传递需要指定数组大小
    int i;                           //a[]储存每个x项的系数值，在函数外初始化
    double p=a[n];//初始化p；
    for(i=n;i>0;i--)
        p=a[i-1]+x*p;
    return p;
}

2.2复数

2.2.1虚数

设 x^2=-1

则 $x=\pm \sqrt {-1}$ ,此方程在实数集上没有解。于是，我们再设置一个复数集。

在复数集上，我们设 i^2=-1 ,于是此方程也有解。其中，我们把叫做虚数

2.2.2复数及其属性

我们把实数+虚数结合体叫做复数，比如： 1+2i

其中，是实部，是虚部。

设复数 Z=a+bi ，则的模长为 $\left | Z \right |=\sqrt{a^2+b^2}$ ,辐角为 $\theta = arctan \frac{b}{a}$

2.2.3复数的运算法则

设 Z_1=a+bi , Z_2=c+di

则 $Z_1 \cdot Z_2= (ac-bd)+(ad+bc)i$

则模长为 $\sqrt{(a^2+b^2)(c^2+d^2)}$ ，辐角为 $\theta=arctan \frac{ad+bc}{ac-bd}=\theta_1+\theta_2$

因此，有结论：两个复数相乘：模长相乘，辐角相加

2.2.4一种特殊的复数

一种特殊的复数： $e^{i\varphi }=cos \varphi +i \cdot sin \varphi$ , 它的模长是 $\sqrt{cos \varphi ^2+sin \varphi^2 }=1$ , 辐角为 $\theta= arctan \frac{b}{a}=\varphi$ 。

这个复数可以用来旋转一个复数.

比如，设 Z=a+bi , 再设 $Z_1=Z \cdot e^{i\varphi}$

由上面的结论，我们得到： $\left | Z_1 \right |=\left | Z \right |=\sqrt{a^2+b^2}$ , $\theta=arctan\frac{b}{a}+\varphi$

即复数顺时针旋转了 $\theta$

2.2.5欧拉公式

令上面的 $\varphi=\pi$ ,得到 $e^{i\pi}+1=0$

这就是著名的欧拉公式。

2.3多项式的两种表示方式

设两个次多项式分别为 $A(x)= \sum_{i=0}^{n-1}a_ix^i$ ， $B(x)= \sum_{i=0}^{n-1}b_ix^i$

该多项式在数学上有两种表示方式，分别为：系数向量法、点值法

注意，在这里，我们用任意一个变量取代，也就是把上面这些式子当作形式幂级数（形式幂级数的概念在母函数中提到过）

2.3.1系数向量法

所谓的系数向量法就是利用多项式中未知数前的系数构成的一组系数向量（我们将它看成是列向量）

即 $(a_0,a_1,...a_{n-1})$ 和 $(b_0,b_1,...b_{n-1})$

若采用系数向量法来表示两个多项式相乘，

则有 $c_k=\sum_{i=0}^{k}a_ib_{k-i}$ $(k \in [0,2n-2])$

时间复杂度为 O(n^2)

2.3.2点值法

所谓的点值法就是用个点值对表示该多项式，每一个点值对由和 a(x) 构成

即 $\{(x_0,a(x_0)),(x_1,a(x_1)),...(x_{n-1},a(x_{n-1})) \}$

点值法的本质就是取个样本，从而构成一个阶方阵。然后就可以唯一确定一组系数向量，从而唯一确定一个多项式

2.3.3两者之间的联系

系数法 $\rightarrow$ 点值法 “求值”运算

点值法 $\rightarrow$ 系数法 “插值”运算

2.4单位根

2.4.1基本概念

设 $\omega$ 为一个复数，令 $\omega ^n=1$ ，我们称其为次单位复数根（就是模长为的复数的次等于的复数集）

2.4.2性质

有以下性质：

性质一：复数集的基数是，也就是有个这样的复数

性质二：复数集中的复数均匀分布在复平面上的单位圆上，将单位圆均匀的分成了个部分

性质三：复数集中的各个复数的模长为

性质四：复数集中的元素的形式为 $\omega =e^{\frac{2\pi ki}{n}}$ ,其中 $k \in[0,n-1]$ 。我们把 $\omega_n=e^{\frac{2\pi i}{n}}$ 叫做主次单位复数根

性质五：由性质四，得到复数集中各元素分别为： $\omega_{n}^{0},\omega_{n}^{1},..\omega_{n}^{n-1}$

性质五：复数集是一个群，满足以下恒等式

一、 $\omega_{n}^{j} \cdot \omega_{n}^{k}=\omega_{n}^{j+k(mod \ n)}$

二、 $\omega_{dn}^{dk} = \omega_{n}^{k}$ （消去引理）

三、（由恒等式二得到本推论） $\omega_{n}^{n/2}=\omega_{2}^{1}=-1$

四、 $\omega_{n}^{0} = \omega_{n}^{n}=1$

五、 $\sum_{i=0}^{n-1}(\omega_{n}^{i})^p=0$ ,其中，是一个常数

六、 $\omega_{n}^{n/2+k}=\omega_{n}^{n/2}\cdot \omega_{n}^{k}=\omega_2\cdot\omega_{n}^{k}=-\omega_{n}^{k}$ （折半引理）

备注：折半引理是FFT得以成立的充要条件

3.FFT的计算过程概要

polynomial-multiplication

（图片摘自图片原文地址）

STEP 1：将系数向量A，B转换为点值法 ------------------ O(nlogn)

STEP 2：将A,B两个点值相乘 --------------------------------- O(n)

STEP 3：将乘好以后的点值转换为新的系数向量C------ O(nlogn)

因此，可以得到 FFT 的时间复杂度为 O(nlogn)

我们将STEP 1 叫做 DFT（即离散傅里叶变换），将STEP 2 叫做 IDFT（即逆离散傅里叶变换）

STEP 1 和 STEP 2互为逆运算

接下来，我将分成了两个部分（DFT和IDFT）来讲解FFT是怎么算的

不过，在此之前，约定下文中的是的幂次（即 n=2^m ，不到 2^m ，就用补齐）

4.DFT（离散傅里叶变换）

是将系数向量转化为点值，是求值运算

设多项式为 $A(x)= \sum_{i=0}^{n-1}a_ix^i$

设多项式的系数向量 $(a_0,a_1,...a_{n-1})$ 为 $\underset{a}{\rightarrow}$ （已知）

设用点值法表示为 $\{(x_0,a(x_0)),(x_1,a(x_1)),...(x_{n-1},a(x_{n-1})) \}$

4.1单位复数根的点值法表示多项式

接下来，用次单位复数根中的个复数和其对应的多项式的值，分别代替点值法中的 x_i 和 A(x_i)

从而得到一组新的点值式： $\{(\omega_{n}^{0},A(\omega_{n}^{0}) ),(\omega_{n}^{1},A(\omega_{n}^{1})),...(\omega_{n}^{n-1},A(\omega_{n}^{n-1})) \}$ ，记作 $\underset{y}{\rightarrow}$ （未知）

则通过 DFT_n 运算： $\underset{y}{\rightarrow}$ = DFT_n $\underset{a}{\rightarrow}$ , 得到 $\underset{a}{\rightarrow}$ 的离散傅里叶变换值 $\underset{y}{\rightarrow}$

4.2 DFT（）

设 $A(\omega _{n}^{k})=\sum_{i=0}^{n-1}a_{i}(\omega_{n}^{k})^i$

再设 $A_1(\omega _{n}^{k})=\sum_{i=0}^{n-1}a_{2i}(\omega_{n}^{k})^i$ , $A_2(\omega _{n}^{k})=\sum_{i=0}^{n-1}a_{2i+1}(\omega_{n}^{k})^i$

那么就有 $A(\omega_{n}^{k})=A_1(\omega_{n}^{2k})+A_2(\omega_{n}^{2k})$

由 2.4.2 的消去引理得到： $A(\omega_{n}^{k})=A_1(\omega_{n/2}^{k})+A_2(\omega_{n/2}^{k})\cdot \omega_{n}^{k}$ ------------------------ (1)

由 2.4.2 的折半引理得到： $A(\omega_{n}^{k+n/2})=A_1(\omega_{n/2}^{k})-A_2(\omega_{n/2}^{k})\cdot \omega_{n}^{k}$ ----------------- (2)

观察 (1) 和 (2) ，我们得到一个结论：每次计算一个 $A(\omega_{n}^{k})$ 之后，就能直接计算得到 $A(\omega_{n}^{k+n/2})$

根据这个结论，就能使得原问题的规模缩小到原来的一半，因此， DFT 的时间复杂度就是 O(nlogn)

讲的通俗一点就是：

$A(\omega_{n}^{0})=A(\omega_{n}^{n/2})$

$A(\omega_{n}^{1})=A(\omega_{n}^{n/2+1})$

$A(\omega_{n}^{n/2-1})=A(\omega_{n}^{n-1})$

4.3 DFT的两种求法

4.3.1递归分治法(自顶向下，常数比较大)

以下是伪代码:

DFT ( $\underset{a}{\rightarrow}$ ):

$\underset{a}{\rightarrow}$ .length

$if(n==1) \ return$ $\underset{a}{\rightarrow}$

$\omega=1$

$\omega_n=e^{\frac{2\pi i}{n}}$

$y_{[0]}=DFT(a_0,a_2,a_4...)$

$y_{[1]}=DFT(a_1,a_3,a_5...)$

$for \ k=0 \ to \ n/2-1$

$y_k=y_{k}^{[0]}+y_{k}^{[1]}\cdot \omega$

$y_k=y_{k}^{[0]}-y_{k}^{[1]}\cdot \omega$

$\omega=\omega \cdot \omega_n$

$return \ y$

4.3.2迭代法（自底向上，常数比较小）

以 n=16 为例子，观察下图中系数向量在分治法中的变化：

bit-reverse

（图片摘自图片原文地址）

STEP 1 : 是全序排列

STEP 2 : 末位为0的在左树，末位为1的在右树

STEP 3 : 倒数第二位为0的在左树，倒数第二位为1的在右树

STEP 4 : 倒数第三位为0的在左树，倒数第三位为1的在右树

考虑用迭代的方法计算 DFT ，那么只要知道 STEP log_2 n 的系数向量的顺序就可以自底向上迭代得到 $\underset{y}{\rightarrow}$

再观察 STEP 4 中个系数的顺序，我们发现：

0	8	4	12	2	10	6	14	1	9	5	13	3	11	7	15
0000	1000	0100	1100	0010	1010	0110	1110	0001	1001	0101	1101	0011	1011	0111	1111
0000	0001	0010	0011	0100	0101	0110	0111	1000	1001	1010	1011	1100	1101	1110	1111
0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15

将 STEP 4 中个系数的二进制翻转后的数是逐一递增的。（我们把这个叫做逆序置换）

我们假设 rev(i) 是将的二进制翻转之后的数，那么我们就可以得到以下的伪代码：

reverse ( $\underset{a}{\rightarrow}$ )

$for \ i =0 \ to$ $\underset{a}{\rightarrow}$ .length

a(i) = rev(i)

return $\underset{a}{\rightarrow}$

DFT ( $\underset{a}{\rightarrow}$ )

$\underset{A}{\rightarrow}$ = reverse ( $\underset{a}{\rightarrow}$ )

$\underset{A}{\rightarrow}$ .length

$for \ s=1 \ to \ log_2n:$

m=2^s

$\omega_{n}=e^{\frac{2\pi i}{m}}$

$for \ j=0 \ to \ n-1\ by \ m:$

$\omega=1$

$for \ k=0 \ to \ m/2-1:$

$t=\omega \cdot A[j+k+m/2]$

u=A[j+k]

A[j+k]=u+t

A[j+k+m/2]=u-t

$\omega=\omega \cdot \omega_n$

return $\underset{A}{\rightarrow}$

reverse （逆序置换）的具体代码：

//采用DP
int rev[N];
void reverse(int n)
{
    memset(rev,0,sizeof(rev));
    for(int i=0;i<n;++i)
        rev[i]=(rev[i/2]>>1)|((i&1)<<(n-1));//(i&1)表示当前的项的奇偶，rev[i]的右n-1位与rev[i/2]的右n-1位相同
}

5. IDFT（逆离散傅里叶变换）

是将点值转化为系数向量，是插值运算

设多项式为 $A(x)= \sum_{i=0}^{n-1}a_ix^i$

设多项式的系数向量 $(a_0,a_1,...a_{n-1})$ 为 $\underset{a}{\rightarrow}$ （未知）

通过 DFT 得到一组点值 $\{ (\omega_{n}^{0},A(\omega_{n}^{0})),(\omega_{n}^{1},A(\omega_{n}^{1}))...(\omega_{n}^{n-1},A(\omega_{n}^{n-1})) \}$ ，记作 $\underset{y}{\rightarrow}$ (已知)

则有 $\underset{a}{\rightarrow}$ = IDFT $\underset{y}{\rightarrow}$ ----------- (*)

5.1 IDFT与逆矩阵

计算 (*) ,本质上就是解下列方程组：

写成矩阵的形式就是：

（以上图片截自图片原文地址）

这样， IDFT 就相当于把 DFT 过程中的 $\omega_{n}^{i}$ 换成 $\omega_{n}^{-i}$ ，然后做一次 DFT ，之后结果除以就可以了。

5.2 伪代码

把上面两种 DFT 伪代码中的 $\underset{a}{\rightarrow}$ 换成 $\underset{y}{\rightarrow}$ ， $\underset{y}{\rightarrow}$ 换成 $\underset{a}{\rightarrow}$ 即可。最后，每一项元素除以。

6.卷积定理

设 a,b 是两个次（是的幂次，不到就用补齐）多项式，那么有：

$a \bigotimes b=IDFT_{2n}(DFT_{2n}(a) \cdot DFT_{2n}(b))$ （ $\bigotimes$ 表示卷积）

注意：这里不是，而是！因为 a,b 这两个次多项式相乘的话，就会多产生 n-1 项

比如：

a=1+2x+3x^2 ------------------------------------------------- (1)

b=1+2x+3x^2 -------------------------------------------------- (2)

那么 $a \bigotimes b=1 + 4 x + 10 x^2 + 12 x^3 + 9 x^4$ ------- (3)

观察得到， (3) 比 (1) 或 (2) 多了项，所以要把原来的扩展为

7.快速数论变换

我们发现，使用次单位复数根可以简化多项式的计算，但是复数的计算难免会产生精度的误差。在对大整数的多项式的计算过程中就会失去优势，因此，我们将采用一种新的方式来代替了原来的次单位复数根。

7.1原根

设一个素数为，根据费马大定理有： $a^{p-1}\equiv 1 \ mod \ p$ $\Leftrightarrow \omega_{n}^{n} = 1$

若某个数是的原根，则有 $\forall \ 1\leq x < y \leq p-1 , \ g^x \neq g^y \ mod \ p$

根据原根的特点，我们设一个素数 $p=a \cdot 2^b+1$ ,

设 $g_n=g^{\frac{p-1}{n}}$ $\Leftrightarrow \omega_n=e^{\frac{2 \pi i}{n}}$ （对应主次单位复数根）

有 $g_{n}^{n} = g^{(\frac{p-1}{n})\cdot n}=1 \ mod \ p$ $\Leftrightarrow \omega_{n}^{n}=1$ （对应次单位复数根性质一）

有 $g_{n}^{0},g_{n}^{0},..g_{n}^{n-1}$ 互不相同 $\Leftrightarrow \omega_{n}^{0},\omega_{n}^{1},...\omega_{n}^{n-1}$ 互不相同（对应次单位复数根性质二）

有 $g_{n}^{n/2}=\sqrt{g_{n}^{n}}$ , 于是 $g_{n}^{n/2}=\pm 1$ ，

又因为 $g_{n}^{0}=1$ 和 $g_{n}^{0}\neq g_{n}^{n/2} \ mod \ p$ ,所以 $g_{n}^{n/2}=-1$ $\Leftrightarrow \omega_{n}^{n/2}=-1$ (对应次单位复数根性质三)

有了上面三条一一对应的性质，我们就可以使用这样的对应的 g_n 代替次单位复数根。

7.2 常用模数和对应原根的表

r⋅2k+1	r	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
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	19
7881299347898369	7	50	6
31525197391593473	7	52	3
180143985094819841	5	55	6
1945555039024054273	27	56	5
4179340454199820289	29	57	3