码迷,mamicode.com
首页 > 其他好文 > 详细

快速傅里叶变换(FFT)

时间:2014-11-26 20:32:19      阅读:285      评论:0      收藏:0      [点我收藏+]

标签:blog   http   io   ar   os   sp   for   on   div   

埋了一天的算导就当我看懂了?

。。。

目前仅限于学到FFT计算多项式的卷积,什么频域什么东西的那些我都不懂。。。。

最让我感到神奇的是复数这个概念,,,好强大。。

复数的话我大概懂得这点?:

复数有实部和虚部,其中虚部的单位是$i=$,定义为$e=a+bi$,带$i$的是虚部

然后当虚部为0时,这个复数就是实数。。。(也就是说实数是复数的子集。。。。)

然后复数的运算是:当$x=a+bi, y=c+di$

$$x+y=(a+c)+(b+d)i$$

$$x-y=(a-c)+(b-d)i$$

$$x \times y=(a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i$$

除法怎么搞。。。。。。(以后补。。

然后本文不详细说fft怎么实现,我觉得如果入门的话看算导一定能看懂的,而且算导说的真的很简单。。。

本文只用来记录我之前的疑问和暂时还有疑问的地方:

为什么是用复数指数幂的形式?我的解释是,,,因为根据复数幂定义能提供良好的相消引理。。这个很简单吧。。

为什么复数幂的形式是$e^{iu}=cos(u)+isin(u)$呢?因为周期性?

为什么复数做x时,赋值到复数的实部?之前解释过了,因为实数是复数的子集。

为什么$DFT^{-1}$后取的也是实部?同上

(晚上回来补。。codevs:3123 高精度练习之超大整数乘法

先放模板吧。。

题目:

#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <iostream>
#include <algorithm>
#include <queue>
#include <set>
#include <map>
using namespace std;
typedef long long ll;
#define rep(i, n) for(int i=0; i<(n); ++i)
#define for1(i,a,n) for(int i=(a);i<=(n);++i)
#define for2(i,a,n) for(int i=(a);i<(n);++i)
#define for3(i,a,n) for(int i=(a);i>=(n);--i)
#define for4(i,a,n) for(int i=(a);i>(n);--i)
#define CC(i,a) memset(i,a,sizeof(i))
#define read(a) a=getint()
#define print(a) printf("%d", a)
#define dbg(x) cout << (#x) << " = " << (x) << endl
#define error(x) (!(x)?puts("error"):0)
inline const int getint() { int r=0, k=1; char c=getchar(); for(; c<‘0‘||c>‘9‘; c=getchar()) if(c==‘-‘) k=-1; for(; c>=‘0‘&&c<=‘9‘; c=getchar()) r=r*10+c-‘0‘; return k*r; }
#define rdm(x, i) for(int i=ihead[x]; i; i=e[i].next)

struct cp {
	double r, i;
	cp(double _r=0.0, double _i=0.0) : r(_r), i(_i) {}
	cp operator+ (const cp &x) const { return cp(r + x.r, i + x.i); }
	cp operator- (const cp &x) const { return cp(r - x.r, i - x.i); }
	cp operator* (const cp &x) const { return cp(r*x.r - i*x.i, r*x.i + i*x.r); }
};
const int N=1000005;
const double PI=acos(-1.0);
int rev[N];
cp A[N];
void DFT(cp *a, int n, int flag) {
	rep(i, n) A[rev[i]]=a[i];
	rep(i, n) a[i]=A[i];
	for(int m=2; m<=n; m<<=1) {
		cp wn(cos(2.0*PI/m*flag), sin(2.0*PI/m*flag));
		for(int i=0; i<n; i+=m) {
			cp w(1.0); int mid=m>>1;
			rep(j, mid) {
				cp u=a[i+j+mid]*w, t=a[i+j];
				a[i+j]=t+u;
				a[i+j+mid]=t-u;
				w=w*wn;
			}
		}
	}
	if(flag==-1) rep(i, n) a[i].r/=n;
}

char s[N];
void readin(cp *a, int &n) {
	scanf("%s", s);
	n=strlen(s);
	rep(i, n) a[i].r=s[n-i-1]-‘0‘;
}
void init(int &n) {
	int k=1, L=0;
	while(k<n) k<<=1, ++L;
	n=k;
	rep(i, n) {
		int t=i, ret=0; k=L;
		while(k--) ret<<=1, ret|=(t&1), t>>=1;
		rev[i]=ret;
	}
}
cp a[N], b[N];
int n, len, ans[N];
int main() {
	readin(a, n); len+=n;
	readin(b, n); len+=n;
	--len;
	init(len);
	DFT(a, len, 1); DFT(b, len, 1);
	rep(i, len) a[i]=a[i]*b[i];
	DFT(a, len, -1);
	rep(i, len) ans[i]=(int)(a[i].r+0.5);
	rep(i, len) ans[i+1]+=ans[i]/10, ans[i]%=10;
	++len;
	while(ans[len]==0 && len) --len;
	for3(i, len, 0) printf("%d", ans[i]);
	return 0;
}

 

 

快速傅里叶变换(FFT)

标签:blog   http   io   ar   os   sp   for   on   div   

原文地址:http://www.cnblogs.com/iwtwiioi/p/4123976.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!