HDU 5446 lucas定理 + 中国剩余定理

http://acm.hdu.edu.cn/showproblem.php?pid=5446

求C(N,M)%MOD

n,m<=1e18

mod《1e18

但是mod由不超过10个不同素数组成，是个合数

因此：

令X=c(n,m)

分别求出 

A1=X%m1

A2=X%m2

....

得到同余方程

x==(A1)%m1

x==(A2)%m2

.....

用CRT解出X

由于数据比较大，会爆LL，需要高精度乘

#include<bits/stdc++.h>
using namespace std;
const int maxn = 100010;
long long fac[maxn];
long long a[20],m[20];
long long n,mm;
int k;

void extend_Euclid(long long a, long long b, long long &x, long long &y)
{
if(b == 0)
{
x = 1;
y = 0;
return;
}
extend_Euclid(b, a % b, x, y);
long long tmp = x;
x = y;
y = tmp - (a / b) * y;
}
long long fast_multi(long long m, long long n, long long mod)//快速乘法
{
long long ans = 0;//注意初始化是0，不是1
while (n)
{
if (n & 1)
ans += m;
m = (m + m) % mod;//和快速幂一样，只不过这里是加
m %= mod;//取模，不要超出范围
ans %= mod;
n >>= 1;
}
return ans;
}
long long CRT(long long a[],long long m[],int n,long long & M)
{
M = 1;
long long ans = 0;
for(int i=1; i<=n; i++)
M *= m[i];
for(int i=1; i<=n; i++)
{
long long x, y;
long long Mi = M / m[i];
extend_Euclid(Mi, m[i], x, y);
x = (x%m[i]+m[i])%m[i];
//ans = (ans + ((Mi % M) * (x % M) % M* a[i]) % M) % M;
long long tmp=fast_multi(Mi,x,M);
ans = (ans + fast_multi(tmp, a[i],M)) % M;

}
if(ans < 0) ans += M;
return ans;
}
void init(long long p)
{
fac[0] =1;
for(int i =1; i <= p; i++)
fac[i] = fac[i-1]*i % p;
}

long long fast_m(long long a, long long b,long long p)
{
long long tmp = a % p, ans =1;
while(b)
{
if(b &1) ans = ans * tmp % p;
tmp = tmp*tmp % p;
b >>=1;
}
return ans;
}
long long Lucas(long long n,long long m,long long p)
{
long long ret=1;
while(n&&m)
{
long long a=n%p,b=m%p;
if(a<b) return 0;
ret=(ret*fac[a]*fast_m(fac[b]*fac[a-b]%p,p-2,p))%p;
n/=p;
m/=p;
}
return ret;
}

int main()
{
//freopen("input.txt","r",stdin);

int tt;scanf("%d",&tt);
while (tt--)
{
scanf("%lld%lld%d",&n,&mm,&k);
for (int i=1; i<=k; i++)
{
scanf("%lld",&m[i]);
init(m[i]);
a[i]=Lucas(n,mm,m[i]);
}
long long M;
long long ans=CRT(a,m,k,M);
printf("%lld\n",ans);
}
return 0;
}