[BZOJ2179]FFT快速傅立叶

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
,B
;
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
,R
,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;
}