[噼昂!]叕是斯特林数
作者:互联网
壹、关于题目 ¶
求 \(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;
}
标签:ddiv,斯特林,text,color,噼昂,ull,logn,64 来源: https://www.cnblogs.com/Arextre/p/15485222.html