C(n+m,m) mod p的一类算法

Lucas定理

　　A、B是非负整数，p是质数。AB写成p进制：A=a[n]a[n-1]...a[0]，B=b[n]b[n-1]...b[0]。

　　则组合数C(A,B)与C(a[n],b[n])*C(a[n-1],b[n-1])*...*C(a[0],b[0]) modp同

　　即：Lucas(n,m,p)=c(n%p,m%p)*Lucas(n/p,m/p,p)

　　以求解n! % p为例，把n分段，每p个一段，每一段求的结果是一样的。但是需要单独处理每一段的末尾p, 2p, ...，把p提取出来，会发现剩下的数正好又是(n / p)!，相当于划归成了一个子问题，这样递归求解即可。这个是单独处理n!的情况，当然C(n,m)就是n!/(m!*(n-m)!)，每一个阶乘都用上面的方法处理的话，就是Lucas定理了，注意这儿的p是素数是有必要的。

Lucas最大的数据处理能力是p在10^5左右，不能再大了，hdu 3037就是10^5级别的！

　　未写main函数

const maxn=;

var n,m:int64;

    fac:array [..maxn] of int64;

function quickmod(a,b,mol:int64):int64;

var ans:int64;

begin

    ans:=;

    while b<> do

        begin

            if (b and )= then ans:=(ans*a) mod mol;

            a:=a*a mod mol;

            b:=b>>;

        end;

    exit(ans);

end;

function get_fact(p:int64):int64;

var i:longint;

begin

    fac[]:=;

    for i:= to p do

    fac[i]:=(fac[i-]*i) mod p;

end;

function lucas(n,m,p:int64):int64;

var ans,aa,bb:int64;

begin

    ans:=;

    while (a>) and (k>) do

        begin

            aa:=a mod p;

            bb:=b mod p;

            if aa<bb then exit();

            ans:=ans*fac[aa]*quickmod(fac[bb]*fac[aa-bb] mod p,p-,p) mod p;

            a:=a div p;

            k:=k div p;

        end;

    exit(ans);

end;

　　当p很大怎么搞？显然这样直接开数组会MLE

　　我们可以这样做。

　　分子相乘取模，分母相乘取模

　　对分母求逆元，分子乘逆元取模

　　给个代码

const mol=;

var n,m:longint;

function inv(a,p:int64):int64;

var b,c,q,k1,k2,k3:int64;

begin

  b:=p;

  c:=a mod b;

  q:=a div b;

  k1:=;

  k2:=;

  k3:=(k1+p-q*k2 mod p) mod p;

  while (c xor )<> do

     begin

       a:=b;

       b:=c;

       c:=a mod b;

       q:=a div b;

       k1:=k2;

       k2:=k3;

       k3:=(k1+p-q*k2 mod p) mod p;

     end;

  exit(k3);

end;

function calc(x,y:int64):int64;

var i,j,ans,up,down:int64;

begin

    i:=y+;

    j:=x;

    up:=;

    down:=;

    while i<=j do

        begin

            up:=up*i mod mol;

            inc(i);

        end;

    i:=x-y;

    while i>= do

        begin

            down:=down*i mod mol;

            dec(i);

        end;

    ans:=inv(down,mol);

    ans:=up*ans mod mol;

    exit(ans);

end;

begin

    read(n,m);

    writeln(calc(n+m,m));

end.