题意:求组合数C(n,m)mod(p1*p2*p3......*pk), pi均为质数,1<=m<=n<=10^18,1<=k<=10,2<=pi<=10^5
思路:以前曾经做过一次n,m范围是10^6的(http://blog.youkuaiyun.com/moringrain/article/details/46793743),做法是对n!,m!,(n-m)!分解质因数。
但是由于这道题n,m的范围是在是太大了,所以此路不通。但是可以观察出这道题的取模很有意思,把模分解质因数给你了,质因数个数不多,而且每个质因数的大小相对来说还是不大的。看来要顺着这个思路去想。
第一步:假设已知C(n,m) mod pi,(1<=i<=k),是否可以求C(n,m)mod(p1*p2*p3......*pk)?想到了这个问题:(摘录自百度)
韩信点兵又称为中国剩余定理,相传汉高祖刘邦问大将军韩信统御兵士多少,韩信答说,每3人一列余1人、5人一列余2人、7人一列余4人、13人一列余6人…….刘邦茫然而不知其数.我们先考虑下列的问题:假设兵不满一万,每5人一列、9人一列、13人一列、17人一列都剩3人,则兵有多少?首先我们先求5、9、13、17之最小公倍数9945(注:因为5、9、13、17为两两互质的整数,故其最小公倍数为这些数的积),然後再加3,得9948(人).中国有一本数学古书「孙子算经」也有类似的问题:「今有物,不知其数,三三数之,剩二,五五数之,剩三,七七数之,剩二,问物几何?」 答曰:「二十三」 术曰:「三三数之剩二,置一百四十,五五数之剩三,置六十三,七七数之剩二,置三十,并之,得二百三十三,以二百一十减之,即得.凡三三数之剩一,则置七十,五五数之剩一,则置二十一,七七数之剩一,则置十五,即得.」 孙子算经的作者及确实着作年代均不可考,不过根据考证,着作年代不会在晋朝之後,以这个考证来说上面这种问题的解法,中国人发现得比西方早,所以这个问题的推广及其解法,被称为中国剩余定理.中国剩余定理(Chinese Remainder Theorem)在近代抽象代数学中占有一席非常重要的地位.
好了,这个问题可以解决了。
第二步:如何在不超时的前提下去计算C(n,m) mod pi?说来惭愧,自己刷题还太少,比赛最后半个小时翻学长留下来的数学模板才发现了“Lucas定理”,可以O(pi)的算出组合数,p要求为质数。
关于Lucas定理,可以参考下面的博客:
http://blog.youkuaiyun.com/wukonwukon/article/details/7341270
最后半个小时虽然搞定了思路,但是还是写挂了,赛后和学长交流才知道,中间过程long long *long long 在取模之前爆掉了,要用快速乘法。
所谓快速乘法,就是把乘法变成移位相加,比如6(D)*5(D)=110(B)*(101)(B)=1(B)*110(B)+0(B)*1100(B)+1(B)*11000(B),在移位和累加的时候取模就不会爆掉了。
上算法(抄的退役学长的板子,自己也该总结下了)
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
using namespace std;
const int MAXP=100005;
typedef long long ll;
ll f[MAXP];
ll mul(ll x, ll y, ll p)
{
ll res=0;
while(y)
{
if(y&1) res=(res+x)%p;
x=(x<<1)%p;
y>>=1;
}
return res;
}
void exgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
if (!b)
{
d=a;
x=1;
y=0;
}
else
{
exgcd(b,a%b,d,y,x);
y-=a/b*x;
}
}
void initp(int p)
{
f[0]=f[1]=1;
for(int i=2;i<p;i++)f[i]=(f[i-1]*i)%p;
}
ll comb(ll n,ll m,ll p)
{
if(m>n) return 0;
ll ans=f[n];
ll g,x,y;
exgcd(f[m],p,g,x,y);
//ans=(ans*(x%p)%p+p)%p;
ans=mul(ans,(x%p)%p+p,p);
exgcd(f[n-m],p,g,x,y);
//ans=(ans*(x%p)%p+p)%p;
ans=mul(ans,(x%p)%p+p,p);
return ans;
}
//n取m,如果可能m>n请特判return 0;
ll lucas(ll n,ll m,int p)
{
ll ans=1;
initp(p);
while(n && m && ans)
{
//ans=comb(n%p,m%p,p)*ans%p;
ans=mul(comb(n%p,m%p,p),ans,p);
n/=p;
m/=p;
}
return ans;
}
//解x=bi(mod ai),注意这里a[]两两互质
ll modx(ll a[], ll b[],int len)
{
ll d, x, y, m, n=1, ret=0;
for(int i=0;i<len;i++) n*=a[i];
for(int i=0;i<len;i++)
{
m=n/a[i];
exgcd(a[i],m,d,x,y);
ret=ret+mul(mul(y,m,n),b[i],n);
}
return (n+ret%n)%n;
}
void work()
{
ll a[20],b[20];
ll n,m,k;
scanf("%I64d%I64d%I64d",&n,&m,&k);
for(int i=0;i<k;i++) scanf("%I64d",&a[i]);
for(int i=0;i<k;i++) b[i]=lucas(n,m,a[i]);
ll ans=0;
ans=modx(a,b,k);
printf("%I64d\n",ans);
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
work();
return 0;
}