【SDOI2013】方程

题意

给定方程
X1+X 2+…+Xn=m
我们对第 1.. n1 个变量 进行一些限制 ：
X1≤A1
X2≤A2
…
Xn1 ≤An1
我们对第 n1+1.. n1+1.. n1+ n2 个变量 进行一些限制 ：
X_(n1+1)≥A_(n1+1)
X_(n1+2)≥A_(n1+2)
…
X_(n1+n2) ≥A_(n1+n2)
求：在满足这些限制的前提下， 该方程正整数解的个数。
答案可能很大，请输出对 p取模 后的答案 ，也即 答案除以 p的余数。

分析

我们先考虑特殊情况，若是没有任何约束的话，就是把m个东西放到n个盒子里，
那么我们就可以用组合数答案就是C(m,n),
接着我们考虑其他情况，若是
X_(n1+1)≥A_(n1+1)
X_(n1+2)≥A_(n1+2)
…
X_(n1+n2) ≥A_(n1+n2)
这种情况的话，其实也很好处理，
我们观察到，X的取值也就是拿的球数是一定要先取A个的，所以我们可以用m-A。
若只有这种情况的话就是：C(m-$\sum ai$,n)
那么X1≤A1，X2≤A2这种情况我们可以用容斥原理解决：
原条件是xi<=ai，那么反过来就是xi>ai，那么我们可以求出所有的方案再减去所有不满足的方案，就是xi>ai的方案，这样又变成我们熟悉的东西了：

procedure dfs(x,y,z:int64);
begin
    if x>n1 then begin
       ans:=ans;
       if z mod 2=0 then ans:=(ans+c(m-y-1,n-1))mod p
       else ans:=(ans-c(m-y-1,n-1)+p)mod p;exit;
    end;
    dfs(x+1,y,z);dfs(x+1,y+a[x],z+1);
end;
dfs(1,0,0);

好了，这样的话，我们就差如何求组合数取mod了。

中国剩余定理

这里引入中国剩余定理

我们可以把mod数看作$mo=m1^{b1}*m2^{b2}*m3^{b3}..*mn^{bn}$
所以答案C(M,N)mod mo就是公式里面的x，
我们知道m1..mn，现在我们只用求出C(M,N)mod mi求出每个ai，然后合并。

阶乘

现在我们考虑如何算一个分数模一个数mo，而且这个$mo=m^{b}$
我们可以先算出分子含有m的个数，和分母的，以及它们的差tot
然后答案就是分子除m外的积除以分母除m外的积乘上$m^{tot}$
所以我们只用考虑如何快速求一个数的阶乘除去含有m的数的积，以及含有的个数。
也就是$a!=b*m^{tot}$
那么我们先观察一个数的阶乘1..a，里面一定含有a/m个m，
也就是$m^{a/m}$，那么每个原本是:$1m，2m，3m..(a/m)*m$
被提取出一个m后，就是$1,2,3..a/m$
那么我们便可以继续递归下去。
那么其他的部分呢？
其他的部分也就是：
$1..m^b-1,m^b+1..2m^b-1,2m^b+1..3m^b-1$
我们发现这些数在mod $m^b$的意义下都是一样的$1..m^b-1$
设fac[i]表示前i个数出去m的数的积mod $m^b$
所以我们可以把这一段的数的积看成是$fac[m^b-1]^{a/m^b}$
再乘上一段剩下的积$fac[a mod (m^b)]$

总的架构就是：

$dfs(a)=dfs(a/m)*fac[a mod (m^b)]*fac[m^b-1]^{a/m^b}$
这里指的是除了有因子m的积

NUM a=NUM a/m+a/m (NUM x指1..x的积的除了因子m的积）

这样便可以解决问题了。

代码

var
    t,p,i,n1,n2,m,x,pn,n,j:longint;ans:int64;
    a:array[1..8] of int64;
    b,pi,ai,pw,mi,lnv,phi:array[1..100] of int64;
    tre:array[1..10,0..10000] of int64;
function power(x,y,z:int64):int64;
begin
    if y=0 then exit(1);
    if y=1 then exit(x)
    else begin
        power:=power(x,y div 2,z);
        power:=(power*power)mod z;
        if y mod 2<>0 then power:=(power*x)mod z;
    end;
end;
procedure make(x,y:int64;var x1,y1:int64);
var x2,y2:int64;
begin
    if x<pi[y] then begin
       x1:=tre[y,x];y1:=0;exit;
    end; y1:=x div pi[y];
    make(x div pi[y],y,x2,y2);
    x1:=((x2*tre[y,x mod pw[y]])mod pw[y])*power(tre[y,pw[y]-1],x div pw[y],pw[y]) mod pw[y];
    y1:=y1+y2;
end;
function get(x,y,z:int64):int64;
var x2,y2,x3,y3,x4,y4,xx:int64;
begin
    make(x,z,x2,y2);
    make(x-y,z,x3,y3);
    make(y,z,x4,y4);
    xx:=((x2*(power(x3,phi[z]-1,pw[z]))mod pw[z])*power(x4,phi[z]-1,pw[z])) mod pw[z];
    if y2-y3-y4>=ai[z] then xx:=0 else xx:=(xx*power(pi[z],y2-y3-y4,pw[z])) mod pw[z];
    exit(xx);
end;
function c(x,y:int64):int64;
var an:int64;i:longint;
begin
    if x=y then exit(1)
    else if x<y then exit(0)
    else if y=0 then exit(1);
    fillchar(b,sizeof(b),0);
    for i:=1 to pn do b[i]:=get(x,y,i);
    an:=0;
    for i:=1 to pn do begin
    an:=(an+(b[i]*mi[i]mod p)*lnv[i]) mod p;
    an:=an;
    end;
    exit(an);
end;
procedure dfs(x,y,z:int64);
begin
    if x>n1 then begin
       ans:=ans;
       if z mod 2=0 then ans:=(ans+c(m-y-1,n-1))mod p
       else ans:=(ans-c(m-y-1,n-1)+p)mod p;exit;
    end;
    dfs(x+1,y,z);dfs(x+1,y+a[x],z+1);
end;
begin
    readln(t,p);
    x:=p;
    for i:=2 to trunc(sqrt(p)) do
    if x mod i=0 then begin
       inc(pn);pi[pn]:=i;ai[pn]:=1;pw[pn]:=i;x:=x div i;
       while x mod i=0 do begin
           x:=x div i;
           pw[pn]:=pw[pn]*i;inc(ai[pn]);
       end;tre[pn,0]:=1;
       for j:=1 to pw[pn] do
       if j mod i=0 then tre[pn,j]:=tre[pn,j-1] else tre[pn,j]:=(tre[pn,j-1]*j)mod pw[pn];
       phi[pn]:=pw[pn]*(pi[pn]-1) div pi[pn];
       mi[pn]:=p div pw[pn];
       lnv[pn]:=power(mi[pn],phi[pn]-1,pw[pn]);
    end;
    if x>1 then begin
       inc(pn);pi[pn]:=x;ai[pn]:=1;pw[pn]:=x;
       phi[pn]:=pw[pn]-1;
       mi[pn]:=p div pw[pn];
       lnv[pn]:=power(mi[pn],phi[pn]-1,pw[pn]);tre[pn,0]:=1;
       for j:=1 to pw[pn] do
       if j mod pi[pn]=0 then tre[pn,j]:=tre[pn,j-1] else tre[pn,j]:=(tre[pn,j-1]*j)mod pw[pn];
    end;
    for t:=1 to t do begin
        readln(n,n1,n2,m);
        for i:=1 to n1 do begin
            read(a[i]);
        end;
        for i:=1 to n2 do begin read(x);dec(m,x-1);end;
        if n2=0 then writeln(c(m-1,n-1))
        else begin
             ans:=0;dfs(1,0,0);writeln(ans);
        end;
    end;
end.