[SDOI2013]方程

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 }

 

posted @ 2017-10-08 21:30  Z-Y-Y-S  阅读(442)  评论(0编辑  收藏  举报