51nod1151 N!的非0最低18位

题面

题目链接

解题思路

我们将$10^{18}$拆成$2^{18}\ast 5^{18}$分别求解，再用CRT合并。
考虑$n!\%p^k$如何求解。
我们首先应将n!中所有的p因子剔除。这一步可以递归求解，令f(n,p,k)表示$n!\%p^k$且剔除了p的所有因子。
则有$f(n,p,k)=\{\begin{array}{l}{\prod }_{i=1}^{n}i,n<p\\ f(n/p,p,k)\ast {\prod }_{i=1,i\nmid p}^{n}i,n>=p\end{array}$
考虑${\prod }_{i=1,i\nmid p}^{n}i$如何计算，将模p同余的数一起处理。
$令g(t,d,p,k)={\prod }_{i=0}^t(ip+d)$。可以将该式以p的多项式的形式展开，当次数大于等于k时，在$p^k$的模系下可以被忽略。
此时有$f(n,p,k)=\{\begin{array}{l}{\prod }_{i=1}^{n}i,n<p\\ f(n/p,p,k)\ast {\prod }_{i=1}^{p-1}g(t,i,p,k),n>=p\end{array}$,其中$t=n/p+(n\%p>=i?0:-1)$。
如此就可以在$O(lo{g}_p^{2}n)$的时间内解决。
然而我比较懒，写了$O(plo{g}_p^{2}n)$的，31ms。将就着看看吧，有兴趣的可以自己去优化。
我再提一些细节，因为模系不是质数，自然数幂和可能必须要使用第二类斯特林数来求。
在求某些系数的时候，因为系数的系数里有p，不能被除掉，所以可以把其他的数都乘p来平衡。

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;

typedef long long ll;
typedef __int128 Int;

Int add(Int a, Int b, Int mod) {
	return (a + b >= mod ? a + b - mod : a + b);
}

Int mul(Int a, Int b, Int mod) {
	return a * b % mod;
}

Int qpow(Int x, Int n, Int mod) {
	Int res = 1;
	for (; n; n >>= 1, x = x * x % mod) if (n & 1) res = res * x % mod;
	return res;
}

Int s2[110][110];

void init(Int k, Int mod) {
	s2[0][0] = 1;
	for (int i = 1; i <= k + 10; i++) {
		for (int j = 1; j <= i; j++) {
			s2[i][j] = add(s2[i - 1][j - 1], mul(j, s2[i - 1][j], mod), mod);
		}
	}
}

Int get(Int n, Int k, Int mod) {
	if (n == 0) return 0;
	Int res = 0;
	for (int i = 1; i <= k && i <= n; i++)
	{
		Int sum = s2[k][i], flag = true;
		for (ll j = n - i + 1; j <= n + 1; j++) {
			if (j % (i + 1) == 0 && flag) sum = mul(sum, j / (i + 1), mod), flag = false;
			else sum = mul(sum, j, mod);
		}
		res = add(res, sum, mod);
	}
	return res;
}

void exgcd(Int a, Int b, Int &x, Int &y, Int &d) {
	if (b == 0) { x = 1; y = 0; d = a; }
	else exgcd(b, a % b, y, x, d), y -= x * (a / b);
}

Int rev(Int a, Int p) {
	Int x, y, d;
	exgcd(a, p, x, y, d);
	return d == 1 ? (x + p) % p : -1;
}

Int g(Int t, Int d, Int p, Int k, Int mod) {
	if (t <= k) {
		Int res = 1;
		for (int i = 0; i <= t; i++) res = mul(res, add(mul(i, p, mod), d, mod), mod);
		return res;
	}
	else {
		static Int dp[110], cd[110], cc[110], cp[110], cnt[110];
		cd[k - 1] = qpow(d, t - k + 2, mod);
		for (int i = k - 2; i >= 0; i--) cd[i] = mul(d, cd[i + 1], mod);
		cc[0] = 1;
		for (int i = 1; i < k; i++) cc[i] = get(t, i, mod);
		cp[0] = 1;
		for (int i = 1; i < k; i++) cp[i] = mul(cp[i - 1], p, mod);
		cnt[0] = 0;
		for (int i = 1; i < k; i++) cnt[i] = 0;
		dp[0] = 1;
		Int ma = 0;
		for (int i = 1; i < k; i++) {
			dp[i] = 0;
			for (int j = 1, f = 1; j <= i; j++, f = -f) {
				Int r = mul(cc[j], mul(dp[i - j], cp[ma - cnt[i - j]], mod), mod);
				if (f > 0) dp[i] = add(dp[i], r, mod);
				else dp[i] = add(dp[i], mod - r, mod);
			}
			Int x = i; cnt[i] = ma;
			while (x % p == 0) x /= p, cnt[i]++, ma++;
			dp[i] = mul(dp[i], rev(x, mod), mod);
		}
		Int res = 0;
		for (int i = 0; i < k; i++) 
			res = add(res, mul(cp[i - cnt[i]], mul(cd[i], dp[i], mod), mod), mod);
		return res;
	}
}

Int f(Int n, Int p, Int k, Int mod) {
	if (n < p) {
		Int res = 1;
		for (int i = 2; i <= n; i++) res = mul(res, i, mod);
		return res;
	}
	else {
		Int res = f(n / p, p, k, mod);
		for (int i = 1; i < p; i++)
			res = mul(res, g(n / p + (n % p >= i ? 0 : -1), i, p, k, mod), mod);
		return res;
	}
}

ll n, r5, r2;

int main() {
	//freopen("0.txt", "r", stdin);
	scanf("%lld", &n);
	ll n5 = qpow(5, 18, 1e18), n2 = qpow(2, 18, 1e18), nn = n2 * n5;
	init(18, n5);
	r5 = f(n, 5, 18, n5);
	init(18, n2);
	r2 = f(n, 2, 18, n2);
	ll c5 = 0, c2 = 0, x5 = n, x2 = n;
	while (x5 > 0) c5 += x5 / 5, x5 /= 5;
	while (x2 > 0) c2 += x2 / 2, x2 /= 2;
	r5 = mul(r5, qpow(rev(2, n5), c2, n5), n5);
	r2 = mul(r2, qpow(rev(5, n2), c5, n2), n2);
	ll t5 = Int(n2) * rev(n2, n5) % nn * r5 % nn;
	ll t2 = Int(n5) * rev(n5, n2) % nn * r2 % nn;
	ll r = (t5 + t2) % nn;
	r = r * qpow(2, c2 - c5, nn) % nn;
	printf("%lld\n", r);
	return 0;
}