FFT之大数乘法

 #include <iostream>

 #include <stdio.h>

 #include <cmath>

 #include <algorithm>

 #include <cstring>

 #include <vector>

 using namespace std;

 #define N 50500*2

 const double PI = acos(-1.0);

 struct Vir

 {

     double re, im;

     Vir(double _re = ., double _im = .) :re(_re), im(_im){}

     Vir operator*(Vir r) { return Vir(re*r.re - im*r.im, re*r.im + im*r.re); }

     Vir operator+(Vir r) { return Vir(re + r.re, im + r.im); }

     Vir operator-(Vir r) { return Vir(re - r.re, im - r.im); }

 };

 void bit_rev(Vir *a, int loglen, int len)

 {

     for (int i = ; i < len; ++i)

     {

         int t = i, p = ;

         for (int j = ; j < loglen; ++j)

         {

             p <<= ;

             p = p | (t & );

             t >>= ;

         }

         if (p < i)

         {

             Vir temp = a[p];

             a[p] = a[i];

             a[i] = temp;

         }

     }

 }

 void FFT(Vir *a, int loglen, int len, int on)

 {

     bit_rev(a, loglen, len);

     for (int s = , m = ; s <= loglen; ++s, m <<= )

     {

         Vir wn = Vir(cos( * PI*on / m), sin( * PI*on / m));

         for (int i = ; i < len; i += m)

         {

             Vir w = Vir(1.0, );

             for (int j = ; j < m / ; ++j)

             {

                 Vir u = a[i + j];

                 Vir v = w*a[i + j + m / ];

                 a[i + j] = u + v;

                 a[i + j + m / ] = u - v;

                 w = w*wn;

             }

         }

     }

     if (on == -)

     {

         for (int i = ; i < len; ++i) a[i].re /= len, a[i].im /= len;

     }

 }

 char a[N * ], b[N * ];

 Vir pa[N * ], pb[N * ];

 int ans[N * ];

 int main()

 {

     while (scanf("%s%s", a, b) != EOF)

     {

         int lena = strlen(a);

         int lenb = strlen(b);

         int n = , loglen = ;

         while (n < lena + lenb) n <<= , loglen++;

         for (int i = , j = lena - ; i < n; ++i, --j)

             pa[i] = Vir(j >=  ? a[j] - '' : ., .);

         for (int i = , j = lenb - ; i < n; ++i, --j)

             pb[i] = Vir(j >=  ? b[j] - '' : ., .);

         for (int i = ; i <= n; ++i) ans[i] = ;

         FFT(pa, loglen, n, );

         FFT(pb, loglen, n, );

         for (int i = ; i < n; ++i)

             pa[i] = pa[i] * pb[i];

         FFT(pa, loglen, n, -);

         for (int i = ; i < n; ++i) ans[i] = pa[i].re + 0.5;

         for (int i = ; i<n; ++i) ans[i + ] += ans[i] / , ans[i] %= ;

         int pos = lena + lenb - ;

         for (; pos> && ans[pos] <= ; --pos);

         for (; pos >= ; --pos) printf("%d", ans[pos]);

         puts("");

     }

     return ;

 }