1. 程式人生 > >【模板】快速傅立葉變換(FFT)(BZOJ2179)

【模板】快速傅立葉變換(FFT)(BZOJ2179)

Description

給出兩個n位10進位制整數x和y,你需要計算x*y。

Input

第一行一個正整數n。 第二行描述一個位數為n的正整數x。 第三行描述一個位數為n的正整數y。

Output

輸出一行,即x*y的結果。

Code

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<complex>
#define cp complex <double>
using namespace std;
const double pi=acos(-1);
int l,n,res[1000010];
cp a[1000010],b[1000010],arr[1000010],inv[1000010];
char str1[1000010],str2[1000010];
void init()
{
    for (int i=0;i<n;i++)
    {
        arr[i]=cp(cos(2*pi*i/n),sin(2*pi*i/n));
        inv[i]=conj(arr[i]);
    }
}
void FFT(cp *a,cp *arr)
{
    int lim=0;
    while ((1<<lim)<n) lim++;
    for (int i=0;i<n;i++)
    {
        int t=0;
        for (int j=0;j<lim;j++)
            if ((i>>j) & 1) t|=1<<(lim-j-1);
        if (i<t) swap(a[i],a[t]);
    }
    for (int l=2;l<=n;l*=2)
    {
        int m=l/2;
        for (cp *buf=a;buf!=a+n;buf+=l)
            for (int i=0;i<m;i++)
            {
                cp t=arr[n/l*i]*buf[i+m];
                buf[i+m]=buf[i]-t;
                buf[i]+=t;
            }
    }
}
int main()
{
    scanf("%d",&l);
    scanf("%s",str1);scanf("%s",str2);
    n=1;while (n<2*l) n<<=1;
    init();
    for (int i=0;i<l;i++) a[l-i-1].real(str1[i]-'0');
    for (int i=0;i<l;i++) b[l-i-1].real(str2[i]-'0');
    FFT(a,arr),FFT(b,arr);
    for (int i=0;i<n;i++) a[i]*=b[i];
    FFT(a,inv);
    for (int i=0;i<n;i++)
    {
        res[i]+=floor(a[i].real()/n+0.5);
        res[i+1]+=res[i]/10;
        res[i]%=10;
    }
    for (int i=res[2*l-1]?2*l-1:2*l-2;i>=0;i--)
        putchar('0'+res[i]);
    puts("");
    return 0;
}