[BZOJ2179]FFT快速傅立叶

来源：互联网发布：知己而知彼论坛编辑：程序博客网时间：2024/04/30 20:08

FFT快速傅立叶

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

Solution
一个n位大整数可以视作一个n项n-1次的多项式代入10的结果
整数乘法符合卷积性质

Code

#include <bits/stdc++.h>using namespace std;const double pi=acos(-1.);const int N=301000;struct Comp{    double x,y;    Comp(double _x=0.,double _y=0.):x(_x),y(_y){}}A[N],B[N];inline Comp operator +(const Comp &a,const Comp &b){return Comp(a.x+b.x,a.y+b.y);}inline Comp operator -(const Comp &a,const Comp &b){return Comp(a.x-b.x,a.y-b.y);}inline Comp operator *(const Comp &a,const Comp &b){return Comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}int n,m,C[N],R[N],L;void read(Comp *A){    for(int i=n-1;i>=0;i--){        char ch; while(!isdigit(ch=getchar()));        A[i].x=ch-'0';    }}void FFT(Comp *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){        Comp wn=Comp(cos(pi/i),f*sin(pi/i));        for(int j=0,p=i<<1;j<n;j+=p){            Comp w=Comp(1.,0.);            for(int k=0;k<i;k++,w=w*wn){                Comp A0=a[j+k],A1=w*a[j+k+i];                a[j+k]=A0+A1;a[j+k+i]=A0-A1;            }        }    }}int main(){    scanf("%d",&n);read(A);read(B);    for(m=2*n,n=1;n<=m;n<<=1,L++);    for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-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<n;i++)C[i]=(int)round(A[i].x/n);    for(int i=0,d=0;i<=m;i++){C[i]+=d;d=C[i]/10;C[i]%=10;}    for(;!C[m];m--);    for(;m>=0;m--)printf("%d",C[m]);puts("");    return 0;    }

0 0