解决问题
求一个的阶乘 n ! m o d P n! \ mod P n! modP。
总述
本文提出一种快速阶乘算法,对于模数P没有质数或是合数的要求,假设
P
P
P对于素数
p
p
p的最大分解是
p
k
p^k
pk,该算法能够做到在
O
(
m
a
x
p
k
∣
P
k
2
p
l
o
g
p
n
)
O(max_{p^k|P}k^2plog_pn)
O(maxpk∣Pk2plogpn)的时间复杂度内将阶乘分解成
r
∗
p
t
r*p^t
r∗pt的形式,其中
r
r
r与
p
p
p互质。
本文在文末指出了优化方向,能够将该算法优化至
O
(
m
a
x
p
k
∣
P
p
l
o
g
2
p
+
k
2
l
o
g
p
n
)
O(max_{p^k|P}\sqrt plog_2p+k^2log_pn)
O(maxpk∣Pplog2p+k2logpn)的时间复杂度。
依赖算法
- 扩展欧几里得算法求逆:通过欧几里得算法求得某个数在某个模数下的逆元。
- 第一类斯特林数求自然数幂和:可以在对模数无任何要求的情况下,在 O ( k 2 ) O(k^2) O(k2)的时间复杂度内求得自然数幂的 1 1 1至 k k k次幂和。
- 扩展卢卡斯定理(exLucas):能够在 O ( m a x p k ∣ P p k ) O(max_{p^k|P}p^k) O(maxpk∣Ppk)的时间复杂度内求得阶乘。
- 中国剩余定理(CRT):通过合并各个因子,得到模P的具体值。
算法综述
首先将模数P质因数分解,考虑
n
!
%
p
k
n!\%p^k
n!%pk如何求解。
我们首先应将n!中所有的p因子剔除。这一步可以递归求解,令f(n,p,k)表示
n
!
%
p
k
n!\%p^k
n!%pk且剔除了p的所有因子。
则有
f
(
n
,
p
,
k
)
=
{
∏
i
=
1
n
i
,
n
<
p
f
(
n
/
p
,
p
,
k
)
∗
∏
i
=
1
,
p
∤
i
n
i
,
n
>
=
p
f(n,p,k)= \left\{ \begin{array} {lc} \prod_{i=1}^ni, n<p \\ f(n/p,p,k)*\prod_{i=1,p\nmid i}^n i ,n>=p \end{array} \right.
f(n,p,k)={∏i=1ni,n<pf(n/p,p,k)∗∏i=1,p∤ini,n>=p
考虑
∏
i
=
1
,
i
∤
p
n
i
\prod_{i=1,i\nmid p}^n i
∏i=1,i∤pni如何计算,将模p同余的数一起处理。
令
g
(
t
,
d
,
p
,
k
)
=
∏
i
=
0
t
(
i
p
+
d
)
g(t,d,p,k)=\prod_{i=0}^t(ip+d)
g(t,d,p,k)=∏i=0t(ip+d)。可以将该式以p的多项式的形式展开,当次数大于等于k时,在
p
k
p^k
pk的模系下为0。
此时有
f
(
n
,
p
,
k
)
=
{
∏
i
=
1
n
i
,
n
<
p
f
(
n
/
p
,
p
,
k
)
∗
∏
i
=
1
p
−
1
g
(
t
,
i
,
p
,
k
)
,
n
>
=
p
f(n,p,k)= \left\{ \begin{array} {lc} \prod_{i=1}^ni, n<p \\ f(n/p,p,k)*\prod_{i=1}^{p-1}g(t,i,p,k) ,n>=p \end{array} \right.
f(n,p,k)={∏i=1ni,n<pf(n/p,p,k)∗∏i=1p−1g(t,i,p,k),n>=p,其中
t
=
n
/
p
+
(
n
%
p
>
=
i
?
0
:
−
1
)
t=n/p+(n\%p>=i ? 0 : -1)
t=n/p+(n%p>=i?0:−1)。
对于函数
g
(
t
,
d
,
p
,
k
)
=
∏
i
=
0
t
(
i
p
+
d
)
g(t,d,p,k)=\prod_{i=0}^t(ip+d)
g(t,d,p,k)=∏i=0t(ip+d),我们将其展开后能够写成
∑
i
=
0
k
−
1
a
i
p
i
\sum_{i=0}^{k-1}a_ip^i
∑i=0k−1aipi,对于
p
c
,
c
>
=
k
p^c,c>=k
pc,c>=k的项,会被模数
p
k
p^k
pk约去。我们发现能够利用自然数幂和配合容斥原理求得每一项系数,最后求得函数
g
g
g的值。
后续优化方向
1.我们发现函数
f
f
f在代码103行处求了一次阶乘,该求阶乘在仅在
n
<
p
n<p
n<p的情况时求一次。对于小于模数的阶乘的快速求解在洛谷P5282处有
O
(
p
l
o
g
2
p
)
O(\sqrt plog_2p)
O(plog2p)的做法。
2.函数
f
f
f在代码108行处调用了函数
g
g
g,而每次调用函数
g
g
g的参数
t
t
t只有两种取值,而参数
d
d
d虽然每次不同但是
d
d
d参数仅在函数
g
g
g的96行处被调用。因此函数
g
g
g的其余部分是相同的可以复用。
结合两处优化可将快速阶乘算法优化至
O
(
m
a
x
p
k
∣
P
p
l
o
g
2
p
+
k
2
l
o
g
p
n
)
O(max_{p^k|P}\sqrt plog_2p+k^2log_pn)
O(maxpk∣Pplog2p+k2logpn)。
例题
代码
#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;
}