【BZOJ2142】礼物（扩展lucas定理，中国剩余定理合并方程）

题意：有n件礼物，m个人，每个人分别需要w[i]件礼物，求分礼物的不同方案数 mod P

提示：设P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi为质数。

1≤n≤10^9，1≤m≤5，1≤pi^ci≤10^5。

P不一定为质数

思路：经推导答案即为n!/(w[i]!)，i=1..n

考虑P不是质数

将P分解为提示中所说的形式，可以发现所有pt^ct都是互质的，所以我们可以用下图的中国剩余定理合并

From http://blog.csdn.net/popoqqq/article/details/39891263

然后对于每个pi^ai，我们进行以下处理：

将分子和分母化为x*pi^y的形式

然后分母的x部分与pi互质，可以求逆元，分子分母的y部分直接相减即可

然后我们处理阶乘

以19为例，将19!化为x*pi^y的形式，其中pi=3，ai=2 则有

19!%9=(1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19) %9

           =(1*2*4*5*7*8*10*11*13*14*16*17*19)*3^6*(1*2*3*4*5*6) %9

式子的左半部分是不为3的倍数的数，存在长度为p^a的循环节 求出一个循环节 快速幂处理 然后处理剩余部分

右半部分是6!%9 递归处理即可

用这样的思路就可以换解决C(n,m) mod p，p不为质数的情况了

var a:array[1..1000]of int64;
    n,m,i:longint;
    ans,s,mo:int64;

function mult(x,y,p:int64):int64;
var tmp:int64;
begin
 mult:=1; tmp:=x;
 while y>0 do
 begin
  if y and 1=1 then mult:=mult*tmp mod p;
  tmp:=tmp*tmp mod p;
  y:=y>>1;
 end;
end;

procedure exgcd(a,b:int64;var d,x,y:int64);
begin
 if b=0 then
 begin
  d:=a; x:=1; y:=0;
 end
  else
  begin
   exgcd(b,a mod b,d,y,x);
   y:=y-(a div b)*x;
  end;
end;

function inv(a,n:int64):int64;
var d,x,y:int64;
begin
 exgcd(a,n,d,x,y);
 if d=1 then exit((x+n) mod n);
 exit(-1);
end;

function fac(n,p,pr:int64):int64;
var i,re,r:int64;
begin
 if n=0 then exit(1);
 re:=1; i:=2;
 while i<=pr do
 begin
  if i mod p>0 then re:=re*i mod pr;
  inc(i);
 end;

 re:=mult(re,n div pr,pr);
 r:=n mod pr;
 i:=2;
 while i<=r do
 begin
  if i mod p>0 then re:=re*i mod pr;
  inc(i);
 end;

 exit(re*fac(n div p,p,pr) mod pr);
end;

function c(n,m,p,pr:int64):int64;
var x,y,z,t,tmp:int64;
begin
 if n<m then exit(0);
 x:=fac(n,p,pr);
 y:=fac(m,p,pr);
 z:=fac(n-m,p,pr);
 c:=0;
 t:=n;
 while t>0 do
 begin
  c:=c+t div p;
  t:=t div p;
 end;
 t:=m;
 while t>0 do
 begin
  c:=c-t div p;
  t:=t div p;
 end;
 t:=n-m;
 while t>0 do
 begin
  c:=c-t div p;
  t:=t div p;
 end;
 tmp:=x*inv(y,pr) mod pr*inv(z,pr) mod pr*mult(p,c,pr) mod pr;
 exit(tmp*(mo div pr) mod mo*inv(mo div pr,pr) mod mo);
end;


function lucas(n,m:int64):int64;
var x,re,i,pr:int64;
begin
 i:=2; x:=mo; re:=0;
 while i<=x do
 begin
  if x mod i=0 then
  begin
   pr:=1;
   while x mod i=0 do
   begin
    x:=x div i; pr:=pr*i;
   end;
   re:=(re+c(n,m,i,pr)) mod mo;
  end;
  inc(i);
 end;
 exit(re);
end;

begin
 assign(input,'bzoj2142.in'); reset(input);
 assign(output,'bzoj2142.out'); rewrite(output);
 readln(mo);
 readln(n,m);
 for i:=1 to m do
 begin
  read(a[i]); s:=s+a[i];
 end;
 if s>n then
 begin
  writeln('Impossible');
  close(input);
  close(output);
  exit;
 end;
 ans:=1;
 for i:=1 to m do
 begin
  ans:=ans*lucas(n,a[i]) mod mo;
  n:=n-a[i];
 end;
 writeln(ans);

 close(input);
 close(output);
end.