HDU 5446 Unknown Treasure（Lucas定理+中国剩余定理）

Description
求C(n,m)%M
Input
第一行一整数T表示用例组数，每组用例输入三个整数n,m,k,之后k个不同的素数p1~pk表示M=p1*p2*…*pk
(T<=20,1<=m<=n<=10^18,1<=k<=10,M<=10^18,pi<=10^5)
Output
对于每组用例，输出C(n,m)%M
Sample Input
1
9 5 2
3 5
Sample Output
6
Solution
设x=C(n,m)，因n和mm都比较大故可以由lucas定理得到C(n,m)%pi(1<=i<=k)的值ai，进而得到一个模方程组x mod pi=ai，因模数互素由中国剩余定理可以得到该方程组的通解，即x=t*M+x0(0<=x0< M)，x0即为答案
Code

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
#define maxn 111111
ll rev[maxn];
ll power[maxn]; 
ll mod_pow(ll a,ll b,ll p)
{
    ll ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
void init(ll n,ll p)//预处理出n的阶乘及其阶乘的逆元 
{
    power[0]=1;
    for(int i=1;i<n;i++)
        power[i]=power[i-1]%p*i%p;
    rev[1]=1;
    for(ll i=0;i<n;i++) 
        rev[i]=mod_pow(power[i],p-2,p);
}
ll C(ll n,ll m,ll p)//求组合数 
{
    if(m>n||m<0||n<0)return 0;
    if(m==0||m==n) return 1;
    return power[n]*rev[m]%p*rev[n-m]%p;
}
//n,m较大时需要用到lucas定理
ll lucas(ll n,ll m,ll p)
{
    ll ans=1;
    while(n&&m)
    {
        ll a=n%p,b=m%p;
        if(a<b) return 0;
        ans=ans*C(a,b,p)%p;
        n/=p;
        m/=p;
    }
    return ans;
} 
ll mul(ll a,ll b,ll p)
{
    if(b<0)b=(b%p+p)%p;
    a%=p;
    if(a<b)swap(a,b);
    ll ans=0;
    while(b)
    {
        if(b&1)ans=(ans+a)%p;
        a=(a+a)%p;
        b>>=1;
    }
    return ans;
}
ll extend_gcd(ll a,ll b,ll &x,ll &y)
{ 
    ll d=a;
    if(b!=0)
    {
        d=extend_gcd(b,a%b,y,x);
        y-=(a/b)*x;
    }
    else 
    {
        x=1;
        y=0;
    }
    return d;
}
ll CRT(ll a[],ll m[],int n)//a数组表示余数,m数组表示模数,n表示方程个数 
{
    ll M=1,Mi,ans=0;
    ll x,y,d;
    for(int i=0;i<n;i++) 
        M*=m[i];
    for(int i=0;i<n;i++)
    {
        Mi=M/m[i];
        d=extend_gcd(Mi,m[i],x,y);
        ans=(ans+mul(mul(Mi,x,M),a[i],M))%M;
    }
    if(ans<0) ans+=M;
    return ans;
}
int T,k;
ll n,m,r[11],p[11];  
int main()
{
    scanf("%d",&T);
    while(T--)
    {
        scanf("%I64d%I64d%d",&n,&m,&k);
        for(int i=0;i<k;i++)
        {
            scanf("%d",&p[i]);
            init(p[i],p[i]);
            r[i]=lucas(n,m,p[i]);
        }
        ll ans=CRT(r,p,k);
        printf("%I64d\n",ans);
    }
    return 0;
}