HDU 5446 Unknown Treasure (Lucas + 中国剩余定理)

题目大意

求 C(m,n) % (p1 * p2 * … *pk)

其中，( 1 ≤ m ≤ n ≤ 10^18 ，1 ≤ k ≤ 10 )，M = p1 ⋅ p2 ⋅⋅⋅ pk ≤ 10^18 and pi ≤ 10^5，pi 是素数。

分析

由于是对大组合数取模，且n 和 m 均大于10^5，所以很容易想到Lucas定理，但Lucas定理的模数p要求是一个素数， 且只能处理p < 10^5.

因此，先用Lucas定理对每一个pi取模得到数组ai，然后在用中国剩余定理合并答案

Lucas定理

，并且p是素数，若有



则有


C(m,n) % p = n! / (m! * (n-m)!) % p ，令 X 为 m! * (n - m)! 的逆元，则

C(m,n) = n! * X (mod p)

因此，对每一项C(mi , ni)采用逆元计算即可。

中国剩余定理

对于方程组 x = ai (mod mi)，若所有模mi两两互素，求x。

令M 为所有mi的乘积，wi = M / mi ，则gcd(wi , mi) = 1。

用扩展欧几里德算法可以找到pi和qi使得wi*pi + mi*qi = 1，然后令 ei = wi*pi，pi为wi对于模mi的逆。

原方程等价于

x = e1*a1 + e2*a2 + ... + en*an (mod M)

因为wi是mj的倍数(i != j)，所以ei % mj == 0，所以x % ai，则只剩下ei*ai这一项，所以等价。

所以只需求e1*a1 + e2*a2 + … + en*an (mod M)，即可求出x.

代码

#include <iostream>

using namespace std;
typedef long long LL;
const int maxn = 100010;
LL fact[maxn];
void getFact(LL p) {
fact[0] = 1;
for(int i = 1; i <= p; i++) fact[i] = fact[i-1] * i % p;
}
LL pow_mod(LL a , LL n , LL m)
{
if(n == 0) return 1;
LL ret = pow_mod(a , n >> 1 , m);
ret = ret * ret % m;
if(n & 1) ret = ret * a % m;
return ret;
}
//Lucas算法依次得到a[i]
LL lucas(LL n , LL m , LL p)
{
getFact(p);
LL ret = 1;
while(n && m) {
LL a = n % p , b = m % p;
if(a < b) return 0;
ret = (ret * fact[a] * pow_mod(fact[b] * fact[a-b] % p , p-2 , p)) % p;
n /= p , m /= p;
}
return ret;
}
//扩展欧几里算法得求逆元
void ex_gcd(LL a , LL b , LL &d , LL &x , LL &y)
{
if(!b) {d = a; x = 1; y = 0;}
else {ex_gcd(b , a % b , d , y , x); y -= x * (a / b);}
}
//a * b % m
LL mul_mod(LL a , LL b , LL n) {
LL res = 0;
while(b) {
if(b & 1) res = (res + a) % n;
a = (a + a) % n;
b >>= 1;
}
return res;
}
//中国剩余定理合并答案
LL china(LL n , LL *a , LL *p) {
LL M = 1 , d  , y , x = 0;
for(int i = 0; i < n; i++) M *= p[i];
for(int i = 0; i < n; i++) {
LL w = M / p[i];
ex_gcd(p[i] , w , d , d , y);
//y * w * a[i]会超LL，所以要mul_mod
x = (x + mul_mod( y * w % M , a[i], M))%M;
}
return (x + M) % M;
}
int main()
{
LL m , n , k , t , p[20] , a[20];
cin >> t;
while(t--)
{
cin >> n >> m >> k;
for(int i = 0; i < k; i++) {
cin >> p[i];
a[i] = lucas(n , m , p[i]);
}
cout << china(k , a , p) << endl;
}
return 0;
}