SDOI2013 方程
传送门
看到这个题想起来熟悉的插板(正整数解,不能为空)。
首先看后一种限制,这个好办,强制先往这些盒子里放一些球,之后再插板。
前一种限制怎么办……? 那我们考虑反着做,用容斥原理计算出不合法的情况,减掉就可以了。
想起来还是很简单的,但是这玩意实现起来真是恶心的要死…… 模数可以分解,于是我全部预处理了阶乘,不然是真的过不去……
看一下暴力的代码。
// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define rep(i,a,n) for(ll i = a;i <= n;i++)
#define per(i,n,a) for(ll i = n;i >= a;i--)
#define enter putchar('\n')
using namespace std;
typedef long long ll;
const int M = 200005;
const int INF = 2147483647;
const ll mod1 = 262203414;
const ll mod2 = 437367875;
const double eps = 1e-8;
ll read()
{
ll ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
ll T,p,n,n1,n2,m,A[20],sum,ans,fac[M],inv[M];
ll bas[10] = {0,2,3,11,397,10007};
ll mul(ll a,ll b,ll t)
{
ll res = a * b - (ll)((long double)a / t * b + eps) * t;
return (res % t + t) % t;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b){x = 1,y = 0;return a;}
ll d = exgcd(b,a%b,y,x);
y -= a / b * x;
return d;
}
ll Inv(ll a,ll b)
{
ll x,y;
exgcd(a,b,x,y);
return (x % b + b) % b;
}
ll CRT(ll b,ll t) {return mul(mul(b,Inv(p/t,t),p),p/t,p);}
ll qpow(ll a,ll b,ll t)
{
ll p = 1;
while(b)
{
if(b & 1) p = mul(p,a,t);
a = mul(a,a,t),b >>= 1;
}
return p;
}
void init(ll p,ll pk)
{
fac[0] = 1;
rep(i,1,pk) (i%p) ? fac[i] = mul(fac[i-1],i,pk) : fac[i] = fac[i-1];
}
ll Fac(ll x,ll pi,ll pk)
{
if(!x) return 1;
ll res = fac[pk];
res = mul(qpow(res,x/pk,pk),fac[x%pk],pk);
return mul(res,Fac(x/pi,pi,pk),pk);
}
ll C(ll n,ll m,ll pi,ll pk)
{
init(pi,pk);
ll d = Fac(n,pi,pk),d1 = Fac(m,pi,pk),d2 = Fac(n-m,pi,pk);
ll k = 0;
for(ll i = n;i;i /= pi) k += i / pi;
for(ll i = m;i;i /= pi) k -= i / pi;
for(ll i = n-m;i;i /= pi) k -= i / pi;
return mul(mul(d,Inv(d1,pk),pk),mul(qpow(pi,k,pk),Inv(d2,pk),pk),pk);
}
ll exlucas(ll n,ll m)
{
if(n < m) return 0;
ll res = CRT(C(n,m,p,p),p);
return res;
}
ll exlucas1(ll n,ll m)
{
if(n < m) return 0;
ll res = 0;
rep(i,1,5) res += CRT(C(n,m,bas[i],bas[i]),bas[i]),res %= p;
return res;
}
ll exlucas2(ll n,ll m)
{
if(n < m) return 0;
ll res = 0;
res += CRT(C(n,m,5,125),125),res %= p;
res += CRT(C(n,m,7,343),343),res %= p;
res += CRT(C(n,m,101,10201),10201),res %= p;
return res;
}
int get(int x)
{
sum = 0;
int cur = 1,tot = 0;
while(x)
{
if(x&1) sum += A[cur],tot++;
x >>= 1,cur++;
}
return tot;
}
int main()
{
T = read(),p = read();
while(T--)
{
n = read(),n1 = read(),n2 = read(),m = read();
rep(i,1,n1+n2) {A[i] = read();if(i > n1) m -= (A[i]-1);}
if(p == 10007)
{
ans = exlucas(m-1,n-1);
rep(i,1,(1<<n1)-1)
{
int d = get(i);
if(!((d-1)&1)) ans += (p - exlucas(m-1-sum,n-1)),ans %= p;
else ans += exlucas(m-1-sum,n-1),ans %= p;
}
printf("%lld\n",ans);
}
if(p == 262203414)
{
ans = exlucas1(m-1,n-1);
rep(i,1,(1<<n1)-1)
{
int d = get(i);
if(!((d-1)&1)) ans += (p - exlucas1(m-1-sum,n-1)),ans %= p;
else ans += exlucas1(m-1-sum,n-1),ans %= p;
}
printf("%lld\n",ans);
}
if(p == 437367875)
{
ans = exlucas2(m-1,n-1);
rep(i,1,(1<<n1)-1)
{
int d = get(i);
if(!((d-1)&1)) ans += (p - exlucas2(m-1-sum,n-1)),ans %= p;
else ans += exlucas2(m-1-sum,n-1),ans %= p;
}
printf("%lld\n",ans);
}
}
return 0;
}
当你意识到,每个上一秒都成为永恒。