对于大k
,我们可以通过利用两个基本事实来显着减少工作:
如果p
是素数,则p
的素数因式分解中的指数n!
由 给出(n - s_p(n)) / (p-1)
,其中是基本表示s_p(n)
中的数字之和(因此对于,它是 popcount)。因此,在 的素数分解中的指数是,特别是,当且仅当在基数中执行加法时没有进位(指数是进位的数量)时,它才为零。n
p
p = 2
p
choose(n,k)
(s_p(k) + s_p(n-k) - s_p(n)) / (p-1)
k + (n-k)
p
威尔逊定理: p
是素数,当且仅当(p-1)! ≡ (-1) (mod p)
。
p
因式分解中的指数n!
通常由下式计算
long long factorial_exponent(long long n, long long p)
{
long long ex = 0;
do
{
n /= p;
ex += n;
}while(n > 0);
return ex;
}
choose(n,k)
对by的可分性检查p
并不是绝对必要的,但首先进行检查是合理的,因为通常是这样,然后工作量就会减少:
long long choose_mod(long long n, long long k, long long p)
{
// We deal with the trivial cases first
if (k < 0 || n < k) return 0;
if (k == 0 || k == n) return 1;
// Now check whether choose(n,k) is divisible by p
if (factorial_exponent(n) > factorial_exponent(k) + factorial_exponent(n-k)) return 0;
// If it's not divisible, do the generic work
return choose_mod_one(n,k,p);
}
现在让我们仔细看看n!
。我们将数字≤ n
分成 的倍数p
和互质的数字p
。和
n = q*p + r, 0 ≤ r < p
p
贡献的倍数p^q * q!
。这些数字互质以p
贡献(j*p + k), 1 ≤ k < p
for0 ≤ j < q
的乘积和 的乘积(q*p + k), 1 ≤ k ≤ r
。
对于互质的数字,p
我们只会对贡献模感兴趣p
。每个完整的运行j*p + k, 1 ≤ k < p
都与(p-1)!
模数一致p
,因此它们总共产生(-1)^q
模数的贡献p
。最后一次(可能)不完整的运行产生r!
模数p
。
所以如果我们写
n = a*p + A
k = b*p + B
n-k = c*p + C
我们得到
choose(n,k) = p^a * a!/ (p^b * b! * p^c * c!) * cop(a,A) / (cop(b,B) * cop(c,C))
其中是cop(m,r)
与 互质的所有数的乘积。p
≤ m*p + r
有两种可能性,a = b + c
和A = B + C
,或a = b + c + 1
和A = B + C - p
。
在我们的计算中,我们已经预先排除了第二种可能性,但这不是必须的。
在第一种情况下,p
取消的明确权力,我们留下
choose(n,k) = a! / (b! * c!) * cop(a,A) / (cop(b,B) * cop(c,C))
= choose(a,b) * cop(a,A) / (cop(b,B) * cop(c,C))
任何p
除法choose(n,k)
都来自choose(a,b)
- 在我们的例子中,不会有,因为我们之前已经消除了这些情况 - 并且,虽然cop(a,A) / (cop(b,B) * cop(c,C))
不必是整数(考虑 eg choose(19,9) (mod 5)
),当考虑表达式 modulo 时p
,cop(m,r)
减少到(-1)^m * r!
, 所以,因为a = b + c
,(-1)
取消,我们只剩下
choose(n,k) ≡ choose(a,b) * choose(A,B) (mod p)
在第二种情况下,我们发现
choose(n,k) = choose(a,b) * p * cop(a,A)/ (cop(b,B) * cop(c,C))
因为a = b + c + 1
. 最后一位的进位意味着A < B
,所以模p
p * cop(a,A) / (cop(b,B) * cop(c,C)) ≡ 0 = choose(A,B)
(我们可以用模逆的乘法代替除法,或者将其视为有理数的同余,这意味着分子可以被 整除p
)。无论如何,我们再次发现
choose(n,k) ≡ choose(a,b) * choose(A,B) (mod p)
现在我们可以重现该choose(a,b)
部分。
例子:
choose(144,6) (mod 5)
144 = 28 * 5 + 4
6 = 1 * 5 + 1
choose(144,6) ≡ choose(28,1) * choose(4,1) (mod 5)
≡ choose(3,1) * choose(4,1) (mod 5)
≡ 3 * 4 = 12 ≡ 2 (mod 5)
choose(12349,789) ≡ choose(2469,157) * choose(4,4)
≡ choose(493,31) * choose(4,2) * choose(4,4
≡ choose(98,6) * choose(3,1) * choose(4,2) * choose(4,4)
≡ choose(19,1) * choose(3,1) * choose(3,1) * choose(4,2) * choose(4,4)
≡ 4 * 3 * 3 * 1 * 1 = 36 ≡ 1 (mod 5)
现在实现:
// Preconditions: 0 <= k <= n; p > 1 prime
long long choose_mod_one(long long n, long long k, long long p)
{
// For small k, no recursion is necessary
if (k < p) return choose_mod_two(n,k,p);
long long q_n, r_n, q_k, r_k, choose;
q_n = n / p;
r_n = n % p;
q_k = k / p;
r_k = k % p;
choose = choose_mod_two(r_n, r_k, p);
// If the exponent of p in choose(n,k) isn't determined to be 0
// before the calculation gets serious, short-cut here:
/* if (choose == 0) return 0; */
choose *= choose_mod_one(q_n, q_k, p);
return choose % p;
}
// Preconditions: 0 <= k <= min(n,p-1); p > 1 prime
long long choose_mod_two(long long n, long long k, long long p)
{
// reduce n modulo p
n %= p;
// Trivial checks
if (n < k) return 0;
if (k == 0 || k == n) return 1;
// Now 0 < k < n, save a bit of work if k > n/2
if (k > n/2) k = n-k;
// calculate numerator and denominator modulo p
long long num = n, den = 1;
for(n = n-1; k > 1; --n, --k)
{
num = (num * n) % p;
den = (den * k) % p;
}
// Invert denominator modulo p
den = invert_mod(den,p);
return (num * den) % p;
}
要计算模逆,可以使用费马(所谓的小)定理
如果p
是素数a
且不能被 整除p
,则a^(p-1) ≡ 1 (mod p)
.
并将逆计算为a^(p-2) (mod p)
,或使用适用于更广泛参数范围的方法、扩展欧几里得算法或连分数展开式,它为您提供任何一对互质(正)整数的模逆:
long long invert_mod(long long k, long long m)
{
if (m == 0) return (k == 1 || k == -1) ? k : 0;
if (m < 0) m = -m;
k %= m;
if (k < 0) k += m;
int neg = 1;
long long p1 = 1, p2 = 0, k1 = k, m1 = m, q, r, temp;
while(k1 > 0) {
q = m1 / k1;
r = m1 % k1;
temp = q*p1 + p2;
p2 = p1;
p1 = temp;
m1 = k1;
k1 = r;
neg = !neg;
}
return neg ? m - p2 : p2;
}
就像计算一样a^(p-2) (mod p)
,这是一种O(log p)
算法,对于某些输入,它的速度要快得多(实际上是O(min(log k, log p))
,因此对于 smallk
和 large p
,它要快得多),而对于其他输入则要慢得多。
总的来说,这种方式我们最多需要计算 O(log_p k) 个二项式系数 modulo p
,其中每个二项式系数最多需要 O(p) 个操作,从而产生 O(p*log_p k) 个操作的总复杂度。当k
显着大于时p
,这比解要好得多O(k)
。对于k <= p
,它简化为O(k)
具有一些开销的解决方案。