[SDOI2013]方程

Posted 2020-10-10 Z-Y-Y-S

Description

给定方程
    X1+X2+. +Xn=M
我们对第l..N1个变量进行一些限制：
Xl < = A
X2 < = A2
Xn1 < = An1
我们对第n1 + 1..n1+n2个变量进行一些限制：
Xn1+l > = An1+1
Xn1+2 > = An1+2
Xnl+n2 > = Anl+n2
求：在满足这些限制的前提下，该方程正整数解的个数。
答案可能很大，请输出对p取模后的答案，也即答案除以p的余数。

Input

    输入含有多组数据，第一行两个正整数T，p。T表示这个测试点内的数据组数，p的含义见题目描述。
    对于每组数据，第一行四个非负整数n，n1，n2，m。
    第二行nl+n2个正整数，表示A1..n1+n2。请注意，如果n1+n2等于0，那么这一行会成为一个空行。

Output

  共T行，每行一个正整数表示取模后的答案。

Sample Input

3 10007
3 1 1 6
3 3
3 0 0 5

3 1 1 3
3 3

Sample Output

3
6
0

【样例说明】
对于第一组数据，三组解为(1，3，2)，(1,4,1),(2,3,1)
对于第二组数据，六组解为(1，1，3)，(1,2,2),(1,3,1),(2,1,2),(2,2,1),(3,1,1)

HINT

n < = 10^9  , n1 < = 8   , n2 < = 8   ,  m < = 10^9  ,p<=437367875

对于l00%的测试数据：  T < = 5，1 < = A1..n1_n2  < = m，n1+n2 < = n

出题人很贱，故意使模数为合数

首先，是简单的排列组合

对于n2个条件，直接m-=ai-1就行了

对于满足n1条件，可以用容斥，枚举一个不满足减去，两个加上，以此类推

不满足条件直接把m-=ai，就是说直接分ai给i使i不符合

用隔版法，把m个物品分成n个非空子集：C(n-1,m-1)

由于n，m很大，所以只能lucas，但因为模数是合数，所以不可以

于是就有了拓展lucas

http://www.cnblogs.com/Przz/p/5409566.html

先分解模数Mod的质因数，Mod=p1^t1*p2^t2*.....

对于每个p=p1^t1取模，这样就可以用中国剩余定理(要求模数互质)

这样是为了方便快速阶乘取模

C(x,y)=y!*(x!)^-1*(y-x!)^-1

说明一下快速阶乘取模，因为p只含一个因数,拿p=3^2,即pi=3，ti=2，n=19举例

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

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

因为1×2×3×4×5×6×7×8≡10×11×12×13×14×15×16×17  (mod 9)

n!=(fac[9-1]^(19/9))*fac[19%9]*((19/3)!%p)

后面部分递归处理，至于红色部分，在阶乘前就算出x!,y!,(y-x)!中pi的含量a,b,c

用f(x)=x/pi+f[x/pi]递归求出

将答案乘pi^(a-b-c)

还有求出了阶乘，还要算逆元，由于p=pi^ti是合数

所以费马和线性模逆元公式都无法使用，但庆幸的是可以用拓展欧几里德

ax≡1(mod p) 的条件是a与p互质

意味着a只要含有pi就没有逆元

很简单，把所有pi的倍数事先挑出来(上面红色部分就是为了现在)

在阶乘前缀积中不予计算，这样算出来的阶乘就可以套用拓展欧几里德求逆元了

至于中国剩余定理这里不做说明

 

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cstring>
  4 #include<algorithm>
  5 using namespace std;
  6 typedef long long ll;
  7 ll fac[100011],tot,B[10001],d[10001],pri[10001],A[10001],c[10001],n,n1,n2,m;
  8 ll a[10001],Mod,ans;
  9 ll qpow(ll x,ll y,ll p)
 10 {
 11   ll res=1;
 12   while (y)
 13     {
 14       if (y&1) res=res*x%p;
 15       x=x*x%p;
 16       y/=2;
 17     }
 18   return res;
 19 }
 20 ll exgcd(ll a,ll b,ll &x,ll &y)
 21 {
 22   if (!b) 
 23     {
 24       x=1;y=0;
 25       return a;
 26     }
 27   ll d=exgcd(b,a%b,x,y);
 28   int t=x;x=y;y=t-(a/b)*y;
 29   return d; 
 30 }
 31 ll reverse(ll a,ll b)//求逆元
 32 {
 33   ll x,y;
 34   exgcd(a,b,x,y);
 35     return (x%b+b)%b;
 36 }
 37 ll calfac(ll x,ll p,ll t)//快速阶乘取模
 38 {
 39   if (x<t) return fac[x];
 40   ll s=qpow(fac[p-1],x/p,p);
 41   s=s*fac[x%p]%p;
 42   s=s*calfac(x/t,p,t)%p;
 43   return s;
 44 }
 45 ll calc(ll x,ll y,ll p,ll t)//组合数取模(拓展lucas)
 46 {int i;
 47   ll ap=0,bp=0,cp=0;
 48   fac[0]=1;
 49   for (i=1;i<=p-1;i++)
 50     if (i%t==0) fac[i]=fac[i-1];
 51     else fac[i]=fac[i-1]*i%p;
 52   for (i=y;i;i/=t) ap+=i/t;
 53   for (i=x;i;i/=t) bp+=i/t;
 54   for (i=y-x;i;i/=t) cp+=i/t;
 55   ap=ap-bp-cp;
 56   ll s=qpow(t,ap,p);
 57   ap=calfac(y,p,t);bp=calfac(x,p,t);cp=calfac(y-x,p,t);
 58   s=(((s*ap%p)*reverse(bp,p)%p)*reverse(cp,p)%p);
 59   return s;
 60 }
 61 ll C(ll x,ll y,ll p)//中国剩余定理
 62 {int i;
 63   if (x>y) return 0;
 64   for (i=1;i<=tot;i++)
 65     B[i]=calc(x,y,d[i],pri[i]);
 66   for (i=1;i<=tot;i++)
 67     A[i]=reverse(c[i],d[i]);
 68   ll s=0;
 69   for (i=1;i<=tot;i++)
 70     {
 71       s+=((B[i]*A[i]%p)*c[i]%p);
 72     s%=p;
 73     }
 74   return s;
 75 }
 76 void dfs(ll id,ll m,ll cnt)//容斥原理
 77 {
 78   if (id>n1)
 79     {
 80       ans=(ans+C(n-1,m-1,Mod)*cnt+Mod)%Mod;//隔版法求方案
 81       return;
 82     }
 83   dfs(id+1,m,cnt);
 84   if (m>=a[id]) dfs(id+1,m-a[id],-cnt);
 85 }
 86 int main()
 87 {ll T,x,i;
 88   cin>>T>>Mod;
 89   x=Mod;
 90   for (i=2;i*i<=Mod;i++)
 91     if (x%i==0)
 92       {int cnt=0;
 93     d[++tot]=1;
 94     pri[tot]=i;
 95       while (x%i==0)
 96     {
 97       ++cnt;
 98       x/=i;
 99       d[tot]*=i;
100     }
101 
102     }
103   if (x!=1)
104     {
105       tot++;
106       d[tot]=x;
107       pri[tot]=x;
108     }
109   for (i=1;i<=tot;i++) c[i]=Mod/d[i];
110   while (T--)
111     {
112       cin>>n>>n1>>n2>>m;
113       for (i=1;i<=n1;i++)
114       scanf("%lld",&a[i]);
115       for (i=1;i<=n2;i++)
116     {
117       scanf("%lld",&x);
118       if (x) m-=x-1;
119     }
120       if (m<n1)
121         cout<<0<<endl;
122       else 
123     {
124       ans=0;
125       dfs(1,m,1);
126       cout<<ans<<endl;
127     }
128     }
129 }