关于FFT快速傅里叶变换

前一段时间荒废掉了……一直在纠结zkw大牛的数学归纳法的论文……结果看不懂……


一直很纠结FFT是怎么实现的

一直以为FFT是很巧妙的利用数论之类的东西做,所以每次看到这方面的资料就蛋疼,因为一般都画了个坐标系,然后就懵了……


昨天晚上无聊翻算导看到FFT这一章,就把它学了,发现算导有的时候讲细一点还是有好处的,至少看得懂……


首先要明确的一点是FFT是用来求多项式相乘的,高精度乘法只是其中一种应用而已(把基看做多项式中的x)

那么大致过程是多项式的系数表示->DFT(O(nlogn))->多项式的点值表示->相乘(O(n))->DFT逆变换(O(nlogn))->得到相乘后的多项式系数表示

一直到省选之前应该都比较忙……具体细节以后再来补好了……


我用了complex类,代码参考自7k+……//今天才知道有complex类……

要注意的是就算complex使用的是double(应该说正因为使用的是double……),基不能选的很大

在贡献了若干次wa以后得到的结论是基最好是100000,再高就会爆精度,哪怕用long double也是一样……

http://61.187.179.132/JudgeOnline/problem.php?id=2179

FFT模板题


//Lib #include<cstdio> #include<cstring> #include<cstdlib> #include<cmath> #include<ctime> #include<cctype> #include<climits> #include<iostream> #include<algorithm> #include<vector> #include<string> #include<queue> #include<complex> using namespace std; //Macro #define rep(i,a,b) for(int i=a,tt=b;i<=tt;++i) #define drep(i,a,b) for(int i=a,tt=b;i>=tt;--i) #define erep(i,e,x) for(int i=x;i;i=e[i].next) #define irep(i,x) for(__typedef(x.begin()) i=x.begin();i!=x.end();i++) #define read() (strtol(ipos,&ipos,10)) #define sqr(x) ((x)*(x)) #define pb push_back #define PS system("pause"); typedef long long ll; typedef pair<int,int> pii; typedef complex<long double> comp; const int oo=~0U>>1; const double inf=1e20; const double eps=1e-6; const double pi=acos(-1.0); string name="fft",in=".in",out=".out"; const int maxn=131073; const comp I=comp(0,1); const int limit=100000; //Var comp A[maxn],B[maxn],C[maxn]; int base[10],len,n,la,lb; long double d; char s[100000]; void Getnum(comp num[],int &ln) { scanf("%s",s); for(int i=len-1;i>=0;i-=5) { num[ln]=s[i]-'0'; rep(j,1,4)if(i>=j)num[ln]=num[ln].real()+(s[i-j]-'0')*base[j]; ln++; } } void FFT(comp num[],long double root) { comp w,x;int i,j,k,m,h; for(m=n;h=m/2,m>=2;m=h,root*=2) for(i=0;i<h;i++) for(w=exp(i*root*I),j=i;j<n;j+=m) k=j+h,x=num[j]-num[k],num[j]+=num[k],num[k]=w*x; for(i=0,j=1;j<n-1;j++) { for(k=n/2;k>(i^=k);k/=2); if(j<i)swap(num[i],num[j]); } } void Init() { base[0]=1;rep(i,1,4)base[i]=base[i-1]*10; scanf("%d",&len); Getnum(A,la);Getnum(B,lb); } void Work() { for(n=1;n<la+lb;n<<=1); FFT(A,2*pi/n); FFT(B,2*pi/n); rep(i,0,n-1)C[i]=A[i]*B[i]; FFT(C,-2*pi/n); rep(i,0,n-1) { C[i]=C[i].real()/n; C[i]=C[i].real()+d; d=(int)(C[i].real()/limit); C[i]=C[i].real()-d*limit; if(C[i].real()>99999.9)C[i]=0.0,d++; if(C[i].real()<0.1)C[i]=0.0; } n--;while(n>=0&&fabs(C[n].real())<eps)n--; printf("%d",(int)(C[n].real()+0.5)); drep(i,n-1,0)printf("%05d",(int)(C[i].real()+0.5)); puts(""); } int main() { // freopen((name+in).c_str(),"r",stdin); // freopen((name+out).c_str(),"w",stdout); Init(); Work(); // PS; return 0; }


  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值