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;
    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);
    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);
    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;
}


Disclaimer: some contents on this website are collected through internet etc. Please notify if violated the original author's copyright and we will delete it immediately.
免责声明：本站部分文章来源于网络等其它媒体，如果侵犯了原作者的版权，请联系我们，本站将立即删除。