首页 > 其他分享 >FFT转载

FFT转载

时间：2024-04-21 17:22:06浏览次数：15  
标签：nk nn FFT kn a2 n2 转载 0n

快速傅里叶变换(FFT)详解

原文链接：快速傅里叶变换(FFT)详解 - 自为风月马前卒 - 博客园 (cnblogs.com)

目录

 

本文只讨论FFT在信息学奥赛中的应用

文中内容均为个人理解，如有错误请指出，不胜感激

回到顶部

前言

先解释几个比较容易混淆的缩写吧

DFT：离散傅里叶变换—>O(n2)O(n2)计算多项式乘法

FFT：快速傅里叶变换—>O(n∗log(n)O(n∗log⁡(n)计算多项式乘法

FNTT/NTT：快速傅里叶变换的优化版—>优化常数及误差

FWT：快速沃尔什变换—>利用类似FFT的东西解决一类卷积问题

MTT：毛爷爷的FFT—>非常nb/任意模数

FMT 快速莫比乌斯变化—>感谢stump提供

回到顶部

多项式

系数表示法

设A(x)A(x)表示一个n−1n−1次多项式

则A(x)=∑ni=0ai∗xiA(x)=∑i=0nai∗xi

例如：A(3)=2+3∗x+x2A(3)=2+3∗x+x2

利用这种方法计算多项式乘法复杂度为O(n2)O(n2)

（第一个多项式中每个系数都需要与第二个多项式的每个系数相乘）

点值表示法

将nn互不相同的xx带入多项式，会得到nn个不同的取值yy

则该多项式被这nn个点(x1,y1),(x2,y2),…,(xn,yn)(x1,y1),(x2,y2),…,(xn,yn)唯一确定

其中yi=∑n−1j=0aj∗xjiyi=∑j=0n−1aj∗xij

例如：上面的例子用点值表示法可以为(0,2),(1,6),(2,12)(0,2),(1,6),(2,12)

利用这种方法计算多项式乘法的时间复杂度仍然为O(n2)O(n2)

（选点O(n)O(n)，每次计算O(n)O(n)）

 

我们可以看到，两种方法的时间复杂度都为O(n2)O(n2)，我们考虑对其进行优化

对于第一种方法，由于每个点的系数都是固定的，想要优化比较困难

对于第二种方法，貌似也没有什么好的优化方法，不过当你看完下面的知识，或许就不这么想了

 

回到顶部

复数

在介绍复数之前，首先介绍一些可能会用到的东西

向量

同时具有大小和方向的量

在几何中通常用带有箭头的线段表示

圆的弧度制

等于半径长的圆弧所对的圆心角叫做1弧度的角，用符号rad表示，读作弧度。用弧度作单位来度量角的制度叫做弧度制

公式:

1∘=π180rad1∘=π180rad

180∘=πrad180∘=πrad

平行四边形定则

 

(好像画的不是很标准。。)

平行四边形定则：AB+AD=AC

 

复数

定义

设a,ba,b为实数，i2=−1i2=−1，形如a+bia+bi的数叫复数，其中ii被称为虚数单位，复数域是目前已知最大的域

在复平面中，xx代表实数，yy轴（除原点外的点）代表虚数，从原点(0,0)(0,0)到(a,b)(a,b)的向量表示复数a+bia+bi

模长：从原点(0,0)(0,0)到点(a,b)(a,b)的距离，即√a2+b2a2+b2

幅角：假设以逆时针为正方向，从xx轴正半轴到已知向量的转角的有向角叫做幅角

运算法则

加法：

因为在复平面中，复数可以被表示为向量，因此复数的加法与向量的加法相同，都满足平行四边形定则（就是上面那个）

乘法：

几何定义：复数相乘，模长相乘，幅角相加

代数定义：(a+bi)∗(c+di)(a+bi)∗(c+di)

=ac+adi+bci+bdi2=ac+adi+bci+bdi2

=ac+adi+bci−bd=ac+adi+bci−bd

=(ac−bd)+(bc+ad)i=(ac−bd)+(bc+ad)i

 

单位根

下文中，默认nn为22的正整数次幂

在复平面上，以原点为圆心，11为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的nn等分点为终点，做nn个向量，设幅角为正且最小的向量对应的复数为ωnωn，称为nn次单位根。

根据复数乘法的运算法则，其余n−1n−1个复数为ω2n,ω3n,…,ωnnωn2,ωn3,…,ωnn

注意ω0n=ωnn=1ωn0=ωnn=1（对应复平面上以xx轴为正方向的向量）

那么如何计算它们的值呢？这个问题可以由欧拉公式解决ωkn=cos k∗2πn+isink∗2πnωnk=cos⁡ k∗2πn+isin⁡k∗2πn

 

 

例如

 

图中向量ABAB表示的复数为88次单位根

单位根的幅角为周角的1n1n

在代数中，若zn=1zn=1，我们把zz称为nn次单位根

单位根的性质

  • ωkn=cosk2πn+isink2πnωnk=cos⁡k2πn+isin⁡k2πn（即上面的公式）
  • ω2k2n=ωknω2n2k=ωnk

证明：

ω2k2n=cos2k∗2π2n+isin2k∗2π2nω2n2k=cos⁡2k∗2π2n+isin⁡2k∗2π2n

=ωkn=ωnk

  • ωk+n2n=−ωknωnk+n2=−ωnk

ωn2n=cosn2∗2πn+isinn2∗2πnωnn2=cos⁡n2∗2πn+isin⁡n2∗2πn

=cosπ+isinπ=cos⁡π+isin⁡π

=−1=−1

  • ω0n=ωnn=1ωn0=ωnn=1

 讲了这么多，貌似跟我们的正题没啥关系啊。。

 OK！各位坐稳了，前方高能！

回到顶部

快速傅里叶变换

我们前面提到过，一个nn次多项式可以被nn个点唯一确定。

那么我们可以把单位根的00到n−1n−1次幂带入，这样也可以把这个多项式确定出来。但是这样仍然是O(n2)O(n2)的呀！

我们设多项式A(x)A(x)的系数为(ao,a1,a2,…,an−1)(ao,a1,a2,…,an−1)

那么A(x)=a0+a1∗x+a2∗x2+a3∗x3+a4∗x4+a5∗x5+⋯+an−2∗xn−2+an−1∗xn−1A(x)=a0+a1∗x+a2∗x2+a3∗x3+a4∗x4+a5∗x5+⋯+an−2∗xn−2+an−1∗xn−1

将其下标按照奇偶性分类

A(x)=(a0+a2∗x2+a4∗x4+⋯+an−2∗xn−2)+(a1∗x+a3∗x3+a5∗x5+⋯+an−1∗xn−1)A(x)=(a0+a2∗x2+a4∗x4+⋯+an−2∗xn−2)+(a1∗x+a3∗x3+a5∗x5+⋯+an−1∗xn−1)

设

 

A1(x)=a0+a2∗x+a4∗x2+⋯+an−2∗xn2−1A1(x)=a0+a2∗x+a4∗x2+⋯+an−2∗xn2−1

A2(x)=a1+a3∗x+a5∗x2+⋯+an−1∗xn2−1A2(x)=a1+a3∗x+a5∗x2+⋯+an−1∗xn2−1

那么不难得到

A(x)=A1(x2)+xA2(x2)A(x)=A1(x2)+xA2(x2)

我们将ωkn(k<n2)ωnk(k<n2)代入得

A(ωkn)=A1(ω2kn)+ωknA2(ω2kn)A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)

=A1(ωkn2)+ωknA2(ωkn2)=A1(ωn2k)+ωnkA2(ωn2k)

同理，将ωk+n2nωnk+n2代入得

A(ωk+n2n)=A1(ω2k+nn)+ωk+n2n(ω2k+nn)A(ωnk+n2)=A1(ωn2k+n)+ωnk+n2(ωn2k+n)

=A1(ω2kn∗ωnn)−ωknA2(ω2kn∗ωnn)=A1(ωn2k∗ωnn)−ωnkA2(ωn2k∗ωnn)

=A1(ω2kn)−ωknA2(ω2kn)=A1(ωn2k)−ωnkA2(ωn2k)

 

大家有没有发现什么规律？

没错！这两个式子只有一个常数项不同！

那么当我们在枚举第一个式子的时候，我们可以O(1)O(1)的得到第二个式子的值

又因为第一个式子的kk在取遍[0,n2−1][0,n2−1]时，k+n2k+n2取遍了[n2,n−1][n2,n−1]

所以我们将原来的问题缩小了一半！

而缩小后的问题仍然满足原问题的性质，所以我们可以递归的去搞这件事情！

直到多项式仅剩一个常数项，这时候我们直接返回就好啦

 

时间复杂度：

不难看出FFT是类似于线段树一样的分治算法。

因此它的时间复杂度为O(nlogn)O(nlogn)

 

回到顶部

快速傅里叶逆变换

不要以为FFT到这里就结束了。

我们上面的讨论是基于点值表示法的。

但是在平常的学习和研究中很少用点值表示法来表示一个多项式。

所以我们要考虑如何把点值表示法转换为系数表示法，这个过程叫做傅里叶逆变换

 

(y0,y1,y2,…,yn−1)(y0,y1,y2,…,yn−1)为(a0,a1,a2,…,an−1)(a0,a1,a2,…,an−1)的傅里叶变换（即点值表示）

设有另一个向量(c0,c1,c2,…,cn−1)(c0,c1,c2,…,cn−1)满足

ck=n−1∑i=0yi(ω−kn)ick=∑i=0n−1yi(ωn−k)i

即多项式B(x)=y0,y1x,y2x2,…,yn−1xn−1B(x)=y0,y1x,y2x2,…,yn−1xn−1在ω0n,ω−1n,ω−2n,…,ω−(n−1)n−1ωn0,ωn−1,ωn−2,…,ωn−1−(n−1)处的点值表示

emmmm又到推公式时间啦

(c0,c1,c2,…,cn−1)(c0,c1,c2,…,cn−1)满足
ck=n−1∑i=0yi(ω−kn)ick=∑i=0n−1yi(ωn−k)i

=n−1∑i=0(n−1∑j=0aj(ωin)j)(ω−kn)i=∑i=0n−1(∑j=0n−1aj(ωni)j)(ωn−k)i

=n−1∑i=0(n−1∑j=0aj(ωjn)i)(ω−kn)i=∑i=0n−1(∑j=0n−1aj(ωnj)i)(ωn−k)i

=n−1∑i=0(n−1∑j=0aj(ωjn)i(ω−kn)i)=∑i=0n−1(∑j=0n−1aj(ωnj)i(ωn−k)i)

=n−1∑i=0n−1∑j=0aj(ωjn)i(ω−kn)i=∑i=0n−1∑j=0n−1aj(ωnj)i(ωn−k)i

=n−1∑i=0n−1∑j=0aj(ωj−kn)i=∑i=0n−1∑j=0n−1aj(ωnj−k)i

=n−1∑j=0aj(n−1∑i=0(ωj−kn)i)=∑j=0n−1aj(∑i=0n−1(ωnj−k)i)

 

设S(x)=∑n−1i=0xiS(x)=∑i=0n−1xi

将ωknωnk代入得

S(ωkn)=1+(ωkn)+(ωkn)2+…(ωkn)n−1S(ωnk)=1+(ωnk)+(ωnk)2+…(ωnk)n−1

当k!=0k!=0时

等式两边同乘ωknωnk得

ωknS(ωkn)=ωkn+(ωkn)2+(ωkn)3+…(ωkn)nωnkS(ωnk)=ωnk+(ωnk)2+(ωnk)3+…(ωnk)n

两式相减得

ωknS(ωkn)−S(ωkn)=(ωkn)n−1ωnkS(ωnk)−S(ωnk)=(ωnk)n−1

S(ωkn)=(ωkn)n−1ωkn−1S(ωnk)=(ωnk)n−1ωnk−1

S(ωkn)=(ωnn)k−1ωkn−1S(ωnk)=(ωnn)k−1ωnk−1

S(ωkn)=1−1ωkn−1S(ωnk)=1−1ωnk−1

观察这个式子，不难看出它分母不为0，但是分子为0

因此，当K!=0K!=0时，S(ωkn)=0S(ωnk)=0

那当k=0k=0时呢？

很显然，S(ω0n)=nS(ωn0)=n

 

继续考虑刚刚的式子

ck=n−1∑j=0aj(n−1∑i=0(ωj−kn)i)ck=∑j=0n−1aj(∑i=0n−1(ωnj−k)i)
当j≠kj≠k时，值为00
当j=kj=k时，值为nn
因此，
ck=nakck=nak
ak=cknak=ckn

这样我们就得到点值与系数之间的表示啦

 

回到顶部

理论总结

至此，FFT的基础理论部分就结束了。

我们来小结一下FFT是怎么成功实现的

 

首先，人们在用系数表示法研究多项式的时候遇阻

于是开始考虑能否用点值表示法优化这个东西。

然后根据复数的两条性质（这个思维跨度比较大）得到了一种分治算法。

最后又推了一波公式，找到了点值表示法与系数表示法之间转换关系。

 

emmmm

其实FFT的实现思路大概就是

系数表示法—>点值表示法—>系数表示法

引用一下远航之曲大佬的图

 

 

当然，在实现的过程中还有很多技巧

我们根据代码来理解一下

 

回到顶部

递归实现

递归实现的方法比较简单。

就是按找我们上面说的过程，不断把要求的序列分成两部分，再进行合并

在c++的STL中提供了现成的complex类，但是我不建议大家用，毕竟手写也就那么几行，而且万一某个毒瘤卡STL那岂不是很ＧＧ？

 

// luogu-judger-enable-o2
// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN = 4 * 1e6 + 10;
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();}
    while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();}
    return x * f;
}
const double Pi = acos(-1.0);
struct complex {
    double x, y;
    complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
} a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);} //不懂的看复数的运算那部分
void fast_fast_tle(int limit, complex *a, int type) {
    if (limit == 1) return ; //只有一个常数项
    complex a1[limit >> 1], a2[limit >> 1];
    for (int i = 0; i <= limit; i += 2) //根据下标的奇偶性分类
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    fast_fast_tle(limit >> 1, a1, type);
    fast_fast_tle(limit >> 1, a2, type);
    complex Wn = complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = complex(1, 0);
    //Wn为单位根，w表示幂
    for (int i = 0; i < (limit >> 1); i++, w = w * Wn) //这里的w相当于公式中的k
        a[i] = a1[i] + w * a2[i],
               a[i + (limit >> 1)] = a1[i] - w * a2[i]; //利用单位根的性质，O(1)得到另一部分
}
int main() {
    int N = read(), M = read();
    for (int i = 0; i <= N; i++) a[i].x = read();
    for (int i = 0; i <= M; i++) b[i].x = read();
    int limit = 1; while (limit <= N + M) limit <<= 1;
    fast_fast_tle(limit, a, 1);
    fast_fast_tle(limit, b, 1);
    //后面的1表示要进行的变换是什么类型
    //1表示从系数变为点值
    //-1表示从点值变为系数
    //至于为什么这样是对的，可以参考一下c向量的推导过程，
    for (int i = 0; i <= limit; i++)
        a[i] = a[i] * b[i];
    fast_fast_tle(limit, a, -1);
    for (int i = 0; i <= N + M; i++) printf("%d ", (int)(a[i].x / limit + 0.5)); //按照我们推倒的公式，这里还要除以n
    return 0;
}
递归版

 

update:递归版我本地是可以AC的，只要开大数组就可以了，但是交到洛谷上会WA0

 

 

woc?  脸好疼。。。。。。

咳咳。

速度什么的才不是关键呢？

关键是我们AC不了啊啊啊

表着急，AC不了不代表咱们的算法不对，只能说这种实现方法太low了

下面介绍一种更高效的方法

 

回到顶部

迭代实现

再盗一下那位大佬的图

 

观察一下原序列和反转后的序列？

聪明的你有没有看出什么显而易见的性质？

没错！

我们需要求的序列实际是原序列下标的二进制反转！

因此我们对序列按照下标的奇偶性分类的过程其实是没有必要的

 

这样我们可以O(n)O(n)的利用某种操作得到我们要求的序列，然后不断向上合并就好了

 

复制代码 
// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN = 1e7 + 10;
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();}
    while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();}
    return x * f;
}
const double Pi = acos(-1.0);
struct complex {
    double x, y;
    complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
} a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);} //不懂的看复数的运算那部分
int N, M;
int l, r[MAXN];
int limit = 1;
void fast_fast_tle(complex *A, int type) {
    for (int i = 0; i < limit; i++)
        if (i < r[i]) swap(A[i], A[r[i]]); //求出要迭代的序列
    for (int mid = 1; mid < limit; mid <<= 1) { //待合并区间的长度的一半
        complex Wn( cos(Pi / mid) , type * sin(Pi / mid) ); //单位根
        for (int R = mid << 1, j = 0; j < limit; j += R) { //R是区间的长度，j表示前已经到哪个位置了
            complex w(1, 0); //幂
            for (int k = 0; k < mid; k++, w = w * Wn) { //枚举左半部分
                complex x = A[j + k], y = w * A[j + mid + k]; //蝴蝶效应
                A[j + k] = x + y;
                A[j + mid + k] = x - y;
            }
        }
    }
}
int main() {
    int N = read(), M = read();
    for (int i = 0; i <= N; i++) a[i].x = read();
    for (int i = 0; i <= M; i++) b[i].x = read();
    while (limit <= N + M) limit <<= 1, l++;
    for (int i = 0; i < limit; i++)
        r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1) ) ;
    // 在原序列中 i 与 i/2 的关系是 ： i可以看做是i/2的二进制上的每一位左移一位得来
    // 那么在反转后的数组中就需要右移一位，同时特殊处理一下奇数
    fast_fast_tle(a, 1);
    fast_fast_tle(b, 1);
    for (int i = 0; i <= limit; i++) a[i] = a[i] * b[i];
    fast_fast_tle(a, -1);
    for (int i = 0; i <= N + M; i++)
        printf("%d ", (int)(a[i].x / limit + 0.5));
    return 0;
}
复制代码

 

作者：自为风月马前卒 个人博客http://attack204.com// 出处：http://zwfymqz.cnblogs.com/ 本文版权归作者和博客园共有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文连接，否则保留追究法律责任的权利。   标签: 快速傅里叶变换FFT 好文要顶 关注我 收藏该文 微信分享 自为风月马前卒
粉丝 - 654 关注 - 104
    +加关注 207 0       « 上一篇： HDU 5536 Chip Factory
» 下一篇： LOJ #108. 多项式乘法 posted @ 2018-02-11 18:53  自为风月马前卒  阅读(72101)  评论(82)  编辑  收藏  举报

标签：nk,nn,FFT,kn,a2,n2,转载,0n
From： https://www.cnblogs.com/zhouruoheng/p/18149200

相关文章

  • 转载Using Domain-Driven Design(DDD)in Golang
    转载自：https://dev.to/stevensunflash/using-domain-driven-design-ddd-in-golang-3ee5UsingDomain-DrivenDesign(DDD)inGolang#go#ddd#redis#postgresDomain-DrivenDesignpatternisthetalkofthetowntoday.Domain-DrivenDesign(DDD)isanapproachtosoft......
  • 【转载】Java函数式编程——为什么new Thread()中要用函数式编程
    面向对象过分强调“必须通过对象的形式来做事情”，而函数式思想则尽量忽略面向对象的复杂语法——强调做什么，而不是以什么形式做。面向对象的思想:做一件事情,找一个能解决这个事情的对象,调用对象的方法,完成事情.函数式编程思想:只要能获取到结果,谁去做的,怎么做的都不重要,......
  • 【转载】WPF中Binding使用StringFormat格式化字符串方法
    原文链接：https://www.cnblogs.com/xuliming/articles/StringFormat.htmlWPF中Binding使用StringFormat格式化字符串方法 货币格式<TextBlockText="{BindingPrice,StringFormat={}{0:C}}"/>//$123.46货币格式，一位小数<TextBoxText="{BindingPrice,Stri......
  • [转载]mtu值怎样设置才网速最快
    推荐看原文链接mtu值怎样设置才网速最快-百度经验背景：什么是mtu答：mtu本地大于网络mtu,传输就会拆包，mtu本地设置小于mtu,就发挥不了网络最大传输能力，合适最好测试最佳mtu值ping-l1500-fwww.baidu.com如果出现超时结果正在Pingwww.baidu.com具有1500字节的数据:......
  • 【转载】清空WPF自带样式，以及透明按钮
    原文：https://www.cnblogs.com/Cindys/archive/2012/09/11/2680501.html空样式按钮<Stylex:Key="EmptyButtonStyle"TargetType="Button">           <SetterProperty="Padding"Value="0"/>           <SetterProper......
  • 华为网络设备配置（转载）
    网络工程师考试案例分析中常用的配置可以分为四大类，如下所示。本文我们一起关注前两部分路由器和交换机常用配置命令。 常用配置命令分类【路由器常见配置】 路由器常见配置访问控制列表配置 路由器ACL配置RIP路由配置路由器RIP路由配置OSPF配置......
  • 网安博士生的收藏夹-转载
    出处：https://blog.csdn.net/u013648063/article/details/122048461计算机领域如何寻找idea：加州大学河滨分校钱志云老师写的计算机领域如何找idea的一篇文章。（需要科*上网才能访问）如何阅读一篇论文：来源于滑铁卢大学Keshav教授。分享了一种“Three-pass”的方法。分三遍去读，每一遍......
  • SQL SERVER锁（LOCK）知识及锁应用(转载）
    SQLSERVER锁（LOCK）知识及锁应用(转载）ChenSir°已于 2023-10-0915:10:48 修改阅读量2.1k收藏8点赞数1文章标签：数据库sqlserver于 2023-10-0915:07:09 首次发布  提示：这里所摘抄的关于锁的知识有的是不同sqlserver版本的，对应于特......
  • FFT 与 NTT 学习笔记
    【概念】点值：给定多项式\(f(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\)和\(x_1\simx_m\)，求\(f(x_1),f(x_2),\dots,f(x_m)\)。（\(m=n\)）求点值的算法一般是\(O(n^2)\)的，但对于特殊的多项式，可以做到更快。插值：给定\((x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\)，其中\(x_0\s......
  • 【转载】冲压过程仿真模拟及优化 —— 冲压仿真的方法分类PPT
    地址：https://www.renrendoc.com/paper/310415051.html......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99