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;
}