bzoj2179: FFT快速傅立叶

http://www.lydsy.com/JudgeOnline/problem.php?id=2179

 

FFT做高精乘

 

#include<cmath>
#include<cstdio>
#include<complex>

using namespace std;

#define N 60001
const int M=(1<<17)+10;

const double pi=acos(-1);

typedef complex<double> E;

int n;
char s[N];

int r[M];

E a[M],b[M];

int res[M];

void fft(E *a,int f)
{
    for(int i=0;i<n;++i) 
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<n;i<<=1)
    {
        E wn(cos(pi/i),f*sin(pi/i));
        for(int p=i<<1,j=0;j<n;j+=p)
        {
            E w(1,0);
            for(int k=0;k<i;++k,w*=wn)
            {
                E x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y; a[j+k+i]=x-y;
            }
        }
    }
}

int main()
{
    int m;
    scanf("%d",&n); 
    n--; m=n;
    scanf("%s",s);
    for(int i=n;i>=0;--i) a[n-i]=s[i]-'0';
    scanf("%s",s);
    for(int i=m;i>=0;--i) b[m-i]=s[i]-'0';
    m+=n; int bit=0;
    for(n=1;n<=m;n<<=1) bit++;
    for(int i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<bit-1);
    fft(a,1);
    fft(b,1);
    for(int i=0;i<n;++i) a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=m;++i) res[i]=(int)(a[i].real()/n+0.5);
    for(int i=0;i<=m;++i) res[i+1]+=res[i]/10,res[i]%=10;
    if(res[m+1]) m++;
    //while(m && !res[m]) --m;
    for(int i=m;i>=0;--i) printf("%d",res[i]);
}