HDU 1402 A * B Problem Plus

Problem Description

Calculate A * B.

Input

Each line will contain two integers A and B. Process to end of file.

Note: the length of each integer will not exceed 50000.

Output

For each case, output A * B in one line.

Sample Input

1

2

1000

2

Sample Output

2

2000

FFT的入门模板题

花了好几天在网上看了无数资料，终于算是理解了FFT到底做了什么，FFT为什么那么神奇了。

先推荐几个不错得资料

http://blog.jobbole.com/70549/ 这篇轻松有趣，讲解FFT的定义

http://tieba.baidu.com/p/2513502552#36771774998l 这篇介绍FFT的过程，简单易懂

http://wlsyzx.yzu.edu.cn/kcwz/szxhcl/kechenneirong/jiaoan/jiaoan3.htm 这个就是纯数学了

这里我就简单介绍一下FFT是如何来解决大数乘法的问题的吧。

首先你要知道有一个叫做卷积的东西，详细的可以看看百度百科什么的

​

然后卷积有线性卷积和循环卷积，所谓线性卷积，在离散的情况下，可以看做是两个关于x的多项式的乘积的系数。

而FFT是DFT的快速算法，是信号分析中的重要内容，可以用来快速计算循环卷积。

然后我们把线性卷积看作是长度很长的循环卷积，那么计算出来的循环卷积可以在数值上等于线性卷积。

接下来讲一下计算过程，需要计算的数组分别是a和b

即两个多项式 

a[0]x^(n-1)+a[1]x^(n-2)+...+a[n-2]x+a[n-1]

b[0]x^(m-1)+b[1]x^(m-2)+...+b[m-2]x+b[m-1]相乘，需要的循环卷积长度最少为n+m-1，证明上面的链接里有      

还有一点是FFT只能算2^k的循环卷积的，所以要把长度加到2的次方去。

然后对a数组和b数组分别求一次FFT，得到数组对应位置再相乘

最后再IFFT回来就是它们的循环卷积了。

这样对于高精度乘法来说，就是x=10的特殊情况了，只要按照上述过程就可以得到结果了。

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int maxn = 100005;
const double pi = acos(-1.0);
const int low(int x){ return x&-x; }//树状数组的lowbit，可以计算二进制最右边的1
int n, m;
char s1[maxn], s2[maxn];

class FFT
{
private:
  const static int maxn = 150000;//要注意长度是2^k方
  class Plural
  {
  public:
    double x, y;
    Plural(double x = 0.0, double y = 0.0) :x(x), y(y){}
    Plural operator +(const Plural &a)
    {
      return Plural(x + a.x, y + a.y);
    }
    Plural operator -(const Plural &a)
    {
      return Plural(x - a.x, y - a.y);
    }
    Plural operator *(const Plural &a)
    {
      return Plural(x*a.x - y*a.y, x*a.y + y*a.x);
    }
    Plural operator /(const double &u)
    {
      return Plural(x / u, y / u);
    }
  };//定义复数的相关运算
  Plural x[maxn], x1[maxn], x2[maxn];
  Plural y[maxn], y1[maxn], y2[maxn];
  int ans[maxn];
  int n, len;
public:
  int reverse(int x)
  {
    int ans = 0;
    for (int i = 1, j = n >> 1; j; i <<= 1, j >>= 1) if (x&i) ans |= j;
    return ans;
  }//数字倒序,FFT的初始步骤
  Plural w(double x, double y)
  {
    return Plural(cos(2 * pi * x / y), -sin(2 * pi * x / y));
  }
  void setx(char *s, char *S)
  {
    int len1 = 0, len2 = 0;
    for (int i = 0; s[i]; i++, len1++) x1[i].x = s[i] - '0', x1[i].y = 0;
    for (int i = 0; S[i]; i++, len2++) x2[i].x = S[i] - '0', x2[i].y = 0;
    for (len = n = len1 + len2 - 1; n != low(n); n += low(n));
    for (int i = len1; i < n; i++) x1[i] = Plural(0, 0);
    for (int i = len2; i < n; i++) x2[i] = Plural(0, 0);
  }
  void fft(Plural*x, Plural*y, int flag)
  {
    for (int i = 0; i < n; i++) y[i] = x[reverse(i)];
    for (int i = 1; i < n; i <<= 1)
    {
      Plural uu = w(flag, i + i);
      for (int j = 0; j < n; j += i + i)
      {
        Plural u(1, 0);
        for (int k = j; k < j + i; k++)
        {
          Plural a = y[k];
          //w(flag*(k - j), i + i) 可以去掉u和uu用这个代替，精度高些，代价是耗时多了
          Plural b = u * y[k + i];
          y[k] = a + b;
          y[k + i] = a - b;
          u = u*uu;
        }
      }
    }
    if (flag == -1) for (int i = 0; i < n; i++) y[i] = y[i] / n;
  }//1是FFT，-1是IFFT，答案数组是y数组
  void solve()
  {
    fft(x1, y1, 1); fft(x2, y2, 1);
    for (int i = 0; i < n; i++) y[i] = y1[i] * y2[i];
    fft(y, x, -1);
    for (int i = 0; i < n; i++) ans[i] = (int)(x[i].x + 0.5);//调整精度
    for (int i = n - 1; i; i--)
    {
      ans[i - 1] += ans[i] / 10;
      ans[i] %= 10;
    }
    for (int i = 0; i < len; i++)
    {
      if (ans[i] != 0)
      {
        for (; i < len; i++) printf("%d", ans[i]);
        break;
      }
      else if (i == len - 1) putchar('0');
    }
    putchar(10);
  }
}fft;

int main()
{
  while (scanf("%s%s", s1, s2) != EOF)
  {
    fft.setx(s1, s2);
    fft.solve();
  }
  return 0;
}