HDU1402 FFT高精度乘法模板题

#include<bits/stdc++.h>
using namespace std;
const int N = 50005*4;
const double PI = acos(-1);
typedef complex <double> cp;
char sa[N], sb[N];
int n = 1, lena, lenb, res[N];
cp a[N], b[N], DFT[N], IDFT[N];
void init()
{
    for(int i = 0; i < n; i++)
    {
        DFT[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n));
        IDFT[i] = conj(DFT[i]);
    }
}
void fft(cp *a, cp *wn)
{
    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]); // i < t 的限制使得每对点只被交换一次（否则交换两次相当于没交换）
    }
    for(int l = 2; l <= n; l *= 2)
    {
        int m = l / 2;
        for(cp *p = a; p != a + n; p += l)
            for(int i = 0; i < m; i++)
            {
                cp t = wn[n / l * i] * p[i + m];
                p[i + m] = p[i] - t;
                p[i] += t;
            }
    }
}
int main()
{

    while(scanf("%s%s",sa,sb)==2)
    {
        n=1;
        memset(res,0,sizeof(res));
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        lena = strlen(sa), lenb = strlen(sb);
        while(n < lena + lenb) n *= 2;
        for(int i = 0; i < lena; i++)
            a[i].real(sa[lena-1-i]-'0');
        for(int i = 0; i < lenb; i++)
            b[i].real(sb[lenb - 1 - i] - '0');
        init();
//求DFT
        fft(a,DFT);
        fft(b,DFT);
        for(int i = 0; i < n; i++)
            a[i] = a[i]*b[i];
//求IDFT
        fft(a,IDFT);

        for(int i = 0; i < n; i++)
        {
            res[i] += (int)(a[i].real() / n + 0.5);
            res[i + 1] += res[i] / 10;
            res[i] %= 10;
        }
        n = lena+lenb-1;
        while(res[n] <= 0 && n > 0)n--;
        for(int i = n; i >= 0; i--)
            printf("%c",res[i]+'0');
        printf("
");
    }
    return 0;
}