组合数取模方法总结(Lucas定理介绍)

2019-04-14 18:04发布

转自:https://www.cnblogs.com/fzl194/p/9095177.html

1.当n,m都很小的时候可以利用杨辉三角直接求。 
C(n,m)=C(n-1,m)+C(n-1,m-1);   2、n和m较大,但是p为素数的时候 Lucas定理是用来求 c(n,m) mod p,p为素数的值C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p   也就是Lucas(n,m)%p=Lucas(n/p,m/p)*C(n%p,m%p)%p     求上式的时候,Lucas递归出口为m=0时返回1 求C(n%p, m%p)%p的时候,此处写成C(n, m)%p(p是素数,n和m均小于p) C(n, m)%p = n! / (m ! * (n - m )!) % p = n! * mod_inverse[m! * (n - m)!, p] % p 由于p是素数,有费马小定理可知,m! * (n - m)! 关于p的逆元就是m! * (n - m)!的p-2次方。   p较小的时候预处理出1-p内所有阶乘%p的值,然后用快速幂求出逆元,就可以求出解。p较大的时候只能逐项求出分母和分子模上p的值,然后通过快速幂求逆元求解。 P较大,不打表: ll pow(ll a, ll b, ll m) { ll ans = 1; a %= m; while(b) { if(b & 1)ans = (ans % m) * (a % m) % m; b /= 2; a = (a % m) * (a % m) % m; } ans %= m; return ans; } ll inv(ll x, ll p)//x关于p的逆元,p为素数 { return pow(x, p - 2, p); } ll C(ll n, ll m, ll p)//组合数C(n, m) % p { if(m > n)return 0; ll up = 1, down = 1;//分子分母; for(int i = n - m + 1; i <= n; i++)up = up * i % p; for(int i = 1; i <= m; i++)down = down * i % p; return up * inv(down, p) % p; } ll Lucas(ll n, ll m, ll p) { if(m == 0)return 1; return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; }   P较小,打表: const int maxn = 1e5 + 10; ll fac[maxn];//阶乘打表 void init(ll p)//此处的p应该小于1e5,这样Lucas定理才适用 { fac[0] = 1; for(int i = 1; i <= p; i++) fac[i] = fac[i - 1] * i % p; } ll pow(ll a, ll b, ll m) { ll ans = 1; a %= m; while(b) { if(b & 1)ans = (ans % m) * (a % m) % m; b /= 2; a = (a % m) * (a % m) % m; } ans %= m; return ans; } ll inv(ll x, ll p)//x关于p的逆元,p为素数 { return pow(x, p - 2, p); } ll C(ll n, ll m, ll p)//组合数C(n, m) % p { if(m > n)return 0; return fac[n] * inv(fac[m] * fac[n - m], p) % p; } ll Lucas(ll n, ll m, ll p) { if(m == 0)return 1; return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; }   3、n,m较大且p不为素数的时候 扩展Lucas定理:   #include using namespace std; typedef long long ll; const int maxn = 1e6 + 10; const int mod = 1e9 + 7; ll pow(ll a, ll b, ll m) { ll ans = 1; a %= m; while(b) { if(b & 1)ans = (ans % m) * (a % m) % m; b /= 2; a = (a % m) * (a % m) % m; } ans %= m; return ans; } ll extgcd(ll a, ll b, ll& x, ll& y) //求解ax+by=gcd(a, b) //返回值为gcd(a, b) { ll d = a; if(b) { d = extgcd(b, a % b, y, x); y -= (a / b) * x; } else x = 1, y = 0; return d; } ll mod_inverse(ll a, ll m) //求解a关于模上m的逆元 //返回-1表示逆元不存在 { ll x, y; ll d = extgcd(a, m, x, y); return d == 1 ? (m + x % m) % m : -1; } ll Mul(ll n, ll pi, ll pk)//计算n! mod pk的部分值 pk为pi的ki次方 //算出的答案不包括pi的幂的那一部分 { if(!n)return 1; ll ans = 1; if(n / pk) { for(ll i = 2; i <= pk; i++) //求出循环节乘积 if(i % pi)ans = ans * i % pk; ans = pow(ans, n / pk, pk); //循环节次数为n / pk } for(ll i = 2; i <= n % pk; i++) if(i % pi)ans = ans * i % pk; return ans * Mul(n / pi, pi, pk) % pk;//递归求解 } ll C(ll n, ll m, ll p, ll pi, ll pk)//计算组合数C(n, m) mod pk的值 pk为pi的ki次方 { if(m > n)return 0; ll a = Mul(n, pi, pk), b = Mul(m, pi, pk), c = Mul(n - m, pi, pk); ll k = 0, ans;//k为pi的幂值 for(ll i = n; i; i /= pi)k += i / pi; for(ll i = m; i; i /= pi)k -= i / pi; for(ll i = n - m; i; i /= pi)k -= i / pi; ans = a * mod_inverse(b, pk) % pk * mod_inverse(c, pk) % pk * pow(pi, k, pk) % pk;//ans就是n! mod pk的值 ans = ans * (p / pk) % p * mod_inverse(p / pk, pk) % p;//此时用剩余定理合并解 return ans; } ll Lucas(ll n, ll m, ll p) { ll x = p; ll ans = 0; for(ll i = 2; i <= p; i++) { if(x % i == 0) { ll pk = 1; while(x % i == 0)pk *= i, x /= i; ans = (ans + C(n, m, p, i, pk)) % p; } } return ans; } int main() { ll n, m, p; while(cin >> n >> m >> p) { cout<