1. 程式人生 > >【HDU1402】 【FFT求大數乘法】

【HDU1402】 【FFT求大數乘法】

描述:

A * B Problem Plus

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 18015    Accepted Submission(s): 3981


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
Author DOOM III
Recommend DOOM III
分析:

這題的資料量是5w, 也就是傳統意義上的n^2演算法是不可取的。這裡就用到了FFT加快計算

FFT一般的作用就是使得多項式乘法的複雜度降到nlogn。利用FFT可以快速求出迴圈卷積。

那麼卷積又是什麼樣一個東西。

訊號處理中的一個重要運算是卷積.初學卷積的時候,往往是在連續的情形,
  兩個函式f(x),g(x)的卷積,是∫f(u)g(x-u)du
  當然,證明卷積的一些性質並不困難,比如交換,結合等等,但是對於卷積運算的來處,初學者就不甚了了。
  
  其實,從離散的情形看卷積,或許更加清楚,
  對於兩個序列f[n],g[n],一般可以將其卷積定義為s[x]= ∑f[k]g[x-k]
  
  卷積的一個典型例子,其實就是初中就學過的多項式相乘的運算,
  比如(x*x+3*x+2)(2*x+5)
  一般計算順序是這樣,
  (x*x+3*x+2)(2*x+5)
  = (x*x+3*x+2)*2*x+(x*x+3*x+2)*5
  = 2*x*x*x+3*2*x*x+2*2*x+ 5*x*x+3*5*x+10
  然後合併同類項的係數,
  2 x*x*x
  3*2+1*5 x*x
  2*2+3*5 x
  2*5
  ----------
  2*x*x*x+11*x*x+19*x+10
  
  實際上,從線性代數可以知道,多項式構成一個向量空間,其基底可選為
  {1,x,x*x,x*x*x,...}
  如此,則任何多項式均可與無窮維空間中的一個座標向量相對應,
  如,(x*x+3*x+2)對應於
  (1 3 2),
  (2*x+5)對應於
  (2,5).
  
  線性空間中沒有定義兩個向量間的卷積運算,而只有加法,數乘兩種運算,而實際上,多項式的乘法,就無法線上性空間中說明.可見線性空間的理論多麼侷限了.
  但如果按照我們上面對向量卷積的定義來處理座標向量,
  (1 3 2)*(2 5)
  則有
  2 3 1
  _ _ 2 5
  --------
      2
  
  
  2 3 1
  _ 2 5
  -----
    6+5=11
  
  2 3 1
  2 5
  -----
  4+15 =19
  
  
  _ 2 3 1
  2 5
  -------
    10
  
   或者說,
  (1 3 2)*(2 5)=(2 11 19 10)
  
  回到多項式的表示上來,
  (x*x+3*x+2)(2*x+5)= 2*x*x*x+11*x*x+19*x+10
  
  似乎很神奇,結果跟我們用傳統辦法得到的是完全一樣的.
  換句話,多項式相乘,相當於係數向量的卷積.
  
  其實,琢磨一下,道理也很簡單,
  卷積運算實際上是分別求 x*x*x ,x*x,x,1的係數,也就是說,他把加法和求和雜合在一起做了。(傳統的辦法是先做乘法,然後在合併同類項的時候才作加法)
  以x*x的係數為例,得到x*x,或者是用x*x乘5,或者是用3x乘2x,也就是
  2 3 1
  _ 2 5
  -----
   6+5=11
  其實,這正是向量的內積.如此則,卷積運算,可以看作是一串內積運算.既然是一串內積運算,則我們可以試圖用矩陣表示上述過程。
  
  [ 2 3 1 0 0 0]
  [ 0 2 3 1 0 0]==A
  [ 0 0 2 3 1 0]
  [ 0 0 0 2 3 1]
  
  [0 0 2 5 0 0]' == x
  
  b= Ax=[ 2 11 19 10]'
  
  採用行的觀點看Ax,則b的每行都是一個內積。
  A的每一行都是序列[2 3 1]的一個移動位置。
  
  ---------
  
  顯然,在這個特定的背景下,我們知道,卷積滿足交換,結合等定律,因為,眾所周知的,多項式的乘法滿足交換律,結合律.在一般情形下,其實也成立.
  
  在這裡,我們發現多項式,除了構成特定的線性空間外,基與基之間還存在某種特殊的聯絡,正是這種聯絡,給予多項式空間以特殊的性質.
  
  在學向量的時候,一般都會舉這個例子,甲有三個蘋果,5個橘子,乙有5個蘋果,三個橘子,則共有幾個蘋果,橘子。老師反覆告誡,橘子就是橘子,蘋果就是蘋果,可不能混在一起。所以有(3,5)+(5,3)=(8,8).是的,橘子和蘋果無論怎麼加,都不會出什麼問題的,但是,如果考慮橘子乘橘子,或者橘子乘蘋果,這問題就不大容易說清了。
  
  又如複數,如果僅僅定義複數為數對(a,b),僅僅線上性空間的層面看待C2,那就未免太簡單了。實際上,只要加上一條(a,b)*(c,d)=(ac-bd,ad+bc)
  則情況馬上改觀,複變函式的內容多麼豐富多彩,是眾所周知的。
  
  另外,回想訊號處理裡面的一條基本定理,頻率域的乘積,相當於時域或空域訊號的卷積.恰好跟這裡的情形完全對等.這後面存在什麼樣的隱態聯絡,需要繼續參詳.
  
  從這裡看,高等的卷積運算其實不過是一種初等的運算的抽象而已.中學學過的數學裡面,其實還蘊涵著許多高深的內容(比如交換代數)。溫故而知新,斯言不謬.
  
  其實這道理一點也不復雜,人類繁衍了多少萬年了,但過去n多年,人們只知道男女媾精,乃能繁衍後代。精子,卵子的發現,生殖機制的研究,也就是最近多少年的事情。
  
  孔子說,道在人倫日用中,看來我們應該多用審視的眼光看待周圍,乃至自身,才能知其然,而知其所以然。

----------------------------------------------------------完畢------------------------------

然後我們就知道卷積大概的作用了。

那麼FFT本來是訊號裡面的東西,而我沒學過訊號。 所以看的也不怎麼懂。

大概就是對離散的訊號,先將其轉變為一些正弦函式,然後這些正弦函式疊加能構成這個離散訊號,但是這些正弦函式易於處理。處理完之後就可以再轉變回來。

兩個過程叫做DFT和IDFT。

如果令上面的x=10

那麼就可以把兩個大整數相乘看做是多項式乘法。

最後求出各系數後再進位即可

程式碼一:
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <map>
#include <queue>
#include <set>
#include <vector>
using namespace std;
#define L(x) (1 << (x))
const double PI = acos(-1.0);
const int Maxn = 133015;
double ax[Maxn], ay[Maxn], bx[Maxn], by[Maxn];
char sa[Maxn/2],sb[Maxn/2];
int sum[Maxn];
int x1[Maxn],x2[Maxn];

int revv(int x, int bits){
  int ret = 0;
  for (int i = 0; i < bits; i++){
    ret <<= 1;
    ret |= x & 1;
    x >>= 1;
  }
  return ret;
}

void fft(double * a, double * b, int n, bool rev){
  int bits = 0;
  while (1 << bits < n) ++bits;
  for (int i = 0; i < n; i++){
    int j = revv(i, bits);
    if (i < j)
      swap(a[i], a[j]), swap(b[i], b[j]);
  }
  for (int len = 2; len <= n; len <<= 1){
    int half = len >> 1;
    double wmx = cos(2 * PI / len), wmy = sin(2 * PI / len);
    if (rev) wmy = -wmy;
    for (int i = 0; i < n; i += len){
      double wx = 1, wy = 0;
      for (int j = 0; j < half; j++){
        double cx = a[i + j], cy = b[i + j];
        double dx = a[i + j + half], dy = b[i + j + half];
        double ex = dx * wx - dy * wy, ey = dx * wy + dy * wx;
        a[i + j] = cx + ex, b[i + j] = cy + ey;
        a[i + j + half] = cx - ex, b[i + j + half] = cy - ey;
        double wnx = wx * wmx - wy * wmy, wny = wx * wmy + wy * wmx;
        wx = wnx, wy = wny;
      }
    }
  }
    if (rev)
    {
        for (int i = 0; i < n; i++)
            a[i] /= n, b[i] /= n;
    }
}

int solve(int a[],int na,int b[],int nb,int ans[]){
  int len = max(na, nb), ln;
  for(ln=0; L(ln)<len; ++ln);
  len=L(++ln);
  for(int i = 0; i < len ; ++i){
    if (i >= na) ax[i] = 0, ay[i] =0;
    else ax[i] = a[i], ay[i] = 0;
  }
  fft(ax, ay, len, 0);
  for (int i = 0; i < len; ++i){
    if (i >= nb) bx[i] = 0, by[i] = 0;
    else bx[i] = b[i], by[i] = 0;
  }
  fft(bx, by, len, 0);
  for (int i = 0; i < len; ++i){
    double cx = ax[i] * bx[i] - ay[i] * by[i];
    double cy = ax[i] * by[i] + ay[i] * bx[i];
    ax[i] = cx, ay[i] = cy;
  }
  fft(ax, ay, len, 1);
  for(int i = 0; i < len; ++i)
    ans[i] = (int)(ax[i] + 0.5);
  return len;
}

int main(){
  int l1,l2,l;
  int i;
  while(gets(sa)){
    gets(sb);
    memset(sum, 0, sizeof(sum));
    l1 = strlen(sa);
    l2 = strlen(sb);
    for(i = 0; i < l1; i++)
      x1[i] = sa[l1-i-1]-'0';
    for(i = 0; i < l2; i++)
      x2[i] = sb[l2-i-1]-'0';
    l = solve(x1, l1, x2, l2, sum);
    for(i = 0; i<l || sum[i] >= 10; i++){ // 進位
      sum[i + 1] += sum[i] / 10;
      sum[i] %= 10;
    }
    l = i;
    while(sum[l] <= 0 && l>0)   l--; // 檢索最高位
    for(i = l; i >= 0; i--)   putchar(sum[i] + '0'); // 倒序輸出
    putchar('\n');
  }
  return 0;
}

乘法其實就是做線性卷積。

用DFT得方法可以求迴圈卷積,但是當迴圈卷積長度LN+M-1,就可以做線性卷積了。

使用FFT將兩個數列轉換成傅立葉域,在這的乘積就是時域的卷積。

給幾個學習的連結吧:

(得好好研究研究_(:зゝ∠)_)

程式碼二:

(bin神模板)

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <math.h>
using namespace std;

const double PI = acos(-1.0);
//複數結構體
struct complex{
  double r,i;
  complex(double _r = 0.0,double _i = 0.0){
    r = _r; i = _i;
  }
  complex operator +(const complex &b){
    return complex(r+b.r,i+b.i);
  }
  complex operator -(const complex &b){
    return complex(r-b.r,i-b.i);
  }
  complex operator *(const complex &b){
    return complex(r*b.r-i*b.i,r*b.i+i*b.r);
  }
};
/*
 * 進行FFT和IFFT前的反轉變換。
 * 位置i和 (i二進位制反轉後位置)互換
 * len必須去2的冪
 */
void change(complex y[],int len){
  int i,j,k;
  for(i = 1, j = len/2;i < len-1; i++){
    if(i < j)swap(y[i],y[j]);
    //交換互為小標反轉的元素,i<j保證交換一次
    //i做正常的+1,j左反轉型別的+1,始終保持i和j是反轉的
    k = len/2;
    while( j >= k){
      j -= k;
      k /= 2;
    }
     if(j < k) j += k;
  }
}
/*
 * 做FFT
 * len必須為2^k形式,
 * on==1時是DFT,on==-1時是IDFT
 */
void fft(complex y[],int len,int on){
  change(y,len);
  for(int h = 2; h <= len; h <<= 1){
    complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
    for(int j = 0;j < len;j+=h){
      complex w(1,0);
      for(int k = j;k < j+h/2;k++){
        complex u = y[k];
        complex t = w*y[k+h/2];
        y[k] = u+t;
        y[k+h/2] = u-t;
        w = w*wn;
      }
    }
  }
  if(on == -1)
  for(int i = 0;i < len;i++)
    y[i].r /= len;
}

const int maxn = 200010;
complex x1[maxn],x2[maxn];
char str1[maxn/2],str2[maxn/2];
int sum[maxn];
int main(){
  while(scanf("%s%s",str1,str2)==2){
    int len1 = strlen(str1);
    int len2 = strlen(str2);
    int len = 1;
    while(len < len1*2 || len < len2*2)len<<=1;
    for(int i = 0;i < len1;i++)
      x1[i] = complex(str1[len1-1-i]-'0',0);
    for(int i = len1;i < len;i++)
      x1[i] = complex(0,0);
    for(int i = 0;i < len2;i++)
      x2[i] = complex(str2[len2-1-i]-'0',0);
    for(int i = len2;i < len;i++)
      x2[i] = complex(0,0);
    //求DFT
    fft(x1,len,1);
    fft(x2,len,1);
    for(int i = 0;i < len;i++)
      x1[i] = x1[i]*x2[i];
    fft(x1,len,-1);
    for(int i = 0;i < len;i++)
      sum[i] = (int)(x1[i].r+0.5);
    for(int i = 0;i < len;i++){   // 進位
      sum[i+1]+=sum[i]/10;
      sum[i]%=10;
    }
    len = len1+len2-1;
    while(sum[len] <= 0 && len > 0)len--;// 檢索最高位
    for(int i = len;i >= 0;i--)  // 倒序輸出
      printf("%c",sum[i]+'0');
    printf("\n");
  }
  return 0;
}