【BZOJ 2179】 FFT快速傅立叶

最新推荐文章于 2018-11-26 15:01:00 发布

原创最新推荐文章于 2018-11-26 15:01:00 发布 · 280 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#fft

本文介绍了一种利用快速傅立叶变换（FFT）实现高精度乘法的方法，该方法可以在O(nlogn)的时间复杂度内完成两个大整数的乘法运算。文章通过将十进制数转换为多项式系数表示，并直接对这些系数应用FFT，从而有效地实现了这一过程。

思路：

$\text{FFT}$ 的模板题，也是重要的用途之一，那就是 $O(n\log n)$ 的高精度乘法。
可以吧一个十进制数写成一个高精度的形式：

x = \sum i = 0 l e n a i \times 10 i

$x=\sum_{i=0}^{len}a_i\times 10^i$

所以可以直接用系数做 $\text{FFT}$ 。

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const double pi = acos(-1);
const int maxn = 60010*3;
int Bit, rx[maxn], res[maxn];

struct Com{
    double x, y;
    Com(double X = 0, double Y = 0){x = X, y = Y;}
    Com(const Com &t){x = t.x, y = t.y;}
    Com operator + (const Com &t)const{return Com(x+t.x, y+t.y);}
    Com operator - (const Com &t)const{return Com(x-t.x, y-t.y);}
    Com operator * (const Com &t)const{return Com(x*t.x-y*t.y, x*t.y+y*t.x);}
};

void fft(Com *A, int n, int f){
    if(rx[n-1] != n-1) for(int i = 0; i < n; i ++) rx[i] = (rx[i>>1]>>1) | ((i&1) << (Bit-1));
    for(int i = 0; i < n; i ++) if(i < rx[i]) swap(A[i], A[rx[i]]);
    for(int i = 1; i < n; i <<= 1){
        const Com Base(cos(pi/i), f*sin(pi/i));
        for(int j = 0; j < n; j += i << 1){
            Com mi(1, 0);
            for(int k = 0; k < i; k ++, mi = mi * Base){
                Com x = A[j+k], y = A[i+j+k] * mi;
                A[j+k] = x+y, A[i+j+k] = x-y;
            }
        }
    }
}

char c1[maxn], c2[maxn];
Com a[maxn], b[maxn];

int main(){
    int n, m;
    scanf("%d", &n), m = n;
    scanf("%s%s", c1, c2);
    n --, m --;
    for(int i = n; i >= 0; i --){
        a[n-i].x = c1[i] - '0';
        b[n-i].x = c2[i] - '0';
    }
    m += n;
    for(n = 1; n <= m; n <<= 1, Bit ++);
    fft(a, n, 1);
    fft(b, n, 1);
    for(int i = 0; i < n; i ++) a[i] = a[i] * b[i];
    fft(a, n, -1);
    for(int i = 0; i <= m; i ++) res[i] = int(a[i].x/n + 0.5);
    for(int i = 0; i <= m; i ++) res[i+1] += res[i] / 10, res[i] %= 10;
    m ++;
    while(res[m] / 10) res[m+1] += res[m] / 10, res[m] %= 10;
    while(m != 0 && res[m] == 0) m --;
    for(int i = m; i >= 0; i --) printf("%d", res[i]);
    return 0;
}