卢卡斯定理
前置知识:阶乘取模
引入
本文讨论大组合数取模的求解。组合数,又称二项式系数,指表达式:
规模不大时,组合数可以通过 递推公式 求解,时间复杂度为
基于 Lucas 定理及其推广,本文讨论一种可以在模数不太大 (
Lucas 定理
首先讨论模数为素数
Lucas 定理
对于素数
其中,当
利用生成函数证明
考虑
所以,当
其中,第三行的同余利用了前文说明的结论,即只有
利用这一结论,考察二项式展开:
等式左侧中,项
转而计算等式右侧中项
令两侧系数相等,就得到 Lucas 定理。
利用阶乘取模的结论证明
此处提供一种基于 阶乘取模 相关结论的证明方法,以方便和后文 exLucas 部分的方法建立联系。已知二项式系数
将阶乘
就得到二项式系数的表达式:
幂次
前者是 Legendre 公式的推论,后者是 Wilson 定理的推论。
将递推公式代入二项式系数的表达式并整理,就得到:
现在考察
所以,利用第一式减去后两式,就得到
等式右侧,前两项的和严格小于
- 如果它是
,那么此时也成立 。因此,上式中的第一个因子的指数为 ,该因子就等于一;第二个因子就是 ;第三个因子则由前文的展开式可知,就等于 。此时,Lucas 公式成立; -
如果它是
,那么第一个因子的指数为 ,该因子就等于零,所以二项式系数的余数为零。同时,Lucas 定理所要证明的等式右侧的 也必然是零,因为此时必然有 ;否则,将有 这显然与余数的定义矛盾。
综合两种情形,就得到了所要求证的 Lucas 定理。这一证明说明,在求解素数模下组合数时,利用 Lucas 定理和利用 exLucas 算法得到的结果是等价的。
Lucas 定理指出,模数为素数
示意
其中,C(n, k, p)
用于计算小规模的组合数。
递归至多进行
参考实现
此处给出的参考实现在
参考实现
#include <iostream>
#include <vector>
class BinomModPrime {
int p;
std::vector<int> fa, ifa;
// Calculate binom(n, k) mod p for n, k < p.
int calc(int n, int k) {
if (n < k) return 0;
long long res = fa[n];
res = (res * ifa[k]) % p;
res = (res * ifa[n - k]) % p;
return res;
}
public:
BinomModPrime(int p) : p(p), fa(p), ifa(p) {
// Factorials mod p till p.
fa[0] = 1;
for (int i = 1; i < p; ++i) {
fa[i] = (long long)fa[i - 1] * i % p;
}
// Inverse of factorials mod p till p.
ifa[p - 1] = p - 1; // Wilson's theorem.
for (int i = p - 1; i; --i) {
ifa[i - 1] = (long long)ifa[i] * i % p;
}
}
// Calculate binom(n, k) mod p.
int binomial(long long n, long long k) {
long long res = 1;
while (n || k) {
res = (res * calc(n % p, k % p)) % p;
n /= p;
k /= p;
}
return res;
}
};
int main() {
int t, p;
std::cin >> t >> p;
BinomModPrime bm(p);
for (; t; --t) {
long long n, k;
std::cin >> n >> k;
std::cout << bm.binomial(n, k) << '\n';
}
return 0;
}
该实现的时间复杂度为
exLucas 算法
Lucas 定理中对于模数
素数幂模的情形
首先考虑模数为素数幂
其中,
式子中的
注意,如果幂次
一般模数的情形
对于
然后,分别计算出模
最后,利用 中国剩余定理 求出模
参考实现
最后,给出模板题目 二项式系数 的参考实现。
参考实现
#include <iostream>
#include <vector>
// Extended Euclid.
void ex_gcd(int a, int b, int& x, int& y) {
if (!b) {
x = 1;
y = 0;
} else {
ex_gcd(b, a % b, y, x);
y -= a / b * x;
}
}
// Inverse of a mod m.
int inverse(int a, int m) {
int x, y;
ex_gcd(a, m, x, y);
return (x % m + m) % m;
}
// Coefficient in CRT.
int crt_coeff(int m_i, int m) {
long long mm = m / m_i;
mm *= inverse(mm, m_i);
return mm % m;
}
// Binominal Coefficient Calculator Modulo Prime Power.
class BinomModPrimePower {
int p, a, pa;
std::vector<int> f;
// Obtain multiplicity of p in n!.
long long nu(long long n) {
long long count = 0;
do {
n /= p;
count += n;
} while (n);
return count;
}
// Calculate (n!)_p mod pa.
long long fact_mod(long long n) {
bool neg = p != 2 || pa <= 4;
long long res = 1;
while (n > 1) {
if ((n / pa) & neg) res = pa - res;
res = res * f[n % pa] % pa;
n /= p;
}
return res;
}
public:
BinomModPrimePower(int p, int a, int pa) : p(p), a(a), pa(pa), f(pa) {
// Pretreatment.
f[0] = 1;
for (int i = 1; i < pa; ++i) {
f[i] = i % p ? (long long)f[i - 1] * i % pa : f[i - 1];
}
}
// Calculate Binom(n, k) mod pa.
int binomial(long long n, long long k) {
long long v = nu(n) - nu(n - k) - nu(k);
if (v >= a) return 0;
auto res = fact_mod(n - k) * fact_mod(k) % pa;
res = fact_mod(n) * inverse(res, pa) % pa;
for (; v; --v) res *= p;
return res % pa;
}
};
// Binominal Coefficient Calculator.
class BinomMod {
int m;
std::vector<BinomModPrimePower> bp;
std::vector<long long> crt_m;
public:
BinomMod(int n) : m(n) {
// Factorize.
for (int p = 2; p * p <= n; ++p) {
if (n % p == 0) {
int a = 0, pa = 1;
for (; n % p == 0; n /= p, ++a, pa *= p);
bp.emplace_back(p, a, pa);
crt_m.emplace_back(crt_coeff(pa, m));
}
}
if (n > 1) {
bp.emplace_back(n, 1, n);
crt_m.emplace_back(crt_coeff(n, m));
}
}
// Calculate Binom(n, k) mod m.
int binomial(long long n, long long k) {
long long res = 0;
for (size_t i = 0; i != bp.size(); ++i) {
res = (bp[i].binomial(n, k) * crt_m[i] + res) % m;
}
return res;
}
};
int main() {
int t, m;
std::cin >> t >> m;
BinomMod bm(m);
for (; t; --t) {
long long n, k;
std::cin >> n >> k;
std::cout << bm.binomial(n, k) << '\n';
}
return 0;
}
该算法在预处理时将模数
习题
创建日期: 2019年7月24日