【快速傅里叶变换】【FFT】【WikiOI】【P3132】【高精度练习之超大整数乘法】

最新推荐文章于 2022-09-13 10:10:28 发布
最新推荐文章于 2022-09-13 10:10:28 发布
阅读量10w+ 收藏 48
点赞数 22
分类专栏： OI 一入OI深似海 文章标签： bzoj FFT 快速傅里叶变换 傅里叶变换
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/iamzky/article/details/22712347
版权

传送门：http://www.wikioi.com/problem/3123/

FFT，快速傅里叶变换，蒟蒻看别人的题解都太深奥，看不懂，好不容易学会，以蒟蒻的理解写给那些想学FFT却又找不到合适的资料的OIer，蒟蒻理解有限，难免有许多错误，请大家多多包涵。


快速傅里叶变换

百度的各种讲解都TM扯什么频率什么的，蒟蒻完全看不懂，后来认真看了看算导，获益匪浅，算导上讲的真心不赖，有很多内容都来自算导。

1.多项式

        多项式的两种表达方式：系数表达和点值表达

系数表达就是大家常用的表达方式，点值表达就像在这个多项式函数上取n个不同的点，这样就可以确定原多项式。

比如说二次函数需要3个点就可以确定,一次函数需要2个点，一个n次多项式需要n个点(n次多项式意思是有0..n-1次幂的多项式)

A(x)=x^2+2*x-1可以被表达为{  ( 0 , -1 ) , ( 1 , 2 ) , ( 2 , 7 )  }

加法和乘法：

        

B(x)=x^2-x+2  { ( 0 , 2 ) , ( 1 , 2 ) , ( 2 , 4 ) }

     C(x)=A(x)+B(x)=2x^2+x+1   { ( 0, 1) , ( 1 , 4 ) , ( 2, 11 ) }

     注意乘法需要2n个点 lz比较懒就不写了……


于是我们得到一个计算多项式的方法：

 2.n次单位复数根

     

有关复数根的性质可以百度到，不再赘述

http://baike.baidu.com/link?url=017EPfseoBwVxWpWPm5aunUn8x9dmRvioav9IubYLSKEGngK8_rDV2bd4PFCM8sJ

3.DFT&&FFT

使用单位根计算点值表达式叫DFT（离散傅里叶变换）复杂度n^2，FFT是其优化版复杂度nlogn


计算FFT的伪代码(好吧用的是python的高亮)

下划线代表的是下标,括号代表上标，for 循环的range是左闭右开的

FFT(a):
	n=a.length()
	if n==1:
		return a
	w_n=e^(pi*i/n)=complex(cos(2*pi/n),sin(2*pi/n))
	w=1
	a(0)=[a0,a2,...a_n-2]
	a(1)=[a1,a3,...a_n-1]
	y(0)=FFT(a(0))
	y(1)=FFT(a(1))
	for k in range(0,n/2):
		y_k=y_k(0)+w*y_k(1)
		y_k+n/2=y_k(0)-w*y_k(1)
		w=w*w_n
	return y

4.递归=>迭代??


FFT的for循环中有两次w_n^k*y_k(1)的计算，于是可以改写成这样

 

for k in range(0,n/2):
	t=w*y_k(1)
	y_k=y_k(0)+t
	y_k+n/2=y_k(0)-t
	w=w*w_n
#这一过程被称蝴蝶操作

观察每次按照奇偶位置分割所形成的树：

每个数和他二进制相反的位置互换！！

伪代码（算导给的真是……）

BIT-REVERSE-COPY(a,A):
	n=a.length()
	for k in range(0,n):
		A[rev(k)]=a_k
#算导说rev函数很好写，就没写……

于是我们给出FFT的迭代实现的伪代码： 

FFT(a):
	BIT-REVERSE-COPY(a,A)
	n=a.length()
	for s in range(1,log2(n)+1):
		m=2^s
		w_m=e^(2*pi*i/m)=complex(cos(2*pi*m),sin(2*pi*m))
		for k in range(0,n,m):
			w=1
			for j in range(0,m/2):
				t=w*A[k+j+m/2]
				u=A[k+j]
				A[k+j]=u+t
				A[k+j+m/2]=u-t
				w=w*w_m
	return A
差不多讲完了,最后给出C++代码，有一大部分是lz借鉴别人的Code，以后附上地址 

#include<bitset>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define N 400005
#define pi acos(-1.0) // PI值
using namespace std;
struct complex
{
	double r,i;
	complex(double real=0.0,double image=0.0){
		r=real;	i=image;
	}
	// 以下为三种虚数运算的定义
	complex operator + (const complex o){
		return complex(r+o.r,i+o.i);
	}
	complex operator - (const complex o){
		return complex(r-o.r,i-o.i);
	}
	complex operator * (const complex o){
		return complex(r*o.r-i*o.i,r*o.i+i*o.r);
	}
}x1[N],x2[N];
char a[N/2],b[N/2];
int sum[N]; // 结果存在sum里
int vis[N];
void brc(complex *a,int l){//原来神犇的二进制平摊反转置换太神看不懂，蒟蒻写了一个O(n)的…… 
	memset(vis,0,sizeof(vis));//O(logn)的在后面 
	for(int i=1;i<l-1;i++){
		int x=i,y=0;
		int m=(int)log2(l)+0.1;
		if(vis[x])continue;
		while(m--){
			y<<=1;
			y|=(x&1);
			x>>=1;
		}
		vis[i]=vis[y]=1;
		swap(a[i],a[y]);
	}	
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
				     		// 其中on==1时为DFT，on==-1为IDFT
{
	register int h,i,j,k;
	complex u,t; 
	brc(y,l); // 调用反转置换
	for(h=2;h<=l;h<<=1) // 控制层数
	{
		// 初始化单位复根
		complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
		for(j=0;j<l;j+=h) // 控制起始下标
		{
			complex w(1,0); // 初始化螺旋因子
			for(k=j;k<j+h/2;k++) // 配对
			{
				u=y[k];
				t=w*y[k+h/2];
				y[k]=u+t;
				y[k+h/2]=u-t;
				w=w*wn; // 更新螺旋因子
			} // 据说上面的操作叫蝴蝶操作…
		}
	}
	if(on==-1)	for(i=0;i<l;i++)	y[i].r/=l; // IDFT
}
/* 
void fft2(complex *a,int s,int t){//蒟蒻自己写的递归版FFT，不保证正确 ,代码内部有未定义变量 
	if((n>>t)==1)return;//s记录起始,t记录深度,调用时应从0开始 
	fft(a,s,t+1);
	fft(a,s+(1<<t),t+1);
	for(int i=0;i<(n>>(t+1));i++){
		p=(i<<(t+1))+s;
		wt=w[i<<t]*a[p+(1<<t)];
		tt[i]=a[p]+wt;
		tt[i+(n>>(t+1))]=a[p]-wt;
	}
	for(i=0;i<(n>>t);i++)a[(i<<t)+s]=tt[i];
}*/
int main(void)
{
	int l1,l2,l;
	register int i;
	while(scanf("%s%s",a,b)!=EOF)
	{
		l1=strlen(a);
		l2=strlen(b);
		l=1;
		while(l<l1*2 || l<l2*2)	l<<=1; // 将次数界变成2^n
					      				// 配合二分与反转置换
		for(i=0;i<l1;i++) // 倒置存入
		{
			x1[i].r=a[l1-i-1]-'0';
			x1[i].i=0.0;
		}
		for(;i<l;i++)	x1[i].r=x1[i].i=0.0;
		// 将多余次数界初始化为0
		for(i=0;i<l2;i++)
		{
			x2[i].r=b[l2-i-1]-'0';
			x2[i].i=0.0;
		}
		for(;i<l;i++)	x2[i].r=x2[i].i=0.0;
		fft(x1,l,1); // DFT(a)
		fft(x2,l,1); // DFT(b)
		for(i=0;i<l;i++)	x1[i]=x1[i]*x2[i]; // 点乘结果存入a
		fft(x1,l,-1); // IDFT(a*b)
		for(i=0;i<l;i++)	sum[i]=x1[i].r+0.5; // 四舍五入
		for(i=0;i<l;i++) // 进位
		{
			sum[i+1]+=sum[i]/10;
			sum[i]%=10;
		}
		l=l1+l2-1;
		while(sum[l]<=0 && l>0)	l--; // 检索最高位
		for(i=l;i>=0;i--)	putchar(sum[i]+'0'); // 倒序输出
		putchar('\n');
	}
	return 0;
}
/*void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
{
	register int i,j,k;
	for(i=1,j=l/2;i<l-1;i++)
	{
		if(i<j)	swap(y[i],y[j]); // 交换互为下标反转的元素
								// i<j保证只交换一次
		k=l/2;
		while(j>=k) // 由最高位检索，遇1变0，遇0变1，跳出
		{
			j-=k;
			k>>=1;
		}
		if(j<k)	j+=k;
	}
}*/

pyc神犇的写法，bzoj3527，zjoi2014 力，无限YM

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=1000010;
int n,N,L;
int rev[maxn];
int dig[maxn];
double p[maxn];
struct cp{
	double r,i;
	cp(double _r=0,double _i=0):
		r(_r),i(_i){}
	cp operator+(cp x){return cp(r+x.r,i+x.i);}
	cp operator-(cp x){return cp(r-x.r,i-x.i);}
	cp operator*(cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}
};
cp a[maxn],b[maxn],c[maxn],A[maxn],x,y;
void FFT(cp a[],int flag){
	for(int i=0;i<N;i++)A[i]=a[rev[i]];
	for(int i=0;i<N;i++)a[i]=A[i];
	for(int i=2;i<=N;i<<=1){
		cp wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));
		for(int k=0;k<N;k+=i){
			cp w(1,0);
			for(int j=0;j<i/2;j++){
				x=a[k+j];
				y=w*a[k+j+i/2];
				a[k+j]=x+y;
				a[k+j+i/2]=x-y;
				w=w*wn;  
			}
		}
	}
	if(flag==-1)for(int i=0;i<N;i++)a[i].r/=N;
}
double anss[maxn];
int main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++)scanf("%lf",&p[i]);
	for(L=0,N=1;N<n;N<<=1,L++);L++;N<<=1;
	for(int i=0;i<N;i++){
		int len=0;
		for(int t=i;t;t>>=1)dig[len++]=t&1;
		for(int j=0;j<L;j++)rev[i]=rev[i]*2+dig[j];
	}
	for(int i=0;i<n;i++)a[i]=cp(p[i],0);
	for(int i=1;i<n;i++)b[i]=cp(1.0/i/i,0);
	FFT(a,1);FFT(b,1);
	for(int i=0;i<N;i++)c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<n;i++)anss[i]=c[i].r;
	memset(a,0,sizeof(a));
	memset(b,0,sizeof(b));
	for(int i=0;i<n;i++)a[i]=cp(p[n-i-1],0);
	for(int i=1;i<n;i++)b[i]=cp(1.0/i/i,0);
	FFT(a,1);FFT(b,1);
	for(int i=0;i<N;i++)c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<n;i++)anss[i]-=c[n-i-1].r;
	for(int i=0;i<n;i++)
		printf("%.9f\n",anss[i]);
	return 0;
}


重新过了一遍高精乘

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=1e6+10;
struct cp{
	double r,i;
	cp(double _r=0,double _i=0):
		r(_r),i(_i){}
	cp operator+(cp x){return cp(r+x.r,i+x.i);}
	cp operator-(cp x){return cp(r-x.r,i-x.i);}
	cp operator*(cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}
};
cp a[maxn],b[maxn],A[maxn],x,y,c[maxn];
char s1[maxn],s2[maxn];
int sum[maxn],a1[maxn],a2[maxn],dig[maxn];
int len1,len2,rev[maxn],N,L;
void FFT(cp a[],int flag){
	for(int i=0;i<N;i++)A[i]=a[rev[i]];
	for(int i=0;i<N;i++)a[i]=A[i];
	for(int i=2;i<=N;i<<=1){
		cp wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));
		for(int k=0;k<N;k+=i){
			cp w(1,0);
			for(int j=k;j<k+i/2;j++){
				x=a[j];
				y=a[j+i/2]*w;
				a[j]=x+y;
				a[j+i/2]=x-y;
				w=w*wn;
			}
		}
	}
	if(flag==-1)for(int i=0;i<N;i++)a[i].r/=N;
}
int main(){
	scanf("%s%s",s1,s2);
	len1=strlen(s1);
	len2=strlen(s2);
	for(N=1,L=0;N<max(len1,len2);N<<=1,L++);N<<=1;L++;
	for(int i=0;i<N;i++){
		int len=0;
		for(int t=i;t;t>>=1)dig[len++]=t&1;
		for(int j=0;j<L;j++)rev[i]=(rev[i]<<1)|dig[j];
	}
	for(int i=0;i<len1;i++)a1[len1-i-1]=s1[i]-'0';
	for(int i=0;i<len2;i++)a2[len2-i-1]=s2[i]-'0';
	for(int i=0;i<N;i++)a[i]=cp(a1[i]);
	for(int i=0;i<N;i++)b[i]=cp(a2[i]);
	FFT(a,1);FFT(b,1);
	for(int i=0;i<N;i++)c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<N;i++)sum[i]=c[i].r+0.5;
	for(int i=0;i<N;i++){
		sum[i+1]+=sum[i]/10;
		sum[i]%=10;
	}
	int l=len1+len2-1;
	while(sum[l]==0&&l>0)l--;
	for(int i=l;i>=0;i--)
	putchar(sum[i]+'0');
	putchar('\n');
	return 0;
}


iamzky
关注
  • 22
    点赞
  • 踩
  • 48
    收藏
    觉得还不错? 一键收藏
  • 10
    评论
专栏目录
十分简明易懂的FFT（快速傅里叶变换）
路人黑的纸巾
08-07 59万+
FFT有什么用 快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。 FFT（Fast Fouri...
数字信号处理-快速傅里叶变换FFT实验报告
05-06
西安交通大学数字信号处理-快速傅里叶变换FFT实验报告
FFT.rar_FFT 变换_c fft_傅立叶变换_快速傅立叶
09-24
实现快速傅立叶变换，此源码最大的特点是简捷明了。易于理解
快速傅立叶变换FFT频谱分析程序.pdf
01-11
快速傅立叶变换FFT频谱分析程序.pdf
快速傅里叶变换FFT的C语言实现及应用_fft_快速傅里叶变换_
09-29
根据图基和库利提出的kfft理论，采用c语言将其编制而成
FFT的详细解释，相信你看了就明白了
teresa_zp的专栏
12-28 11万+
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道F.........
FFT（快速傅里叶变换）
weixin_42638946的博客
10-17 1万+
FFT 1. FFT原理 原理 FFT（Fast Fourier Transformation），中文名快速傅里叶变换，用来加速多项式乘法。 朴素高精度乘法时间复杂度是O(n2)O(n ^ 2)O(n2)，n是数据位数；但FFT能在O(n×log(n))O(n \times log(n))O(n×log(n))的时间内解决。 假设给定一个多项式： A(x)=a0+a1x+a2x2+...an−1xn−1 A(x) = a_0 + a_1 x + a_2 x^2 + ... a_{n-1}
FFT简介小结
知识搬运工的博客
04-18 8611
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 ...
FFT算法讲解——麻麻我终于会FFT了！
热门推荐
谁能代表肯德基
03-12 15万+
FFT——快速傅里叶变换 这块不写东西空荡荡的，我决定还是把FFT的定义给贴上吧 FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。 这三段话其实一点用也没有 FFT是干什么的 FFT在算法竞赛中就有一个用途：加速多项式乘法 ...
实验4 快速傅立叶变换(FFT).doc
07-06
DSP实验4 快速傅立叶变换(FFT)
GDKOI 2017 参赛总结
lkb 的小屋
02-23 2652
陆陆续续也写了一个星期，希望此次比赛的失利和对问题的分析能帮助我们以后在类似的赛事中避免可能出现的问题。
【LOJ】#3119. 「CTS2019 | CTSC2019」随机立方体
weixin_30588675的博客
05-25 71
题解 用容斥，算至少K个极大值的方案数 我们先钦定每一维的K个数出来，然后再算上排列顺序是 \(w_{k} = \binom{n}{k}\binom{m}{k}\binom{l}{k}(k!)^3\) 然后有\((n - k)(m - k)(l - k)\)是可以随便填的 设\(all = nml,v_k = nml - (n - k)(m - k)(l - k)\) 设剩下的数填的方案...
bzoj2111-dp/Lucas定理
abby6688的博客
08-15 71
1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #include<string> 5 #include<cstring> 6 #include<algorithm> 7 #include<iomanip&g...
FFT(傅里叶快速变换，详细讲解+推导) 每日一遍，算法再见！
TheWayForDream
02-22 1万+
FFT详细推导FFT(傅里叶快速变换)一.前置知识1.复数和单位根2.单位根的三个引理3.多项式二.FFT(快速傅里叶变换推导)三.IFFT四.FFT求解多项式乘积模板代码1.递归版2.非递归版(这个更快，省去了递归时间)五.视频资源 FFT(傅里叶快速变换) FFT在实际工程中有着非常的广泛，尤其是在信号领域，在ACM算法竞赛领域主要可以用来快速计算多项式的乘积 一.前置知识 1.复数和单位根 有人会觉得FFT怎么会用到复数的知识，想必是有点古怪，后面在推导过程中会介绍到它的用处。 2.单位根的三个引
快速傅里叶变换——FFT
upupyon996deqing的博客
09-13 2万+
2，用一个公式来表示一下，这个公式的意义就是用来分解信号的，因为之前说FFT使用一个个正弦波来表示时域信号的，所以我们要确定正弦信号和时域信号的相关性，并且这个正弦波也要匹配块长度。上面的这个图就是加上传函数之后出现的FFT频谱图，但是可以看到激波的Bin比较宽，出现这个线性的原因是因为和快长度有关，如果提高块长度，那么分辨率就会变高，那么就会出现一个比较窄的Bin。其中，块长度也就是样本的数量跟FFT采集的周期有直接的关系，块长度对频率分辨率有直接的影响，假设采样的频率是48kHZ,块长度是512。
快速傅里叶变换 （FFT）基础
weixin_46118768的博客
07-15 9813
FFT
BZOJ2179: FFT快速傅立叶 FFT实现高精度乘法
最新发布
05-14
这道题目是要求用 FFT 实现高精度乘法。我们可以先将两个高精度数转换为多项式，然后对这两个多项式进行相乘，最后将结果转换为高精度数输出即可。 具体来说，将两个高精度数 $a$ 和 $b$ 分别转换为多项式 $A(x)$ 和 $B(x)$： $$ A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} \\ B(x) = b_0 + b_1x + b_2x^2 + \cdots + b_{n-1}x^{n-1} $$ 其中 $n$ 是 $a$ 和 $b$ 的位数，为了方便计算，我们要将 $n$ 扩展到 $2^k$，并在低位补 0。即： $$ A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} + 0x^n + \cdots + 0x^{2^k-1} \\ B(x) = b_0 + b_1x + b_2x^2 + \cdots + b_{n-1}x^{n-1} + 0x^n + \cdots + 0x^{2^k-1} $$ 接下来我们需要用 FFT 对 $A(x)$ 和 $B(x)$ 进行多项式乘法，得到多项式 $C(x) = A(x) \times B(x)$。然后将 $C(x)$ 转换为高精度数输出即可。 注意，由于 FFT 实现的是循环卷积，而我们要求的是线性卷积，因此在进行 FFT 前需要对 $A(x)$ 和 $B(x)$ 进行点值扩展，即令 $A'(x) = A(x) \times x^k$，$B'(x) = B(x) \times x^k$，其中 $k = 2^k - n$。这样得到的 $C'(x) = A'(x) \times B'(x)$ 在进行 FFT 后，其前 $n$ 项即为 $C(x)$ 的系数项，可以直接转换为高精度数输出。 代码实现如下：

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 10
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值