[噼昂!]叕是斯特林數
阿新 • • 發佈:2021-10-30
\[\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}
\]
,那麼,除去 \(2\) 的部分的答案實際上就是:
\[Ans=\prod_{p=0}^{\log_2 n}f\bra{\ddiv{n}{2^p}}\times 2^{\ddiv{n}{2^p}} \pmod{2^{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}
\] 很小,可以遞推,特殊處理一下遞推遇到分母為偶數的情況即可。關於程式碼,我寫得不是很好,複雜度有四個 \(\rm log\).
壹、關於題目 ¶
求 \(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^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])\)
注意到後面的東西的 \(j+1\)
叄、參考程式碼 ¶
#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;
}