ACM学习历程—51NOD1028 大数乘法V2（FFT）

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <set>
#include <map>
#include <queue>
#include <string>
#define LL long long

using namespace std;

//多项式乘法运算
//快速傅里叶变换
//FFT
//用于求两个多项式的卷积,但有精度损失
//时间复杂度nlogn
const int maxN = 400005;
const double PI = acos(-1.0);

struct Complex
{
    double r, i;

    Complex(double rr = 0.0, double ii = 0.0)
    {
        r = rr;
        i = ii;
    }

    Complex operator+(const Complex &x)
    {
        return Complex(r+x.r, i+x.i);
    }

    Complex operator-(const Complex &x)
    {
        return Complex(r-x.r, i-x.i);
    }

    Complex operator*(const Complex &x)
    {
        return Complex(r*x.r-i*x.i, i*x.r+r*x.i);
    }
};

//雷德算法--倒位序
//Rader算法
//进行FFT和IFFT前的反转变换。
//位置i和 （i二进制反转后位置）互换
void Rader(Complex y[], int len)
{
    int j = len>>1;
    for(int i = 1; i < len-1; i++)
    {
        if (i < j) swap(y[i], y[j]);
        int k = len >> 1;
        while (j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if (j < k) j += k;
    }
}

//FFT实现
//len必须为2^k形式，
//on==1时是DFT，on==-1时是IDFT
void FFT(Complex y[], int len, int on)
{
    Rader(y, len);
    for (int h = 2; h <= len; h<<=1)//分治后计算长度为h的DFT
    {
        Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h)); //单位复根e^(2*PI/m)用欧拉公式展开
        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;
}

//求卷积
void Conv(Complex a[], Complex b[], int ans[], int len)
{
    FFT(a, len, 1);
    FFT(b, len, 1);
    for (int i = 0; i < len; i++)
        a[i] = a[i]*b[i];
    FFT(a, len, -1);
    //精度复原
    for(int i = 0; i < len; i++)
        ans[i] = a[i].r+0.5;
}

//进制恢复
//用于大数乘法
void turn(int ans[], int len, int unit)
{
    for(int i = 0; i < len; i++)
    {
        ans[i+1] += ans[i]/unit;
        ans[i] %= unit;
    }
}

char str1[maxN], str2[maxN];
Complex za[maxN],zb[maxN];
int ans[maxN];
int len;

void init(char str1[], char str2[])
{
    int len1 = strlen(str1);
    int len2 = strlen(str2);
    len = 1;
    while (len < 2*len1 || len < 2*len2) len <<= 1;

    int i;
    for (i = 0; i < len1; i++)
    {
        za[i].r = str1[len1-i-1]-‘0‘;
        za[i].i = 0.0;
    }
    while (i < len)
    {
        za[i].r = za[i].i = 0.0;
        i++;
    }
    for (i = 0; i < len2; i++)
    {
        zb[i].r = str2[len2-i-1]-‘0‘;
        zb[i].i = 0.0;
    }
    while (i < len)
    {
        zb[i].r = zb[i].i = 0.0;
        i++;
    }
}

void work()
{
    Conv(za, zb, ans, len);
    turn(ans, len, 10);
    while (ans[len-1] == 0) len--;
    for (int i = len-1; i >= 0; i--)
        printf("%d", ans[i]);
    printf("\n");
}

int main()
{
    //freopen("test.in", "r", stdin);
    while (scanf("%s%s", str1, str2) != EOF)
    {
        init(str1, str2);
        work();
    }
    return 0;
}