51Nod-1059-N的阶乘 V2

来源：互联网发布：比特币交易软件编辑：程序博客网时间：2024/05/20 17:28
ACM模版
描述

题解

正如评论区那个初一的大佬所说的那样……
初一就玩转 FFT 感觉有些丧心病狂啊……内心溅射一万点伤害！我现在连模版都用不好……
代码

#include <iostream>#include <stdio.h>#include <cmath>#include <algorithm>using namespace std;typedef unsigned long long ull;struct complex{    double real, imag;    complex(double real1 = 0.0, double imag1 = 0.0)    {        real = real1;        imag = imag1;    }    void operator == (const complex &a)    {        real = a.real;        imag = a.imag;    }};complex operator * (const complex &n1, const complex &n2){    return complex(n1.real * n2.real - n1.imag * n2.imag, n1.real * n2.imag + n1.imag * n2.real);}complex operator + (const complex &n1, const complex &n2){    return complex(n1.real + n2.real, n1.imag + n2.imag);}complex operator - (const complex &n1, const complex &n2){    return complex(n1.real - n2.real, n1.imag - n2.imag);}const static int WORDSIZE = 500000;const static int SIZEN = 262144;const static int MOD = 10000;const static int N = 100000;int cnta1, cntb1;int a1[SIZEN];int b1[SIZEN];int rev[SIZEN];char word[WORDSIZE];complex e1[SIZEN >> 1];complex e2[SIZEN >> 1];complex A[SIZEN];complex w1[SIZEN];complex w2[SIZEN];struct BigInteger{    int l, r;};int cntp1, cntq1;BigInteger p1[N];BigInteger q1[N];void Reverse(int N, int bitlength){    for (int i = 0; i < N; i++)    {        int temp = i, len = bitlength;        while (temp)        {            rev[i] |= (temp & 1) << (len--);            temp >>= 1;        }    }}void DFT(complex *A, complex *w1, complex *e, int N){    int k = SIZEN / N;    int bit = 0;    while (k)    {        bit++;        k >>= 1;    }    for (int i = 0; i < N; i++)    {        w1[i] = A[rev[i] >> bit];    }    for (int h = 2; h <= N; h <<= 1)    {        int h1 = h, bit1 = -1;        while (h1)        {            bit1++;            h1 >>= 1;        }        for (int i = 0; i < N; i += h)        {            int k = h >> 1, q = SIZEN >> bit1;            for (int j = i; j < i + k; j++)            {                complex u = w1[j];                complex t = e[q*(j - i)] * w1[j + k];                w1[j] = u + t;                w1[j + k] = u - t;            }        }    }}void Multiplication(int l1, int r1, int l2, int r2, int *a, int *b, int &cntb){    int length1 = r1 - l1, length2 = r2 - l2;    if (length1 <= 40 && length2 <= 40)    {        int tempcntb = cntb + length1 + length2 + 1;        for (int i = cntb; i <= tempcntb; i++)        {            b[i] = 0;        }        tempcntb = cntb;        for (int i = l1; i < r1; i++)        {            int number = a[i];            for (int j = l2; j < r2; j++)            {                int pos = cntb + i + j - l1 - l2;                b[pos] += a[j] * number;                if (b[pos] >= MOD)                {                    b[pos + 1] += b[pos] / MOD;                    b[pos] %= MOD;                }            }        }        int length = length1 + length2 + 1 + tempcntb;        for (int k = tempcntb; k <= length; k++)        {            if (b[k] != 0)            {                cntb = k + 1;            }        }        return ;    }    int tempLength = max(length1, length2);    int N = 1, bitlength = 0;    while (N < tempLength)    {        N <<= 1;        bitlength++;    }    N <<= 1;    for (int i = 0; i < length1; i++)    {        A[i] = complex(a[i + l1], 0);    }    for (int i = length1; i < N; i++)    {        A[i] = complex(0, 0);    }    DFT(A, w1, e1, N);    for (int i = 0; i < length2; i++)    {        A[i] = complex(a[i + l2], 0);    }    for (int i = length2; i < N; i++)    {        A[i] = complex(0, 0);    }    DFT(A, w2, e1, N);    for (int i = 0; i < N; i++)    {        A[i] = w2[i] * w1[i];    }    DFT(A, w1, e2, N);    ull temp = 0, temp1 = 0;    int anslength = 0, lastcntb = cntb;    double ND = N;    for (int i = 0; i < N; i++)    {        temp += w1[i].real / ND + 0.1;        temp1 = temp / MOD; temp %= MOD;        b[cntb++] = (int)temp;        if (temp != 0)        {            anslength = i + 1;        }        temp = temp1;    }    cntb = lastcntb + anslength;}void Print(int len, int *a){    char *out = word + WORDSIZE;    *(--out) = '\n';    int linesize = 1000;    int cnt = len * 4;    if (a[len - 1] <= 9)    {        cnt -= 3;    }    else if (a[len - 1] <= 99)    {        cnt -= 2;    }    else if (a[len - 1] <= 999)    {        cnt--;    }    cnt %= linesize;    cnt = linesize - cnt;    for (int i = 0; i < len; i++)    {        int number = a[i];        *(--out) = number % 10 + '0';        number /= 10;        if (++cnt == linesize)        {            *(--out) = '\n'; cnt = 0;        }        *(--out) = number % 10 + '0';        number /= 10;        if (++cnt == linesize)        {            *(--out) = '\n'; cnt = 0;        }        *(--out) = number % 10 + '0';        number /= 10;        if (++cnt == linesize)        {            *(--out) = '\n'; cnt = 0;        }        *(--out) = number % 10 + '0';        number /= 10;        if (++cnt == linesize)        {            *(--out) = '\n'; cnt = 0;        }    }    while (*out == '0' || *out == '\n')    {        out++;    }    fwrite(out, 1, word + WORDSIZE - out, stdout);}void factorial(){    BigInteger *p = p1, *q = q1;    int cntp = cntp1, cntq = cntq1;    int cnta = cnta1, cntb = cntb1;    int *a = a1, *b = b1;    while (1)    {        if (cntp == 1)        {            Print(p[0].r, a);            break;        }        for (int i = 0; i < (cntp >> 1); i++)        {            q[cntq].l = cntb;            Multiplication(p[i].l, p[i].r, p[cntp - i - 1].l, p[cntp - i - 1].r, a, b, cntb);            q[cntq++].r = cntb;        }        if ((cntp & 1) == 1)        {            q[cntq].l = cntb;            int l = p[(cntp >> 1)].l;            int r = p[(cntp >> 1)].r;            for (int i = l; i < r; i++)            {                b[cntb++] = a[i];            }            q[cntq++].r = cntb;        }        swap(p, q);        swap(a, b);        swap(cntp, cntq);        swap(cnta, cntb);        cntq = 0;        cntb = 0;    }}void init(int n){    int l = 0;    int n1 = min(n, 99);    int number = 1;    for (int i = 1; i <= n1; i++)    {        if (number * i >= MOD)        {            a1[cnta1++] = number;            p1[cntp1].l = l;            p1[cntp1++].r = cnta1;            l = cnta1;            number = i;        }        else number *= i;    }    if (number > 1)    {        a1[cnta1++] = number;        p1[cntp1].l = l;        p1[cntp1++].r = cnta1;        l = cnta1;    }    for (int i = n1 + 1; i <= n; i++)    {        if (i >= MOD)        {            a1[cnta1++] = i % MOD;            a1[cnta1++] = i / MOD;        }        else a1[cnta1++] = i;        p1[cntp1].l = l;        p1[cntp1++].r = cnta1;        l = cnta1;    }    int bitlength = 0, N = SIZEN >> 1;    while (N)    {        bitlength++;        N >>= 1;    }    Reverse(SIZEN, bitlength);    double pin2 = 2.0 * acos(-1.0) / SIZEN;    int N1 = SIZEN >> 2, N2 = SIZEN >> 1;    for (int i = 0; i <= N1; i++)    {        e1[i] = complex(cos(pin2 * i), sin(pin2 * i));    }    for (int i = N1 + 1; i < N2; i++)    {        e1[i] = complex(-e1[N2 - i].real, e1[N2 - i].imag);    }    for (int i = 0; i < N2; i++)    {        e2[i] = complex(e1[i].real, -e1[i].imag);    }}int n;int main(){    scanf("%d", &n);    init(n);    factorial();    return 0;}
阅读全文
0 0