您的位置:首页 > 其它

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

2015-09-14 19:21 507 查看
题目地址:HDU 5446

题意:求C(N,M)%P的结果,P为k个素数的乘积。(1≤m≤n≤10^18,1≤k≤10)

思路:先用Lucas定理求出a[]数组,再用CRT求出最终的结果x。比赛时这道题坑了俩小时,还是用的不互质的板子过的,不知道为什么用互质的CRT板子一直Wa,sad。

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <iostream>
#include <sstream>
#include <algorithm>
#include <set>
#include <queue>
#include <stack>
#include <map>
using namespace std;
typedef __int64 LL;
LL prime[30];
LL a[30];
LL modxp(LL a,LL b,LL mod)
{
LL res=1;
while(b>0) {
if(b&1) res=res*a%mod;
b=b>>1;
a=a*a%mod;
}
return res;
}
LL C(LL n, LL m, LL mod)
{
if(m>n) return (LL)0;
if(m==n) return (LL)1;
LL ans=1;
for(LL i=1; i<=m; i++) {
LL a=(n+i-m)%mod;
LL b=i%mod;
ans=ans*(a*modxp(b,mod-2,mod)%mod)%mod;
}
return ans;
}
LL Lucas(LL n, LL m,LL mod)
{
if(m==0) return (LL)1;
return C(n%mod,m%mod,mod)*Lucas(n/mod,m/mod,mod)%mod;
}
void exgcd(LL a, LL b, LL &x, LL &y)
{
if(b==0) {
x=1;
y=0;
return ;
}
exgcd(b, a % b, y, x);
y -= x*(a/b);
}
/*LL CRT(LL a[], LL m[], int n)
{
LL M=1;
LL ans=0;
for(int i=0; i<n; i++)
M*=m[i];
for(int i=0; i<n; i++) {
LL x, y;
LL Mi=M/m[i];
extend_Euclid(Mi,m[i],x,y);
ans=(ans+x%M*a[i]%M*Mi%M)%M;
}
if(ans<0) ans+=M;
return ans;
}*/
LL d,x,y;
void gcd(LL a,LL b)
{
if(!b) {
d = a ;
x = 1 ;
y = 0 ;
} else {
gcd(b,a%b);
swap(x,y);
x = -x ;
y = -y ;
y += (a/b)*x ;
}
return ;
}
int main()
{
int T,k;
LL n,m;
LL aa , bb ,flag, i;
LL res;
scanf("%d",&T);
while(T--) {
scanf("%I64d %I64d %d",&n,&m,&k);
if(m<=n) m=n-m;
for(int i=0; i<k; i++) {
scanf("%I64d",&prime[i]);
a[i]=Lucas(n,m,prime[i]);
}
aa=prime[0] ;
bb=a[0] ;
flag=1 ;
for(i=1; i<k; i++) {
gcd(aa,prime[i]);
if((a[i]-bb)%d != 0)
flag=0;
if(flag) {
x=(a[i]-bb)/d*x ;
y=prime[i] / d ;
x=( x%y + y )%y ;
bb=bb + x * aa ;
aa=aa*prime[i]/d ;
}
}
gcd(1,aa);
if( bb%d!=0 )
flag=0;
if( flag) {
x=(bb/d)*x ;
y=aa/d ;
x=(x%y+y)%y ;
}
printf("%I64d\n",x);
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: