FFT实现的高精度乘法

这里不赘述离散傅里叶变换DFT，FFT的一种解释见

https://zhuanlan.zhihu.com/p/567738394.

这篇笔记很好地解释了FFT用于高精度乘法的原理，首先阐述FFT如何处理多项式，再把大整数看成多项式。

#include <cmath>
#include <iostream>
using namespace std;

#define rep(i, a, b) for (int i = (a); i <= (b); i++)
#define lop(i, a, b) for (int i = (a); i < (b); i++)
#define dwn(i, a, b) for (int i = (a); i >= (b); i--)
#define ceil(a, b) (a + (b - 1)) / b
#define db double

constexpr int N = 4e6 + 10, M = 2e6 + 10;
const double PI = acos(-1), eps = 1e-8;

int T, n, m;

struct Complex
{
    double x, y;
    Complex operator+(const Complex &t) const
    {
        return {x + t.x, y + t.y};
    }
    Complex operator-(const Complex &t) const
    {
        return {x - t.x, y - t.y};
    }
    Complex operator*(const Complex &t) const
    {
        return {x * t.x - y * t.y, x * t.y + y * t.x};
    }

} a[N], b[N];

int rev[N], bit, tot, res[N];
void fft(Complex a[], int inv)
{
    for (int i = 0; i < tot; i++)
    {
        if (i < rev[i])
            swap(a[i], a[rev[i]]); //只需要交换一次就行了，交换两次等于没有换
    }
    for (int mid = 1; mid < tot; mid <<= 1)
    {
        auto w1 = Complex({cos(PI / mid), inv * sin(PI / mid)});
        for (int i = 0; i < tot; i += mid * 2)
        {
            auto wk = Complex({1, 0});     //初始为w(0,mid),定义为w(k,mid)
            for (int j = 0; j < mid; j++, wk = wk * w1) //单位根递推式
            {
                auto x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

void workFFT(int n, int m)
{// a[0, n], b[0, m]
    while ((1 << bit) < n + m + 1)
        bit++;
    tot = 1 << bit;
    for (int i = 0; i < tot; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    //递推(bit<<1)在bit之前，就已经被算出rev,最后一位是否为1
    fft(a, 1), fft(b, 1);
    for (int i = 0; i < tot; i++)
        a[i] = a[i] * b[i]; //点表示法直接运算
    fft(a, -1);//逆变换，点表示法转换为多项式表示法
    for (int i = 0; i <= n + m; i++)
        res[i]  = (int)(a[i].x / tot + 0.5); //向上去整
}

// 一个大数看成多项式乘法，第i位 ai * 10 ^ i
int main()
{
    string sa, sb;
    cin >> sa >> sb;
    n = sa.length() - 1;
    m = sb.length() - 1;
    rep(i, 0, n)
        a[i].x = sa[n - i] - '0';
    rep(i, 0, m)
        b[i].x = sb[m - i] - '0';     
    workFFT(n, m);
    int top = n + m;
    int i = 0;
    while(i <= top)
    {
        if(res[i] >= 10)
        {
            res[i + 1] += res[i] / 10;//进位
            top += i == top;
            res[i] %= 10;
        }
        i++;
    }
    while(!res[top] && top > 0) //去除掉最高位的0，同时确保答案=0的时候不会被抹除掉
        top --;
    dwn(i, top, 0)
        cout << res[i];   
}

考虑workFFT()中语句res[i] = (int)(a[i].x / tot + 0.5);事实上这里除以tot这一步在一般教科书中是定义在逆变换中的。

当认为n有界时，fft()时间复杂度是O(tot * (log tot))，所以总的时间复杂度是O( (n+m) * log(n+m) ).

但是n太大时，关于n的语句都要改用高精度运算(包括代码中的mid, tot都会成为大整数)。例如

for( Bigint i = 1; i<=n; ++i)时间复杂度为O(n log n). for (Bigint mid = 1; mid < tot; mid <<= 1)运算可以取对数转化为for (Bigint logmid = 1; logmid < logtot; ++ logmid )，所以时间复杂度为O(log (tot) * log(log(tot)) ). 于是fft()时间复杂度为O(tot * log (tot) * log(log(tot)) ).