快速傅里叶变换(FFT),常用于解答多项式乘法相关内容。
吐槽时间:
见了网上太多太多的资料都是物理的QAQ
不给OIers活路这样真的好么……
都是什么频率啊振幅啊什么乱七八糟的东西TAT
欺负我什么都不懂是吧QAQ
还是自己写一个比较好QAQ
背景故事:
在我们平时计算多项式乘法的时候,我们把第一个多项式的每一项都和第二个多项式的每一项相乘,复杂度为O(n ^ 2),此时我们所使用的表示法就是系数表示法。
现在我们可以有一种比较强大的方式:
点值表示法:
众所周知,假设f(x)的最高次数为
n−1
,即次数界为
n
,那么我们只要知道了n个不相同的x及f(x)值,就能确定出f(x)的多项式。
有一种算法叫做秦九韶算法,它能得出的结论是:
对于一个n次多项式,至多做n次乘法和n次加法。
所以我们知道了n个点的x值之后,我们在O(n^2)的时间内就能计算出所有的y值,然而如果通过快速傅里叶变换的方法,可以在O(nlogn)的时间内求出所有的y值。
引入新定义:
求值:通过多项式的系数表示法求其点值表示法。
插值:通过多项式的点值表示法求其系数表示法。
显然上面两个定义是互逆的关系,FFT就是用来解决求值的过程的。
引入新定义:
n次单位复根:在复平面内,n次单位复根能把整个平面分成n块。它的严格定义是:满足
由复数幂的定义,可知:
(ω1n)k=ωkn
它有很多性质:
1.相消引理:
这是显然的,你可以考虑这个复平面等分成6个平面取第2个根和把这个复平面等分成3个平面取第一个根没什么区别。
2.折半引理:
证明:左边 = ω2k+nn=ω2kn∗ωnn=ω2kn=右边
3.几何级数:
对于在复数域上的x ≠1 ,有下列等式成立:
嗯其实就是等比数列求和。
4.求和引理:
对于任意:n >= 1 且 n不能整除非零整数k ,有:
证明:
由几何级数的计算过程可知:
证毕。
接下来我们分析一下如何得到最终的多项式吧。
1.求A的n个单位根的点值,求B的n个单位根的点值
2.点值相乘,得到C的点值。
3.计算C的多项式。
我们先考虑第一步。
我们希望计算多项式:
在n个单位根处的值。
设多项式A以系数形式给出,系数向量 a⃗ =(a0,a1,...,an−1) ,定义其结果
则向量 y⃗ =y0,y1,...,yn−1 就是系数向量 a⃗ 的 离散傅里叶变换(DFT),记作:y = DFTn(a) .
也就是说,点值的结果就是DFT,那么我们只需要快速计算DFT的值即可。
如果按照正常的算法,时间复杂度显然是
O(n2)
的,所以有快速傅里叶变换(FFT),它采取的是分治的思想。
我们考虑原来的多项式A(x),定义两个新的次数界为n / 2的多项式:
则有
我们如果只是这样的话是没办法减少计算量的,因为这只是相当于展开。
我们考虑一个显然的事实:
对于能够计算 A(x) 的那两个多项式,实际上它们也是可以用FFT求值的, 它们的次数界减少了一半。
我们仔细回忆下单位根的性质,发现:
化简一下第二个式子,发现:
嗯……
这个玩意叫做 蝴蝶变换。
减少了一半的计算量。
我们来计算下复杂度:
T(n) = 2 * T(n / 2) + n ;
T(n / 2) = 2 * T(n / 4) + n / 2;
…
显然,复杂度为 nlog2n 。
我们终于讲完了如何求出n个单位根的点值,现在来说第二步:
把点值相乘。
显而易见的一点是,我们如果有两个多项式A(x),B(x),它们相乘得到多项式C(x),则:
那么我们把n个点值相乘的复杂度自然是O(n)的。
只需要记录c[i] = a[i] * b[i]即可。
接下来我们要做的就是插值了。[忘了定义的请自觉往回翻]
因为我不知道拉格朗日插值公式是什么鬼(划掉),所以我决定还是不讲了。
//别打我别打我……我讲…………
嗯我们考虑把DFT写成矩阵的形式:
终于挤下来了。
我们发现这个X的矩阵叫做范德蒙德矩阵。
发现我们现在要求的就是 a⃗ ,即 DFT−1n(y) 。
那么怎么求系数矩阵呢?
则:若有矩阵AB = C,已知C,A则有:B = A−1C
证明: A−1∗AB=B=A−1∗C
证毕。
显然的一点是:
范德蒙德矩阵的逆是存在的,我们想求 a⃗ 也不是很难,但是,矩阵乘法在这里的复杂度是 O(n2) 的。
对于中间那个矩阵,我们设其为 Vn ,则有:
V−1n的(j,k)的元素为ω−kjn/n
证明:
……不想写了。
我们给定逆矩阵 V−1n ,可以推导出 DFT−1n(y) :
然而对比之前的式子……?
长得一模一样……
用a数组直接存储a数组 * b数组的点值,对当前的a数组做FFT,对其结果除以n就能得到答案。
但是我们现在的东西还停留在递归阶段,我们知道,可以对这个数组进行元素上的调整,这样我们可以先计算那些子节点,省去了递归的过程,这个过程叫做:
雷(二)德(进)算(制)法
大概意思就是把每一个数的二进制头尾翻转过来:
即100 - > 001这样子……
我们可以知道的一件事情就是:0 - > 7这个序列,在递归到叶子节点之后,是:0,4,2,6,1,5,3,7。我们发现把二进制转过来之后按照顺序排序之后,再反转回来就能得到原来的序列了。
嗯……发现自己还是介绍一下雷德算法比较好,它有一个O(n)求每个数的二进制的反转的方法,但是请记住,这个O(n)并不如平时我们写的O(nlogn)要快,但是极其优美。
我们考虑一个数x的二进制反转和x / 2的二进制反转,发现如果x是偶数的话,那么[x / 2]的实际二进制和x的实际二进制只差了在二进制中移一次位。用实际二进制数组表达,那么设二进制数组为A[],则A[x] = A[x >> 1] << 1;
在反转的二进制中,设其数组为Rev[]:
发现,一个数如果是奇数,那么也仅仅是改变了右移之后的最左端的一位。因为:A[x + 1] = A[x] | 1,那么我们就能同样改变反转数组:
左移的位数大家仔细考虑一下发现,我们只需要使用 < n 的那些数字,也就意味着,在0 - > 3区间内:
Rev[0] = 00,Rev[1] = 10,Rev[2] = 01,Rev[3] = 11。
发现如果我们左移Log次好像也是……
可以的……吗?
显然是不能的,因为我们希望的是:
在将其二进制反转之后,用O(n)的方法对其调整计算的顺序。
我们如果二进制每次左移log位,看样子是可行的,但是我们只用到了0 - > n的,而左移log位就会使得除0以外的二进制反转>=n,因为用快速排序会使得整个程序再次陷入 O(nlog2n) ,所以我们每次显然不能左移log位。
我们发现: 左移(log - 1)位之后,O(n)扫一遍整个序列就能排序。
void Rader()
{
//m 是原来两个多项式中较长的那个的长度的2倍。
for(n = 1;n <= m;n <<= 1,LG ++);
for(int i = 1;i < n;i ++)
Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (LG - 1));
}
void FFT(Virt *F,int ty = 1)
{
Rep_0(i,n)if(i < Rev[i])swap(F[i],F[Rev[i]]);//每次FFT的过程中根据二进制的反转调整运算顺序。
...
}
这样的话我们合并相邻的两个,就像线段树一样,这样问题就成功解决了。
请用手写复数类并重载运算符。
请用手写复数类并重载运算符。
请用手写复数类并重载运算符。
题目:
BZOJ2179:给出两个n位10进制整数x和y,你需要计算x*y。
n <= 60000
//据说FFT不开二倍数组长度会跑得更快
//注意:以上内容不保证其正确性与真实性。
//注意:以上内容不保证其正确性与真实性。
//注意:以上内容不保证其正确性与真实性。
我们可以将两个数相乘写成这种形式:
然后直接输出就可以了。
bzoj2179,FFT模版,据说我这样不加注释会运行的更快。
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define Rep(i,n) for(int i = 1;i <= n;i ++)
#define RepG(i,x) for(int i = head[x];~ i;i = edge[i].next)
#define RD(i,x,n) for(int i = x;i <= n;i ++)
#define Rep_0(i,n) for(int i = 0;i < n;i ++)
#define Rep0n(i,n) for(int i = 0;i <= n;i ++)
#define PI M_PI
using namespace std;
struct Virt
{
double r,i;
Virt (double r = 0.0,double i = 0.0)
{
this -> r = r;
this -> i = i;
}
Virt operator+(const Virt &x)
{
return Virt(r + x.r,i + x.i);
}
Virt operator-(const Virt &x)
{
return Virt(r - x.r,i - x.i);
}
Virt operator*(const Virt &x)
{
return Virt(r * x.r - i * x.i,x.r * i + r * x.i);
}
};
const int N = 131072;
int n,m,LG,c[N],Rev[N];
Virt a[N],b[N];
char s[N >> 1];
void Rader()
{
for(n = 1;n <= m;n <<= 1,LG ++);
for(int i = 1;i < n;i ++)
Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (LG - 1));
}
void FFT(Virt *F,int ty = 1)
{
Rep_0(i,n)if(i < Rev[i])swap(F[i],F[Rev[i]]);
for(int i = 1;i < n;i <<= 1)
{
Virt wn(cos(PI / i),ty * sin(PI / i));
for(int j = 0;j < n;j += (i << 1))
{
Virt w(1,0);
for(int k = 0;k < i;k ++,w = w * wn)
{
Virt u = F[j + k],v = F[j + k + i] * w;
F[j + k] = u + v,F[j + k + i] = u - v;
}
}
}
if(ty == -1)Rep_0(i,n)F[i].r = F[i].r / n;
}
int main ()
{
scanf("%d",&n);
m = (n - 1) << 1;
scanf("%s",s);
Rep_0(i,n)a[i] = s[n - i - 1] - '0';
scanf("%s",s);
Rep_0(i,n)b[i] = s[n - i - 1] - '0';
Rader();
FFT(a),FFT(b);
Rep0n(i,n)a[i] = a[i] * b[i];
FFT(a,-1);
Rep0n(i,m)c[i] = (int)(a[i].r + 0.5);
Rep0n(i,m)
{
if(c[i] >= 10)
{
c[i + 1] += c[i] / 10,c[i] %= 10;
if(i == m)m ++;
}
}
for(int i = m;i >= 0;i --)printf("%d",c[i]);
return 0;
}
992ms应该也不算特别慢了吧……
再来看一道模版题:
快速傅立叶之二
请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n ,并且有 n < = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。
Input
第一行一个整数N,接下来N行,第i+2..i+N-1行,每行两个数,依次表示a[i],b[i] (0 < = i < N)。
Output
输出N行,每行一个整数,第i行输出C[i-1]。
原式 = > ?
让我来冷静一下……
FFT解决的卷积不是:
话说这两个到底有什么相同点啊!
于是我重新冷静了一下:
我们一开始解决的 Fk 在第k位存储的是能满足把a和b的第i位的元素看作 ai∗xi 这种形式,实际上根本没有所谓的 xi ,但是这给了我们一个提示……
从0 - > k: Fk 存储的是a和b乘起来之后x的次数等于k的项的系数。
这次的题目k - > n: Ck 存储的是a和b乘起来之后x的次数等于n + k的系数。
那就直接输出n - 1 - > 2 * (n - 1)位的结果即可。
其实这个题也给我提了个醒:
不要关心卷积的过程是什么,你只需要知道:
这个东西是一个卷积。
所以你只要求一下对应的次数的那些卷积就可以了。
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define Rep(i,n) for(int i = 1;i <= n;i ++)
#define RepG(i,x) for(int i = head[x];~ i;i = edge[i].next)
#define RD(i,x,n) for(int i = x;i <= n;i ++)
#define Rep_0(i,n) for(int i = 0;i < n;i ++)
#define Rep0n(i,n) for(int i = 0;i <= n;i ++)
#define PI M_PI
using namespace std;
struct Virt
{
double r,i;
Virt (double r = 0.0,double i = 0.0)
{
this -> r = r;
this -> i = i;
}
Virt operator+(const Virt &x)
{
return Virt(r + x.r,i + x.i);
}
Virt operator-(const Virt &x)
{
return Virt(r - x.r,i - x.i);
}
Virt operator*(const Virt &x)
{
return Virt(r * x.r - i * x.i,x.r * i + r * x.i);
}
};
const int N = 131072 * 2;
int n,m,LG,c[N],Rev[N];
Virt a[N],b[N];
char s[N >> 1];
void Rader()
{
for(n = 1;n <= m;n <<= 1,LG ++);
for(int i = 1;i < n;i ++)
Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (LG - 1));
}
void FFT(Virt *F,int ty = 1)
{
Rep_0(i,n)if(i < Rev[i])swap(F[i],F[Rev[i]]);
for(int i = 1;i < n;i <<= 1)
{
Virt wn(cos(PI / i),ty * sin(PI / i));
for(int j = 0;j < n;j += (i << 1))
{
Virt w(1,0);
for(int k = 0;k < i;k ++,w = w * wn)
{
Virt u = F[j + k],v = F[j + k + i] * w;
F[j + k] = u + v,F[j + k + i] = u - v;
}
}
}
if(ty == -1)Rep_0(i,n)F[i].r = F[i].r / n;
}
int main ()
{
scanf("%d",&n);
m = (n - 1) << 1;
Rep_0(i,n)
scanf("%lf%lf",&a[i].r,&b[n - i - 1].r);
Rader();
FFT(a),FFT(b);
Rep0n(i,n)a[i] = a[i] * b[i];
FFT(a,-1);
for(int i = m / 2;i <= m;i ++)printf("%d\n",c[i] = (int)(a[i].r + 0.5));
return 0;
}
参考资料:
ACdreamers的博客
http://blog.csdn.net/acdreamers/article/details/39005227
zky的博客
http://blog.csdn.net/iamzky/article/details/22712347
算法导论
//请扫描以下二维码关注本人微信。