大整数乘法

最新推荐文章于 2025-05-27 09:05:14 发布

gpwoaibiancheng

最新推荐文章于 2025-05-27 09:05:14 发布

阅读量429

点赞数

CC 4.0 BY-SA版权

文章标签： fft 编程

本文链接：https://blog.csdn.net/pandengkeji/article/details/6488081

大整数乘法（FFT版）

10^10000 * 10^10000
long double 必须不少于80位

#include <cmath>
#include <complex>
#include <cstdio>
#include <cstdlib>
using namespace std;

const long double PI = 3.1415926535897932384626433832795L;

int BitRev(int x, int n)
{
    int res = 0;
    for (; n != 1; n /= 2)
    {
        res = res*2+x%2;
        x /= 2;
    }
    return res;
}

void FFT(complex<long double> y[], complex<long double> x[], int n)
{
    for (int i = 0; i < n; ++i)
        y = x[BitRev(i, n)];
    for (int k = 2; k <= n; k *= 2)
    {
        const complex<long double> omiga_unit(cosl(2*PI/k), sinl(2*PI/k));
        for (int i = 0; i < n; i += k)
        {
            complex<long double> omiga(1, 0);
            for (int j = 0; j < k/2; ++j)
            {
                complex<long double> t = omiga*y[i+j+k/2];
                y[i+j+k/2] = y[i+j]-t;
                y[i+j] += t;
                omiga *= omiga_unit;
            }
        }
    }
}

void IFFT(complex<long double> y[], complex<long double> x[], int n)
{
    for (int i = 0; i < n; ++i)
        y = x[BitRev(i, n)];
    for (int k = 2; k <= n; k *= 2)
    {
        const complex<long double> omiga_unit(cosl(-2*PI/k), sinl(-2*PI/k));
        for (int i = 0; i < n; i += k)
        {
            complex<long double> omiga(1, 0);
            for (int j = 0; j < k/2; ++j)
            {
                complex<long double> t = omiga*y[i+j+k/2];
                y[i+j+k/2] = y[i+j]-t;
                y[i+j] += t;
                omiga *= omiga_unit;
            }
        }
    }
    for (int i = 0; i < n; ++i)
        y /= n;
}

int main()
{
    static char str[10001];
    static complex<long double> a[8192], b[8192], f1[8192], f2[8192];
    int a_size = 0, b_size = 0, s_size;

/*a_size为第一数组长度，b_size为第二个数组长度，s_size为总的长度*/

//获取被乘数，存入数组A中，每3位数一个单元
gets(str);

    for (int size = (int)strlen(str); size > 3; str[size -= 3] = '/0')
        a[a_size++] = atof(str+size-3);
    a[a_size++] = atof(str);

//获取乘数，存入数组B中，每3位数一个单元
    gets(str);
    for (int size = (int)strlen(str); size > 3; str[size -= 3] = '/0')
        b[b_size++] = atof(str+size-3);
    b[b_size++] = atof(str);

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

// atof(str)意思是将字符串编程浮点数 //

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    s_size = a_size+b_size;
    for (a_size = s_size; a_size != (a_size&-a_size); a_size += a_size&-a_size);
    FFT(f1, a, a_size);
    FFT(f2, b, a_size);
    for (int i = 0; i < a_size; ++i)
        f1 *= f2;
    IFFT(a, f1, a_size);
    for (int i = 0; i < s_size-1; ++i)
    {
        a[i+1] += a.real()/1000.;
        a = fmodl(a.real(), 1000.L);
    }
    if (!(int)a[s_size-1].real())
        --s_size;
    printf("%d", (int)a[s_size-1].real());
    for (int i = s_size-2; i >= 0; --i)
        printf("%03d", (int)a.real());
    putchar('/n');
    return 0;
}