1. 程式人生 > 其它 >[噼昂!]叕是斯特林數

[噼昂!]叕是斯特林數

\[\color{red}{\text{校長者,真神人也,左馬桶,右永神,會執利筆破邪炁,何人當之?}} \\ \begin{array}{|} \hline \color{pink}{\text{The principal is really a god}} \\ \color{pink}{\text{with a closestool on the left and Yongshen on the right}} \\ \color{pink}{\text{holding a sharp pen to pierce the truth}} \\ \color{pink}{\text{Who can resist him? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{green}{\text{校長は本當に神であり、左側にトイレ、右側にヨンシェンがあり}} \\ \color{green}{\text{鋭いペンを持って真実を突き刺している。誰が彼に抵抗できるだろうか? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{lightblue}{\text{Le principal est vraiment un dieu}} \\ \color{lightblue}{\text{avec des toilettes à gauche et Yongshen à droite}} \\ \color{lightblue}{\text{tenant un stylo pointu pour percer la vérité}} \\ \color{lightblue}{\text{Qui peut lui résister ? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{purple}{\text{Der Direktor ist wirklich ein Gott}} \\ \color{purple}{\text{mit einer Toilette links und Yongshen rechts}} \\ \color{purple}{\text{der einen spitzen Stift hält}} \\ \color{purple}{\text{um die Wahrheit zu durchdringen.}} \\ \color{purple}{\text{Wer kann ihm widerstehen? }} \\ \hline \end{array} \\ \begin{array}{|} \hline \color{cyan}{\text{Principalis deus est, Yongshen a dextris cum latrina}} \\ \color{cyan}{\text{acuto stylo ad perforandum veritatem: quis resistet ei? }} \\ \hline \end{array} \\ \color{red}{\text{對曰:“無人,狗欲當之,還請賜教!”}} \\ \newcommand\bra[1]{\left({#1}\right)} \newcommand\Bra[1]{\left\{{#1}\right\}} \newcommand\dx[0]{\text{dx}} \newcommand\string[2]{\genfrac{\{}{\}}{0pt}{}{#1}{#2}} \newcommand\down[2]{{#1}^{\underline{#2}}} \newcommand\ddiv[2]{\left\lfloor\frac{#1}{#2}\right\rfloor} \newcommand\udiv[2]{\left\lceil\frac{#1}{#2}\right\rceil} \newcommand\lcm[0]{\operatorname{lcm}} \newcommand\set[1]{\left\{{#1}\right\}} \newcommand\ceil[1]{\left\lceil{#1}\right\rceil} \newcommand\floor[1]{\left\lfloor{#1}\right\rfloor} \]

壹、關於題目 ¶

  求 \(n!\)\(16\) 進位制下的的最後 \(16\) 位(去除末尾 \(0\) 與前導 \(0\)).

貳、關於題解 ¶

  如果我們能夠正確處理末尾 \(0\),實際上我們想要求的就是 \(n!\pmod{16^{16}}=n!\pmod{2^{64}}\).

  由於最後我們要去掉末尾 \(0\),於是不難想到把所有數字的 \(2\) 全部拿出來,如果 \(n!\)\(2^k\) 這個因子,事實上我們只需要 \(2^{k\pmod 4}\),因此我們可以把這個部分單獨拿出來做,對於剩餘部分,實際上就是多個奇數相乘,若記 \(f(x)=\prod_{i=0}^{\ddiv{x-1}{2}}(2i+1)\)

,那麼,除去 \(2\) 的部分的答案實際上就是:

\[Ans=\prod_{p=0}^{\log_2 n}f\bra{\ddiv{n}{2^p}}\times 2^{\ddiv{n}{2^p}} \pmod{2^{64}} \]

  後面 \((2^p)^{\ddiv{n-1}{2^p}}\) 就是關於 \(2\) 的部分,我們可以單獨處理,關於 \(f\) 的部分才是問題所在的地方,我們對這個部分再進行一些分析:

\[\begin{aligned}{} f(n)&=\prod_{i=0}^{\ddiv{n-1}{2}}(2i+1)\pmod {2^{64}} \\ \textbf{If}\qquad F(n,x)&\overset\Delta=\prod_{i=0}^{\ddiv{n-1}{2}}(ix+1)\pmod{2^{64}} \\ \textbf{Then}\qquad f(n)&=\sum_{i=0}^{63}2^i[x^i]F(n,x)\pmod{2^{64}} \\ \textbf{If}\qquad g(n,i)&\overset\Delta=[x^i]F(n,x) \\ \Rightarrow g(n,i)&=\frac{1}{i}\sum_{j=0}^{i-1}(-1)^jg(n,i-1-j)\bra{\sum_{k=0}^{\ddiv{n-1}{2}}k^{j+1}}\pmod{2^{64}} \\ \textbf{If}\qquad h(n,p)&\overset\Delta=\sum_{k=0}^{\ddiv{n-1}{2}}k^p\pmod{2^{64}} \\ \textbf{Then}\qquad g(n,i)&=\frac{1}{i}\sum_{j=0}^{i-1}(-1)^jg(n,i-1-j)h(n,j+1)\pmod{2^{64}} \end{aligned} \]

  現在看看我們需要些什麼:插出 \(f(n,k)(k\in [1,64])\)

,求得 \(g(n,i)(i\in [0,63])\),然後把他們加起來,另外還有 \(2^{p\pmod{4}}\) 的部分,然後就做完了。對於 \(h(n,k)\) 的部分,本來我是想插值插出來的,然後發現它好像是對 \(2^{64}\) 取模,求任意數的逆元挺麻煩的(因為是 \(2^{64}\) 作模數,似乎不大可做),於是我把他下降冪展開了:

\[\begin{aligned} h(n,p)&=\sum_{i=0}^{\ddiv{n-1}{2}}i^p \\ &=\sum_{i=0}^{\ddiv{n-1}{2}}\sum_{j=0}^{\min(i,j)}\down{i}{j}{p\brace j} \\ &=\sum_{j=0}^p{p\brace j}j!\sum_{i=0}^{\ddiv{n-1}{2}}{i\choose j} \\ &=\sum_{j=0}^p{p\brace j}j!{\ddiv{n-1}{2}+1\choose j+1} \end{aligned} \]

  注意到後面的東西的 \(j+1\)

很小,可以遞推,特殊處理一下遞推遇到分母為偶數的情況即可。關於程式碼,我寫得不是很好,複雜度有四個 \(\rm log\).

叄、參考程式碼 ¶

#include <bits/stdc++.h>
using namespace std;

#define USING_FREAD

namespace Elaina {

#define rep(i, l, r) for (int i = l, i##_end_ = r; i <= i##_end_; ++i)
#define drep(i, l, r) for (int i = l, i##_end_ = r; i >= i##_end_; --i)
#define fi first
#define se second
#define Endl putchar('\n')
#define bitcnt(s) (__builtin_popcount(s))
#define bitcntll(s) (__builtin_popcountll(s))

    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int, int> pii;

    template<class T> inline T fab(T x) { return x < 0? -x: x; }
    template<class T> inline void getmax(T& x, const T& rhs) { x = max(x, rhs); }
    template<class T> inline void getmin(T& x, const T& rhs) { x = min(x, rhs); }

#ifdef USING_FREAD
# define CHARRECEI qkgetc()
    inline char qkgetc() {
# define BUFFERSIZE 1 << 18
        static char BUF[BUFFERSIZE], *p1 = BUF, *p2 = BUF;
        return (p1 == p2 && (p2 = (p1 = BUF) + fread(BUF, 1, BUFFERSIZE, stdin), p1 == p2))? EOF : *p1++;
# undef BUFFERSIZE
    }
#else
# define CHARRECEI getchar()
#endif

    template<class T> inline T readret(T x) {
        x = 0; char c; bool f = false;
        while (!isdigit(c = CHARRECEI)) if (c == '-') f = true;
        for (x = (c ^ 48); isdigit(c = CHARRECEI); x = (x << 1) + (x << 3) + (c ^ 48));
        return f? -x: x;
    }
    template<class T> inline void readin(T& x) {
        x = 0; char c; bool f = false;
        while (!isdigit(c = CHARRECEI)) if (c == '-') f = true;
        for (x = (c ^ 48); isdigit(c = CHARRECEI); x = (x << 1) + (x << 3) + (c ^ 48));
        if (f) x = -x;
    }
    template<class T, class... Args> inline void readin(T& x, Args&... args) {
        readin(x), readin(args...);
    }
    template<class T> inline void writln(T x, char c = '\n') {
        static int writ_stk[55], writ_ed;
        if (x < 0) putchar('-'), x = -x;
        do writ_stk[++writ_ed] = x % 10, x /= 10; while (x);
        while (writ_ed) putchar(writ_stk[writ_ed--] ^ 48);
        putchar(c);
    }

} using namespace Elaina;

const int logn = 64;
const ull phiForMod = 1ull << 63;
inline ull qkpow(ull a, ull n) {
    ull ret = 1;
    for (; n; n >>= 1, a *= a)
        if (n & 1) ret *= a;
    return ret;
}
/** @b inv[] only for odd numbers */
ull S[logn + 5][logn + 5], fac[logn + 5], inv[logn + 5];
inline void prelude() {
    S[0][0] = 1;
    rep(i, 1, logn) rep(j, 1, i)
        S[i][j] = S[i - 1][j - 1] + 1ull * j * S[i - 1][j];
    fac[0] = 1; rep(i, 1, logn) fac[i] = fac[i - 1] * i;
    inv[1] = 1;
    for (int i = 3; i <= logn + 1; i += 2)
        inv[i] = qkpow(i, phiForMod - 1);
}

/** @warning @p k should be little enough */
inline ull C(ull n, int k) {
    if(n < (ull)k) return 0;
    // printf("Calc :> n == %llu, k == %d\n", n, k);
    ull ret = 1;
    for (int i = 1; i <= k; ++i) {
        ret *= (n - i + 1);
        if (i & 1) ret *= inv[i];
        else {
            int cur = i;
            do {
                assert(!(ret & 1));
                cur >>= 1, ret >>= 1;
            } while (!(cur & 1));
            ret *= inv[cur];
        }
    }
    // printf("C(%llu, %d) == %llu\n", n, k, ret);
    return ret;
}

ull g[logn + 5], h[logn + 5];
#define sign(i) (((i) & 1)? ((ull)(-1)): 1ull)
inline void solveFunc(ull n) {
    // printf("<-----------  Now n == %llu ----------->\n", n);
    /** solve @b h[] first */
    ull fake_n = (n - 1) / 2;
    rep(p, 0, logn) {
        h[p] = 0;
        rep(j, 0, p) {
            h[p] += S[p][j] * fac[j] * C(fake_n + 1, j + 1);
            // printf("h[%d] += S[%d, %d] * fac[%d] * C(%llu, %d)\n", p, p, j, j, fake_n + 1, j + 1);
        }
    }

    // rep(p, 0, logn) printf("h[%d] == %llu\n", p, h[p]);
    
    /** Then solve @b g[] */
    g[0] = 1;
    rep(i, 1, logn - 1) {
        g[i] = 0;
        for (int j = 0; j < i; ++j)
            g[i] += sign(j) * g[i - 1 - j] * h[j + 1];
        if(i & 1) g[i] *= inv[i];
        else {
            int cur = i;
            do {
                assert(!(g[i] & 1));
                g[i] >>= 1, cur >>= 1;
            } while (!(cur & 1));
            g[i] *= inv[cur];
        }
        // printf("g[%d] == %llu\n", i, g[i]);
    }
}

inline ull solve(ull n) {
    if(n == 0) return 1;
    solveFunc(n);
    ull f = 0;
    for (int j = 0; j < logn; ++j)
        f += g[j] << j;
    // printf("When n == %llu, f == %llu\n", n, f);
    return f;
}

inline void trans(ull f) {
    static int num[20], len = 0;
    while(f) num[++len] = f % 16, f >>= 4;
    while(len) {
        if(num[len] < 10) putchar(num[len] ^ 48);
        else putchar('A' + num[len] - 10);
        --len;
    }
    Endl;
}

ull n;
inline void calcAns() {
    ull ans = 1;
    rep(p, 0, logn - 1) {
        ans *= solve(n >> p);
        // printf("Now ans == %llu\n", ans);
    }
    int cnt = 0;
    rep(p, 1, logn - 1) (cnt += (n >> p) % 4) %= 4;
    // printf("cnt == %d, ans == %llu\n", cnt, ans);
    ans = ans * (1ull << cnt);
    trans(ans);
}

inline void work() {
    readin(n);
    calcAns();
}

signed main() {
    // freopen("multiplication.in", "r", stdin);
    // freopen("multiplication.out", "w", stdout);
    prelude();
    rep(_, 1, readret(1)) work();
    return 0;
}