HDU 5446 Unknown Treasure(Lucas定理+CRT)

题目链接：传送门 

题意：

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

分析：

用Lucas求 A1 = C(n,m)%p1,A2 = C(n,m)%p2,...,Ak
= C(n,m)%pk

则用同余方程x = A1 (mod p1) ,x = A2
(mod p2) ,...,x = Ak(mod pk) ,然后用CRT解一下就好了。

代码如下：

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <cmath>
using namespace std;
typedef long long LL;

const int maxn = 20;

LL a[maxn];
LL M[maxn] ;

LL mul(LL a,LL b,LL mod){
LL ans = 0;
while(b){
if(b&1) ans=(ans+a)%mod;
b>>=1;
a=(a+a)%mod;
}
return ans;
}

LL quick_mod(LL a,LL b,LL m)
{
LL ans = 1;
a %= m;
while(b)
{
if(b & 1)
{
ans = ans * a % m;
b--;
}
b >>= 1;
a = a * a % m;
}
return ans;
}

LL C(LL n, LL m,int cur)
{
LL p = M[cur];
if(m > n) return 0;
if(m>n-m)
m= n-m;
LL ans = 1;
for(LL i=1; i<=m; i++)
{
LL a = (n + i - m) % p;
LL b = i % p;
ans = mul(ans , mul(a , quick_mod(b, p-2,p) , p), p); //p为素数，i对p的逆元可以不用扩展欧几里得进行求解 re=i^(p-2)
}
return ans%p;
}

LL Lucas(LL n,LL k,int cur)
{
LL p = M[cur];
if(k == 0) return 1%p;
return mul(C(n % p, k % p, cur) ,Lucas(n / p, k / p, cur) ,p);
}

void extend_Euclid(LL a, LL b, LL &x, LL &y)
{
if(b == 0)
{
x = 1;
y = 0;
return;
}
extend_Euclid(b, a % b,x, y);
LL tmp = x;
x = y;
y = tmp - a / b * y;
}

LL RemindChina(LL a[],LL m[],int k)
{
LL M = 1;
LL ans = 0;
for(int i=0; i<k; i++)
M *= m[i];
for(int i=0; i<k; i++)
{
LL x, y;
LL Mi = M / m[i];
extend_Euclid(Mi, m[i], x, y);
LL tmp;
if(x<0){
x=-x;
tmp = mul(Mi,x,M);
tmp = mul(tmp,a[i],M);
tmp = -tmp;
}
else {
tmp = mul(Mi,x,M);
tmp = mul(tmp,a[i],M);
}
ans = (ans + tmp) % M;
}
while(ans<0)
ans += M;
return ans;
}

int main()
{
int t;
scanf("%d",&t);
while(t--)
{
LL n,m;
int k;
scanf("%I64d%I64d%d",&n,&m,&k);
for(int i=0;i<k;i++){
scanf("%I64d",&M[i]);
}
for(int i=0;i<k;i++){
a[i] = Lucas(n,m,i)%M[i];
}
printf("%I64d\n",RemindChina(a,M,k));
}
return 0;
}