HDOJ-5446/2015 ACM/ICPC Asia Regional Changchun Online 1010(数论)

题意：求组合数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.csdn.net/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.csdn.net/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
;
    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;
}