前一段时间荒废掉了……一直在纠结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;
}