对于大的
k
,我们可以通过利用两个基本事实来显著减少工作量:
如果p
是一个质数,则n!
的质因数分解中p
的指数由(n - s_p(n)) / (p-1)
给出,其中s_p(n)
是在以p
为底的表示下n
的各位数字之和(所以对于p = 2
,它是popcount)。因此,在choose(n,k)
的质因数分解中p
的指数是(s_p(k) + s_p(n-k) - s_p(n)) / (p-1)
,特别地,当且仅当在以p
为底时加法k + (n-k)
没有进位时,它为零(指数是进位的数量)。
威尔逊定理:如果p
是一个质数,则(p-1)! ≡ (-1) (mod p)
。
n!
的质因数分解中p
的指数通常通过以下方式计算:
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)
是否可以被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
,其中0 ≤ j < q
,以及(q*p + k), 1 ≤ k ≤ r
的乘积。
对于与p
互质的数字,我们只关心它们在模p
下的贡献。每个完整的序列j*p + k, 1 ≤ k < p
在模p
下都等于(p-1)!
,因此它们总共产生了一个贡献(-1)^q
在模p
下。最后(可能)不完整的序列在模p
下等于r!
。
因此,如果我们写成:
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))
任何能够整除
choose(n,k)
的
p
的幂都来自于
choose(a,b)
,在我们的情况下,由于我们已经在之前排除了这些情况,所以不会有这种情况出现。虽然
cop(a,A) / (cop(b,B) * cop(c,C))
不一定是整数(例如考虑
choose(19,9) (mod 5)
),但是当考虑模
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)
现在是实施阶段:
long long choose_mod_one(long long n, long long k, long long p)
{
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);
choose *= choose_mod_one(q_n, q_k, p);
return choose % p;
}
long long choose_mod_two(long long n, long long k, long long p)
{
n %= p;
if (n < k) return 0;
if (k == 0 || k == n) return 1;
if (k > n/2) k = n-k;
long long num = n, den = 1;
for(n = n-1; k > 1; --n, --k)
{
num = (num * n) % p;
den = (den * k) % 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))
,所以对于小的 k
和大的 p
,速度会更快),但对于其他情况则会慢一些。
总体而言,我们需要计算最多 O(log_p k) 个二项式系数模 p
,其中每个二项式系数最多需要 O(p) 次操作,从而得到总复杂度为 O(p*log_p k) 的操作次数。
当 k
明显大于 p
时,这比 O(k)
的解决方案要好得多。对于 k <= p
,它会带来一些额外开销,但仍然可以归结为 O(k)
的解决方案。