【bzoj3129】【SDOI2013】方程 【exLucas】【容斥】

洛谷連結

bzoj連結

題目大意

給一個不定方程 ∑ i = 1 n x i = m \sum_{i=1}^{n}x_i=m ∑i=1n​xi​=m( x i > 1 x_i>1 xi​>1， x i ∈ N + x_i\in N^+ xi​∈N+)，求方程的解的數量，并滿足：

當 0 &lt; i ≤ n 1 0&lt; i\leq n1 0<i≤n1時， x i ≤ a i x_i\leq a_i xi​≤ai​

當 n 1 &lt; i ≤ n 1 + n 2 n1&lt; i\leq n1+n2 n1<i≤n1+n2時， x i ≥ a i x_i\geq a_i xi​≥ai​

( n 1 + n 2 ≤ n ≤ 1 0 9 n1+n2\leq n\leq 10^9 n1+n2≤n≤109， m ≤ 1 0 9 m\leq10^9 m≤109， 0 ≤ n 1 , n 2 ≤ 8 0\leq n1,n2\leq 8 0≤n1,n2≤8)

答案對 P P P取模， P ∈ { 10007 , 262203414 , 437367875 } P\in \{10007,262203414,437367875\} P∈{10007,262203414,437367875}，多組詢問

解法

練習容斥，exLucas好題

若隻有第二種限制就很好辦，插闆法即可，答案是 ( m − ∑ i = n 1 n 1 + n 2 ( a i − 1 ) − 1 n − 1 ) \binom{m-\sum_{i=n1}^{n1+n2}(a_i-1)-1}{n-1} (n−1m−∑i=n1n1+n2​(ai​−1)−1​)

考慮到n1很小，可以容斥，用一共的答案-大于 a i a_i ai​的答案+兩兩重複部分-三三重複部分…(具體見代碼)

出題人毒瘤，模數不是質數，要用exLucas

btw： 262203414 = 2 ∗ 3 ∗ 11 ∗ 397 ∗ 1007 , 437367875 = 5 3 ∗ 7 3 ∗ 10 1 2 262203414=2*3*11*397*1007,437367875=5^3*7^3*101^2 262203414=2∗3∗11∗397∗1007,437367875=53∗73∗1012

代碼

//#pragma GCC optimize(3)
#include<iostream>
#include<iomanip>
#include<stack>
#include<queue>
#include<list>
#include<vector>
#include<set>
#include<string>
#include<algorithm>
#include<complex>
#include<cmath>
#include<cstdio>
#include<cstring>
#define ull unsigned long long
#define ll long long
#define db long double
#define inf 100002
#define infm 1000009
#define INF (int)1e9
//#define mod 100000007LL
#define pi acos(-1)
#define rd(n) {n=0;char ch;int f=0;do{ch=getchar();if(ch=='-'){f=1;}}while(ch<'0'||ch>'9');while('0'<=ch&&ch<='9'){n=(n<<1)+(n<<3)+ch-48;ch=getchar();}if(f)n=-n;}
using namespace std;

ll Qpow(ll x,ll k,ll mod){
    ll ans=1LL;
    while(k){
        if (k&1){
            ans=ans*x%mod;
        }
        x=x*x%mod;
        k/=2;
    }
    return ans;
}

ll exgcd(ll n,ll m,ll &X,ll &Y){
    if (!m){
        X=1,Y=0;
        return n;
    }
    ll g=exgcd(m,n%m,Y,X);
    Y-=X*(n/m);
    return g;
}

ll getinv(ll x,ll mod){
	ll X,Y;
	exgcd(x,mod,X,Y);
    return (X%mod+mod)%mod;
}

ll f[inf];

ll fact(ll x,ll base,ll mod){
    if (x<base){
    	return f[x];
	}
    ll ans=Qpow(f[mod-1],x/mod,mod)*f[x%mod]%mod*fact(x/base,base,mod)%mod;
    return ans;
}

ll C(ll n,ll m,ll base,ll mod){
	f[0]=1;
	for (int i=1;i<=mod;i++){
		if (i%base!=0){
			f[i]=f[i-1]*i%mod;
		}
		else{
			f[i]=f[i-1];
		}
	}
    ll cnt=0;
    for (int i=n;i;i/=base){
        cnt+=i/base;
    }
    for (int i=m;i;i/=base){
        cnt-=i/base;
    }
    for (int i=n-m;i;i/=base){
        cnt-=i/base;
    }
    return Qpow(base,cnt,mod)*fact(n,base,mod)%mod*getinv(fact(m,base,mod),mod)%mod*getinv(fact(n-m,base,mod),mod)%mod;
}

ll mod;
ll ele[10],tms[10];
int ecnt;
ll w[10];

ll Cal(ll n,ll m){
	if (n<m){
		return 0;
	}
    for (int i=1;i<=ecnt;i++){
        ll P=Qpow(ele[i],tms[i],mod+1);
        w[i]=C(n,m,ele[i],P);
    }
    ll ans=0;
    for (int i=1;i<=ecnt;i++){
        ll P=Qpow(ele[i],tms[i],mod+1);
        ll tmp=mod/P,X,Y;
        exgcd(tmp,P,X,Y);
        ans=(ans+X*tmp%mod*w[i]%mod)%mod;
    }
    return (ans+mod)%mod;
}

ll n,m;
ll a[10];
int cnt;

ll IEP(int pos,ll now,int f){
    if (pos==cnt+1){
        return Cal(now-1,n-1)*f;
    }
    ll ans=0;
    ans=IEP(pos+1,now,f);
    if (now>=a[pos]){
    	ans=(ans+IEP(pos+1,now-a[pos],-f))%mod;
	}
    return ans;
}

int main(){
    int task;
    rd(task) rd(mod)
    if (mod==10007){
        ecnt=1;
        ele[1]=mod;
        tms[1]=1;
    }
    else if (mod==262203414){
        ecnt=5;
        ele[1]=2,ele[2]=3,ele[3]=11,ele[4]=397,ele[5]=10007;
        tms[1]=1,tms[2]=1,tms[3]=1,tms[4]=1,tms[5]=1;
    }
    else{
        ecnt=3;
        ele[1]=5,ele[2]=7,ele[3]=101;
        tms[1]=3,tms[2]=3,tms[3]=2;
    }
    ll n1,n2;
    for (int t=1;t<=task;t++){
        rd(n) rd(n1) rd(n2) rd(m)
        ll x;
        for (int i=1;i<=n1;i++){
            rd(a[i])
        }
        for (int i=1;i<=n2;i++){
            rd(x)
            if (x){
            	m-=x-1;
			}
        }
        if (m<n1){
        	puts("0");
        	return 0;
		}
        cnt=n1;
        printf("%lld\n",(IEP(1,m,1)+mod)%mod);
    }
    return 0;
}