2
4
2016
0

【bzoj2179】FFT快速傅里叶

2179: FFT快速傅立叶

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 2196  Solved: 1098
[Submit][Status][Discuss]

Description

给出两个n位10进制整数x和y,你需要计算x*y。

Input

第一行一个正整数n。 第二行描述一个位数为n的正整数x。 第三行描述一个位数为n的正整数y。

Output

输出一行,即x*y的结果。

Sample Input

1
3
4
 

Sample Output

12
数据范围:
n<=60000
 

HINT

我们发现高精度乘法的每一位是满足卷积的形式的,所以我们先fft搞出c数组,然后再一位一位进位上去就可以了。

代码:

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
double pi=acos(-1);
int n,m,l,c[500010];
char s[100010];
struct cp{double x,y;};
cp a[500010],b[500010],cur[500010];
cp operator *(cp x,cp y){return (cp){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x};}
cp operator +(cp x,cp y){return (cp){x.x+y.x,x.y+y.y};}
cp operator -(cp x,cp y){return (cp){x.x-y.x,x.y-y.y};}
void fft(cp *a,int n,int fl)
{
    for (int i=n>>1,j=1;j<n;j++)
    {
        if (i<j) swap(a[i],a[j]);
        int k=n>>1;
        for (;k&i;i^=k,k>>=1);
        i^=k;
    }
    for (int m=2;m<=n;m<<=1)
    {
        cp w=(cp){cos(2*pi*fl/m),sin(2*pi*fl/m)};
        cur[0]=(cp){1,0};
        for (int i=1;i<m;i++) cur[i]=cur[i-1]*w;
        for (int i=0;i<n;i+=m)
            for (int j=i;j<i+(m>>1);j++)
            {
                cp u=a[j],v=a[j+(m>>1)]*cur[j-i];
                a[j]=u+v;
                a[j+(m>>1)]=u-v;
            }
    }
}
int main()
{
    scanf("%d",&n);
    scanf("%s",s);
    for (int i=0;i<n;i++)
        a[i].x=s[n-i-1]-'0';
    scanf("%s",s);
    for (int i=0;i<n;i++)
        b[i].x=s[n-i-1]-'0';
    m=1;while (m<=n*2) m*=2;
    fft(a,m,1);fft(b,m,1);
    for (int i=0;i<=m;i++) a[i]=a[i]*b[i];
    fft(a,m,-1);
    for (int i=0;i<=m;i++)
        c[i]=(int)(a[i].x/m+0.5);
    for (int i=0;i<=m;i++)
    {
        c[i+1]+=c[i]/10;
        c[i]%=10;
    }
    l=m;while (!c[l]) l--;
    for (int i=l;i>=0;i--) printf("%d",c[i]);printf("\n");
}
Category: FFT | Tags: | Read Count: 489

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com