组合数求法讲解

　　首先对于c(n,m)，如果n,m比较大的话，这个值会很大，为了简化变成复杂度，也是为了更好的求解，都要求求c(n,m) mod p的值，我们由最简单的问题慢慢提高难度。

　　要求求解c(n,m) mod p，p为质数，且n,m<=10^5，这时我们可以存储fac[i]=i! mod p.这样就可以直接用公式求了。

　　当n,m的数据范围为10^9，且p为质数的时候且<=10^5，显然没法用数组存储10^9的信息，但是我们可以记录10^5的阶乘信息。然后我们引入lucas定理。

c(n,m)=lucas(n,m)=lucas(n div p,m div p)*c(n mod p,m mod p)。这样我们可以递归的将10^9的数据范围缩小到每次的10^5。

lucas：

function lucas(x,y,p:int64):int64;
var
    a, b                    :int64;
begin
    if y=0 then exit(1);
    a:=x mod p;
    b:=y mod p;
    if a<b then exit(0) else lucas:=lucas(x div p,y div p,p)*combine(a,b,p);
end;   

　　当p不是质数，但是p为互不相同的质数的乘积，且每个质数小于10^5候，我们可以将p分解质因数。分解成k个质数相乘。那么ans=c(n,m) mod p，我们可以得到c(n,m) mod pi=ansi这样的k个方程组，每个方程组解的方法与上一问题相同。利用crt（中国剩余定理）可以将这k个方程组合并，得出原问题的答案。

　　当p为若干质数的乘积，且pi^c<10^5时，我们可以由于各个pi^c互质，所以我们假设能求出c(n,m) mod pi^c这k个等式的解，就可以类似上一问题用crt合并出当前问题的解。那么问题就转化成了求解c(n,m) mod pi^c 因为mod的数是非质数，所以无法类似第一个方法求解（因为没有办法求逆元），观察答案形式，ans=c(n,m)=n!/(m!*(n-m)!)，我们可以提出来分子和分母中含有pi因子的项的pi因子，上下可以消掉，剩下的就是和pi^ci互质的数，就可以求逆元了。

　　最后一个问题bzoj 2142 http://61.187.179.132/JudgeOnline/problem.php?id=2142

//By BLADEVIL
var
    m, n                                 :longint;
    pj, c                               :array[0..110] of longint;
    pi, s, a                            :array[0..110] of int64;
    p                                   :int64;
    tot                                 :longint;
 
procedure divide(p:int64);
var    
    i, j                                :longint;
     
begin
    tot:=0;
    for i:=2 to trunc(sqrt(p)) do
    if p mod i=0 then
    begin
        inc(tot);
        pj[tot]:=i;
        while p mod i=0 do
        begin
                inc(c[tot]);
                p:=p div i;
        end;
    end;
    if p>1 then
    begin
        inc(tot);
        pj[tot]:=p;
        c[tot]:=1;
    end;
    for i:=1 to tot do
    begin
        pi[i]:=1;
        for j:=1 to c[i] do pi[i]:=pi[i]*pj[i];
    end;
end;
 
function ex_gcd(a,b:int64;var x,y:int64):int64;
var    
    t                                   :int64;
begin
    if (b=0) then
    begin
        x:=1;y:=0;
        exit(a);
    end;
    ex_gcd:=ex_gcd(b,a mod b,x,y);
    t:=x;
    x:=y;
    y:=t-(a div b)*y;
end;
 
function gcd(a,p:int64):int64;
var    
    x, y                                :int64;
     
begin
    x:=0;y:=0;
    ex_gcd(a,p,x,y);
    x:=(x mod p+p)mod p;
    exit(x);
end;
 
function mi(x,y,q:int64):int64;
var
    rec                                 :int64;
     
begin
    rec:=1;
    while (y>0) do
    begin
        if y and 1=1 then rec:=rec*x mod q;
        x:=x*x mod q;
        y:=y shr 1;
    end;
    exit(rec);
end;
 
function fac(n,p,q:int64):int64;
var
    cnt                                 :int64;
    i                                    :longint;
begin
    cnt:=1;
    for i:=1 to n do
        if (i mod p>0) then
            cnt:=cnt*i mod q;
    exit(cnt);
end;
 
function fact(n:int64;var sum:int64;p,q:int64):int64;
var
    cnt, rec                            :int64;
     
begin
    rec:=1;
    cnt:=fac(q,p,q);
    while n>=p do
    begin
        sum:=sum+n div p;
        if (n div q>0) then rec:=rec*(mi(cnt,n div q,q) mod q)mod q;
        if (n mod q>0) then rec:=rec*(fac(n mod q,p,q)mod q) mod q;
        n:=n div p;
    end;
    if n>1 then rec:=rec*fac(n,p,q) mod q;
    exit(rec);
end;
 
function combine(n,m,p,q:int64):int64;
var    
    ans1, ans2, ans3, ans               :int64;
    a, b, c                             :int64;
     
begin
    a:=0;b:=0;c:=0;
    ans1:=fact(n,a,p,q);
    ans2:=fact(m,b,p,q);
    ans3:=fact(n-m,c,p,q);
    a:=a-(b+c);
    ans:=mi(p,a,q);
    ans:=ans*ans1 mod q;
    ans:=ans*gcd(ans2,q) mod q;
    ans:=ans*gcd(ans3,q) mod q;
    exit(ans);
end;
 
function doit(n,m:longint):int64;
var
    i                                   :longint;
    x, y, sum                           :int64;
         
begin
    sum:=0;
    for i:=1 to tot do
            a[i]:=combine(n,m,pj[i],pi[i]);
    for i:=1 to tot do
    begin
        x:=0;y:=0;
        ex_gcd(s[i],pi[i],x,y);
        x:=(x mod pi[i]+pi[i])mod pi[i];
        sum:=(sum+((x*s[i] mod p)*a[i])mod p)mod p;
    end;
    exit(sum mod p);
end;
 
procedure main;
var    
    i                                   :longint;
    w                                   :array[0..100] of longint;
    ans                                 :int64;
    sum                                 :int64;
     
begin
    readln(p);
    divide(p);
    for i:=1 to tot do s[i]:=p div pi[i];
    readln(n,m);
    sum:=0;
    for i:=1 to m do
    begin
        readln(w[i]);
        inc(sum,w[i]);
    end;
    if sum>n then
    begin
        writeln('Impossible');
        exit;
    end;
    ans:=1;
    for i:=1 to m do
    begin
        ans:=ans*doit(n,w[i]) mod p;
        n:=n-w[i];
    end;
    writeln(ans);
end;
 
begin
    main;
end.