关于阶乘的快速算法的研究

解决问题

求一个的阶乘$n!modP$。

总述

本文提出一种快速阶乘算法，对于模数P没有质数或是合数的要求，假设$P$对于素数$p$的最大分解是$p^k$，该算法能够做到在$O(ma{x}_{p^k|P}{k}^{2}plo{g}_{p}n)$的时间复杂度内将阶乘分解成$r\ast p^t$的形式，其中$r$与$p$互质。
本文在文末指出了优化方向，能够将该算法优化至$O(ma{x}_{p^k|P}\sqrt{p}lo{g}_{2}p+k^{2}lo{g}_{p}n)$的时间复杂度。

依赖算法

1. 扩展欧几里得算法求逆：通过欧几里得算法求得某个数在某个模数下的逆元。
2. 第一类斯特林数求自然数幂和：可以在对模数无任何要求的情况下，在$O(k^2)$的时间复杂度内求得自然数幂的$1$至$k$次幂和。
3. 扩展卢卡斯定理(exLucas)：能够在$O(ma{x}_{p^k|P}{p}^k)$的时间复杂度内求得阶乘。
4. 中国剩余定理(CRT)：通过合并各个因子，得到模P的具体值。

算法综述

首先将模数P质因数分解，考虑$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,p\nmid i}^{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$的模系下为0。
此时有$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)$。
对于函数$g(t,d,p,k)={\prod }_{i=0}^t(ip+d)$，我们将其展开后能够写成${\sum }_{i=0}^{k-1}{a}_i{p}^i$，对于$p^c,c>=k$的项，会被模数$p^k$约去。我们发现能够利用自然数幂和配合容斥原理求得每一项系数，最后求得函数$g$的值。

后续优化方向

1.我们发现函数$f$在代码103行处求了一次阶乘，该求阶乘在仅在$n<p$的情况时求一次。对于小于模数的阶乘的快速求解在洛谷P5282处有$O(\sqrt{p}lo{g}_{2}p)$的做法。
2.函数$f$在代码108行处调用了函数$g$，而每次调用函数$g$的参数$t$只有两种取值，而参数$d$虽然每次不同但是$d$参数仅在函数$g$的96行处被调用。因此函数$g$的其余部分是相同的可以复用。
结合两处优化可将快速阶乘算法优化至$O(ma{x}_{p^k|P}\sqrt{p}lo{g}_{2}p+k^{2}lo{g}_{p}n)$。

例题

51nod1151N!的非0最低18位

代码

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

typedef long long ll;
typedef __int128 Int;

const int N = 1e5 + 100;

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 r = 1;
	while (n > 0) {
		if (n & 1) r = mul(r, x, mod);
		x = mul(x, x, mod); n /= 2;
	}
	return r;
}

namespace Fac {
	Int mod;
	Int add(Int a, Int b) { return a + b >= mod ? a + b - mod : a + b; }
	Int mul(Int a, Int b) { return a * b % mod; }
	Int qpow(Int x, Int n) {
		Int r = 1;
		while (n > 0) {
			if (n & 1) r = mul(r, x);
			x = mul(x, x); n /= 2;
		}
		return r;
	}

	Int s1[110][110];
	void init() {
		s1[0][0] = 1;
		for (int i = 1; i <= 50; i++) {
			for (int j = 1; j <= i; j++) {
				s1[i][j] = add(s1[i - 1][j - 1], mul(i - 1, s1[i - 1][j]));
			}
		}
	}
	Int sum[110];
	void get(Int n, Int k) {
		sum[0] = n % mod; sum[1] = n * (n + 1) / 2 % mod;
		for (int i = 2; i <= k; i++) {
			sum[i] = 1;
			for (int j = 0; j < i + 1; j++) {
				if ((n - j + 1) % (i + 1) == 0) sum[i] = mul(sum[i], (n - j + 1) / (i + 1));
				else sum[i] = mul(sum[i], n - j + 1);
			}
			for (int j = 1; j < i; j++) {
				if ((i - j) % 2 == 0) sum[i] = add(sum[i], mod - mul(s1[i][j], sum[j]));
				else sum[i] = add(sum[i], mul(s1[i][j], sum[j]));
			}
		}
	}
	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 -= a / b * x;
		}
	}
	Int rev(Int a, Int p) {
		Int b, c, d; exgcd(a, p, b, c, d);
		return d == 1 ? (b + p) % p : -1;
	}
	Int g(Int t, Int d, Int p, Int k) {
		if (t <= k) {
			Int res = 1;
			for (int i = 0; i <= t; i++) res = mul(res, add(mul(i, p), d));
			return res;
		}
		else {
			static Int dp[110], cd[110], cp[110], cnt[110];
			cd[k - 1] = qpow(d, t - k + 2);
			for (int i = k - 2; i >= 0; i--) cd[i] = mul(d, cd[i + 1]);
			get(t, k); sum[0] = 1;
			cp[0] = 1;
			for (int i = 1; i < k; i++) cp[i] = mul(cp[i - 1], p);
			cnt[0] = 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(sum[j], mul(dp[i - j], cp[ma - cnt[i - j]]));
					if (f > 0) dp[i] = add(dp[i], r);
					else dp[i] = add(dp[i], mod - r);
				}
				Int x = i; cnt[i] = ma;
				while (x % p == 0) x /= p, cnt[i]++, ma++;
				dp[i] = mul(dp[i], rev(x, mod));
			}
			Int res = 0;
			for (int i = 0; i < k; i++)
				res = add(res, mul(cp[i - cnt[i]], mul(cd[i], dp[i])));
			return res;
		}
	}
	Int f(Int n, Int p, Int k) {
		if (n < p) {
			Int res = 1;
			for (int i = 2; i <= n; i++)  res = mul(res, i);
			return res;
		}
		else {
			Int res = f(n / p, p, k);
			for (int i = 1; i < p; i++) res = mul(res, g(n / p + (n % p >= i ? 0 : -1), i, p, k));
			return res;
		}
	}


	void solve(Int n, Int p, Int k, Int mod, Int &res, Int &cnt) {
		Fac::mod = mod; init();
		res = f(n, p, k);
		cnt = 0;
		while (n > 0) cnt += n / p, n /= p;
	}
}
using Fac::rev;

int main() {
	ll n; scanf("%lld", &n);
	ll n5 = qpow(5, 18, 1e18), n2 = qpow(2, 18, 1e18), nn = n2 * n5;
	Int r2, c2, r5, c5;
	Fac::solve(n, 2, 18, n2, r2, c2); Fac::solve(n, 5, 18, n5, r5, c5);
	r2 = mul(r2, qpow(rev(5, n2), c5, n2), n2);
	r5 = mul(r5, qpow(rev(2, n5), c2, n5), n5);
	Int t5 = mul(mul(rev(n2, n5), n2, nn), r5, nn);
	Int t2 = mul(mul(rev(n5, n2), n5, nn), r2, nn);
	ll r = add(t2, t5, nn);
	r = mul(r, qpow(2, c2 - c5, nn), nn);
	printf("%lld\n", r);
	return 0;
}