bzoj2179 FFT快速傅立叶

来源:互联网 发布:数龙网络 编辑:程序博客网 时间:2024/04/30 20:46

Description 给出两个n位10进制整数x和y,你需要计算x*y。 Input 第一行一个正整数n。
第二行描述一个位数为n的正整数x。 第三行描述一个位数为n的正整数y。 Output 输出一行,即x*y的结果。

FFT。

#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>using namespace std;const double pi=acos(-1.0);int rd(){    char c=getchar();    while (c<'0'||c>'9') c=getchar();    return c-'0';}struct com{    double a,b;    com operator + (const com &c) const    {        return (com){a+c.a,b+c.b};    }    com operator - (const com &c) const    {        return (com){a-c.a,b-c.b};    }    com operator * (const com &c) const    {        return (com){a*c.a-b*c.b,a*c.b+b*c.a};    }}a[200010],b[200010],p[200010];int n,l,t,ans[200010];int rev(int x){    int i,ret=0;    for (i=0;i<t;i++)        ret|=((x>>i)&1)<<t-i-1;    return ret;}void fft(com *a,int c){    int i,j,k,w;    com x,y;    for (i=0;i<l;i++)        if (rev(i)>i)            swap(a[i],a[rev(i)]);    for (i=2;i<=l;i<<=1)        for (j=0;j<l;j+=i)        {            w=0;            for (k=j;k<j+(i>>1);k++)            {                x=a[k];                y=a[k+(i>>1)]*p[w];                a[k]=x+y;                a[k+(i>>1)]=x-y;                w+=c*(l/i);                if (w<0) w+=l;            }        }}int main(){/*  freopen("in.txt","r",stdin);*/    int i;    scanf("%d",&n);    for (i=1;i<=n;i++)        a[n-i].a=rd();    for (i=1;i<=n;i++)        b[n-i].a=rd();    l=0;    while ((1<<t)<n*2) t++;    l=1<<t;    p[0]=(com){1,0};    p[1]=(com){cos(2*pi/l),sin(2*pi/l)};    for (i=2;i<l;i++)    {        p[i]=p[i/2]*p[i/2];        if (i&1) p[i]=p[i]*p[1];    }    fft(a,1);    fft(b,1);    for (i=0;i<l;i++)        a[i]=a[i]*b[i];    fft(a,-1);    for (i=0;i<l;i++)        a[i].a/=l;    for (i=0;i<l;i++)    {        ans[i]=ans[i]+a[i].a+0.5;        ans[i+1]=ans[i]/10;        ans[i]%=10;    }    for (i=l-1;i>=0&&!ans[i];i--);    for (;i>=0;i--) putchar(ans[i]+'0');}
1 0
原创粉丝点击