快速傅里叶变换学习笔记

(如果读到任何有问题的地方，请先阅读文章末尾)

（这里做一个约定，下文中的$\mod n$表示对$n$取模的剩余）

　　快速傅里叶变换(Fast Fourier Transform，简称FFT)常用于加速卷积运算(例如多项式乘法、大整数乘法)。

虚数和复数

当学习平方根的时候会不会出现这样一个疑问等于多少？

虚数就是这么一个数。那么形如的数，也都可以表示成的形式。

那么如果一个实数和一个虚数相加呢？很遗憾它们并不能合并到一起，所以只能写成形如的形式，这就是复数。其中被称为实部，被称为虚部。

复数的直角坐标系

考虑将直角坐标系的横轴与所有实数一一对应，纵轴与所有虚数一一对应，因此我们就可以将一个复数当成一个向量表示在这个平面上，这个平面又被称为复平面

因此定义复数的模长

复数的加减法

复数的加减法很简单，就是将实部与实部相加减，虚部与虚部相加减(和向量的加减法相同)。

例如设，则。

复数的乘法

复数的乘法法则和多项式的乘法法则相同，但是可以稍作化简:

设，则

共轭复数

两个复数实部相等，虚部互为相反数，则它们互为共轭复数(Conjugate complex number)，复数的共轭复数记作。关于共轭复数有一些神奇的性质。

(1)

(2)

复数的除法

复数的除法等于将两个多项式相除，不过通常希望将分母有理化，就可以将分子和分母同时乘分母的共轭复数即可。

复数乘法的几何意义

设的模长为，和轴正半轴的夹角为，的模长为，和轴正半轴的夹角为。

则有，，

因此得出结论，复数的乘法就是将两个复数的模长相乘，与轴正半轴的夹角相加。

单位复数根

次单位复数根是指满足的所有，我们将次单位复数根记作。不同于在实数域内次单位根不会超过2个（$n$是奇数的时候有$1$个，$n$是偶数的时候有$2$个），在复数域内次单位根有个。

现在我们思考如何求出这些单位根。

首先可以确定。现在来思考它的次幂后的它要和实数轴重合。我们设它原本和轴正半轴的夹角为，则有，因此可以解得。

可以得到。

　　再根据一个定义，简化一下我们的结论，得到。特别地，我们发现所有次单位根所在的直线恰好均分整个复平面(可以自己画一下图)。所以这些次单位根构成了一个循环群，对乘法运算具有封闭性。

对于我们记为次主单位根，其他的都可以表示成它的次幂。

多项式乘法

对于两个多项式和，如果把它们相乘，那么实质上是对它们的系数向量做了一次卷积:。

如果采用传统的系数表示法(用系数去表示一个多项式)，时间复杂度大概是，很难再进行优化。

但是考虑如果用点值表示法(一个次多项式就用对点，如果每一个都不相同，则这个多项式是唯一确定的)，那么乘法只需要将对应的点值乘起来，这么做的话时间复杂度是的，但是对于朴素的转换成点值表达法和进行插值(还原系数表示的多项式)是。

但是对于某些特殊的取值再加上一些黑科技优化可以使得这两步变成的。不过首先给出一些引理。

消去引理 对于给定整数，那么有

证明。

特别地，当且为正整数时，有。

折半引理 当且为正整数时有

证明，然后就可以证出了。

求和引理 对于满足的，则有

证明使用等比数列求和公式:

因为，所以保证了分母不为0。

离散傅里叶变换

对于一个多项式

将我们取好个次单位根带入求值，设系数向量，结果向量。其中。结果向量就是系数向量的离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），也记作。

快速傅里叶变换

由于两个次多项式相乘，最终会变成次多项式，所以我们需要多取一些点，并使得总点数是2的倍数(如果项数不是2的倍数就把它拓展成2的倍数，多出来的系数补零)。

现在考虑计算，如果我们把每个单位根都代入，用秦九韶算法进行求值，总时间复杂度也是的，现在来考虑用基于折半引理的黑科技分治来优化。

首先将多项式分成两个有项的多项式

然后不难得到。

这样就可以将一个有项的多项式分割成两个只有项的多项式。那这两个多项式如何选点呢？

这个担心是多余，因为我们要计算这两个多项式在

的取值时，其实是等价于计算这两个多项式计算在下面

这些点的取值，根据消去引理，这里实际上只有个点，它们分别是

然后这个问题就可以递归处理。

　　既然是递归，那应当有边界，最后当系数向量中只有一个元素的时候就返回(因为)。

　　现在来考虑下一层已经计算出下一层的结果向量了，现在来考虑求值。

对于被扔进左边递归的结果有

而被扔进右边递归的结果有

　　通过上面的步骤，我们能够在$n\log n$的时间内把多项式对应的点值求出来。

但是每次按下标奇偶对系数进行分组，而且又是递归，特别耗时间，所以考虑将递归改成非递归。另外为了更好地处理并且减小常数，考虑预处理出最终系数的序列。

这个有一个现成的算法名为雷德算法，可以快速完成位置的调整。感兴趣的话可以出门左转百度(反正网上讲雷德算法的似乎比讲FFT的还多)。这里只给出代码，不再继续讨论。

 1 void Rader(Complex<double> *f, int len) {
 2     for(int i = 1, j = len >> 1, k; i < len - 1; i++) {
 3         if(i < j)
 4             swap(f[i], f[j]);
 5         
 6         k = len >> 1;
 7         while(j >= k) {
 8             j -= k;
 9             k >>= 1;
10         }
11         
12         if(j < k)
13             j += k;
14     }
15 }

　　有了上面的预处理，就可以用非递归，从解答树的叶节点开始计算。

离散傅里叶逆变换

现在通过计算对应的点值，再将其乘起来，得到了新的一堆点值，现在要做的事是将它还原成系数表示法。

回顾我们得到结果向量的过程，本质可以看作进行了一次矩阵乘法:

现在我们知道结果向量，希望求系数向量，那么可以用结果向量乘左边这个范德蒙德矩阵的逆矩阵就好了。可是求一个矩阵的逆矩阵并不是什么容易的事，直接套公式的话，应该时间复杂度和高斯消元是一个级别的吧。但是幸运的是这里有这么一个结论，就是。

证明的思路是证这个矩阵和范德蒙德矩阵的乘积是单位矩阵。

则根据矩阵乘法法则有(请自行把这个当成我们需要求证的这个矩阵，不能直接把它当成的逆矩阵)：

当时，显然结果等于1。当时，原式等价于，再根据求和引理得到它的值为0。所以这个矩阵是的逆矩阵。

然后我们再用结果向量乘上，得到

再和讲DFT的时候，有关的定义的式子展开作比较

然后发现区别并不大，所以可以用FFT去解决只不过，加入一个参数sign就好了。最后再将所有的除以，于是大功告捷。

　　如果不太清楚实现可以看代码，是uoj的第34题，多项式乘法。

  1 /**
  2  * UOJ
  3  * Problem#34
  4  * Accepted
  5  * Time:3402ms
  6  * Memory:17632k
  7  */
  8 #include <bits/stdc++.h>
  9 typedef bool boolean;
 10 
 11 //Complex Temp Begin
 12 namespace yyf {
 13 
 14     template<typename T>
 15     class Complex {
 16         public:
 17             T r;
 18             T v;
 19             Complex(const T r = 0, const T v = 0):r(r), v(v) {        }
 20             
 21             Complex operator * (Complex b) {
 22                 return Complex(r * b.r - v * b.v, r * b.v + v * b.r);
 23             }
 24             
 25             Complex operator + (Complex b) {
 26                 return Complex(r + b.r, v + b.v);
 27             }
 28             
 29             Complex operator - (Complex b) {
 30                 return Complex(r - b.r, v - b.v);
 31             }
 32             
 33             Complex operator / (T a) {
 34                 return Complex(r / a, v / a);
 35             }
 36             
 37             Complex operator / (Complex b) {
 38                 Complex c = b.conj();
 39                 return Complex(r, v) * c / ((b * c).r);
 40             }
 41             
 42             //Conjugate complex
 43             inline Complex conj() {
 44                 return Complex(r, -v);
 45             }
 46             
 47             boolean operator == (Complex b) {
 48                 return r == b.r && v == b.v;
 49             }
 50     };
 51     
 52     template<typename T>
 53     boolean operator != (Complex<T> a, Complex<T> b) {
 54         return !(a == b);
 55     }
 56     
 57     template<typename T>
 58     Complex<T> operator += (Complex<T> &a, Complex<T> b) {
 59         return (a = a + b);
 60     }
 61     
 62     template<typename T>
 63     Complex<T> operator -= (Complex<T> &a, Complex<T> b) {
 64         return (a = a - b);
 65     }
 66     
 67     template<typename T>
 68     Complex<T> operator *= (Complex<T> &a, Complex<T> b) {
 69         return (a = a * b);
 70     }
 71     
 72     template<typename T>
 73     Complex<T> operator /= (Complex<T> &a, T x) {
 74         return (a = a / x);
 75     }
 76     
 77     template<typename T>
 78     Complex<T> operator /= (Complex<T> &a, Complex<T> b) {
 79         return (a = a / b);
 80     }
 81     
 82 }
 83 // Complex Temp Ends
 84 
 85 using namespace yyf;
 86 using namespace std;
 87 
 88 const double PI = acos(-1);
 89 const double eps = 1e-4;
 90 
 91 int n, m;
 92 Complex<double> a[262150], b[262150];
 93 
 94 inline void init() {
 95     scanf("%d%d", &n, &m);
 96     n++, m++;
 97     for(int i = 0, x; i < n; i++)
 98         scanf("%d", &x), a[i].r = x;
 99     for(int i = 0, x; i < m; i++)
100         scanf("%d", &x), b[i].r = x;
101 }
102 
103 void fft(Complex<double> *f, int len, int sign) {
104     if(len == 1)    return;
105     Complex<double> ql[len >> 1], qr[len >> 1];
106     for(int i = 0; i < (len >> 1); i++)
107         ql[i] = f[i << 1], qr[i] = f[(i << 1) | 1];
108     fft(ql, len >> 1, sign);
109     fft(qr, len >> 1, sign);
110     Complex<double> wn(cos(2 * PI / len), sin(sign * 2 * PI / len)), w(1, 0), t;
111     for(int i = 0; i < (len >> 1); i++, w *= wn) {
112         t = w * qr[i];
113         f[i] = ql[i] + t, f[i + (len >> 1)] = ql[i] - t;
114     }
115 }
116 
117 inline void solve() {
118     m += n, n = 1;
119     while(n < m)    n <<= 1;
120     fft(a, n, 1);
121     fft(b, n, 1);
122     for(int i = 0; i < n; i++)
123         a[i] = a[i] * b[i];
124     fft(a, n, -1);
125     for(int i = 0; i < m - 1; i++)
126         printf("%d ", (int)(a[i].r / n + 0.5));
127 }
128 
129 int main() {
130     init();
131     solve();
132     return 0;
133 }

FFT(recursive)

  1 /**
  2  * UOJ
  3  * Problem#34
  4  * Accepted
  5  * Time:709ms
  6  * Memory:9432k
  7  */
  8 #include <bits/stdc++.h>
  9 typedef bool boolean;
 10 
 11 //Complex Temp Begin
 12 namespace yyf {
 13 
 14     template<typename T>
 15     class Complex {
 16         public:
 17             T r;
 18             T v;
 19             Complex(const T r = 0, const T v = 0):r(r), v(v) {        }
 20             
 21             Complex operator * (Complex b) {
 22                 return Complex(r * b.r - v * b.v, r * b.v + v * b.r);
 23             }
 24             
 25             Complex operator + (Complex b) {
 26                 return Complex(r + b.r, v + b.v);
 27             }
 28             
 29             Complex operator - (Complex b) {
 30                 return Complex(r - b.r, v - b.v);
 31             }
 32             
 33             Complex operator / (T a) {
 34                 return Complex(r / a, v / a);
 35             }
 36             
 37             Complex operator / (Complex b) {
 38                 Complex c = b.conj();
 39                 return Complex(r, v) * c / ((b * c).r);
 40             }
 41             
 42             //Conjugate complex
 43             inline Complex conj() {
 44                 return Complex(r, -v);
 45             }
 46             
 47             boolean operator == (Complex b) {
 48                 return r == b.r && v == b.v;
 49             }
 50     };
 51     
 52     template<typename T>
 53     boolean operator != (Complex<T> a, Complex<T> b) {
 54         return !(a == b);
 55     }
 56     
 57     template<typename T>
 58     Complex<T> operator += (Complex<T> &a, Complex<T> b) {
 59         return (a = a + b);
 60     }
 61     
 62     template<typename T>
 63     Complex<T> operator -= (Complex<T> &a, Complex<T> b) {
 64         return (a = a - b);
 65     }
 66     
 67     template<typename T>
 68     Complex<T> operator *= (Complex<T> &a, Complex<T> b) {
 69         return (a = a * b);
 70     }
 71     
 72     template<typename T>
 73     Complex<T> operator /= (Complex<T> &a, T x) {
 74         return (a = a / x);
 75     }
 76     
 77     template<typename T>
 78     Complex<T> operator /= (Complex<T> &a, Complex<T> b) {
 79         return (a = a / b);
 80     }
 81     
 82 }
 83 // Complex Temp Ends
 84 
 85 using namespace yyf;
 86 using namespace std;
 87 
 88 const double PI = acos(-1);
 89 const double eps = 1e-4;
 90 
 91 int n, m;
 92 Complex<double> a[262150], b[262150];
 93 
 94 inline void init() {
 95     scanf("%d%d", &n, &m);
 96     n++, m++;
 97     for(int i = 0, x; i < n; i++)
 98         scanf("%d", &x), a[i].r = x;
 99     for(int i = 0, x; i < m; i++)
100         scanf("%d", &x), b[i].r = x;
101 }
102 
103 void Rader(Complex<double> *f, int len) {
104     for(int i = 1, j = len >> 1, k; i < len - 1; i++) {
105         if(i < j)
106             swap(f[i], f[j]);
107         
108         k = len >> 1;
109         while(j >= k) {
110             j -= k;
111             k >>= 1;
112         }
113         
114         if(j < k)
115             j += k;
116     }
117 }
118 
119 void fft(Complex<double> *f, int len, int sign) {
120     Rader(f, len);
121     for(int l = 2; l <= len; l <<= 1) {
122         Complex<double> wn(cos(2 * PI / l), sin(sign * 2 * PI / l)), u, v;
123         for(int s = 0; s < len; s += l) {
124             Complex<double> w(1, 0);
125             for(int i = s; i < s + (l >> 1); i++, w *= wn) {
126                 u = f[i];
127                 v = w * f[(i + (l >> 1))];
128                 f[i] = u + v, f[(i + (l >> 1))] = u - v;
129             }
130         }
131     }
132     
133     if(sign == -1)
134         for(int i = 0; i < len; i++)
135             f[i] = f[i] / len;
136 }
137 
138 inline void solve() {
139     m += n, n = 1;
140     while(n < m)    n <<= 1;
141     fft(a, n, 1);
142     fft(b, n, 1);
143     for(int i = 0; i < n; i++)
144         a[i] = a[i] * b[i];
145     fft(a, n, -1);
146     for(int i = 0; i < m - 1; i++)
147         printf("%d ", (int)(a[i].r + 0.5));
148 }
149 
150 int main() {
151     init();
152     solve();
153     return 0;
154 }

FFT(non-recursive)

在后面乱七八糟的东西

　　这次还是先写在word上。所以为了方便复制公式，就用LaTeX在线公式生成作成图片，然而Word 2013可以直接复制图片，但是，竟然不能批量复制。。。所以可能还有地方粘漏了。如果读到什么地方突然缺东西，多半是因为公式粘漏了。如果发现了请一定要在在评论区指出。

　　如果遗漏的地方影响阅读了，可以试试下载[快速傅里叶变换.rar]，只不过这个就没有代码了。

　　由于本人蒟蒻一发，所以如果此文有错误(我想多半是有)，请一定要在评论区指出！欢迎交流或者提问！

参考资料

1.《算法导论》

2.《数学一本通》

3. 快速傅里叶变换FFT 易懂版

posted @ 2017-07-31 19:38 阿波罗2003 阅读(1048) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部