首先我们设 Bi=max(A−Di,1) ∑mi=1∑Xj≤Hj≤Yj[∑nj=1[⌈HjBj⌉≤i]≥k] 这让我们很不好搞，但是我们可以容斥，枚举若干个，然后要求这些一定要击败，其余则随意。容斥系数见博客置顶文章容斥的原理及广义应用中的方法求出。接下来我们用f表示容斥系数。 ∑mi=1∑Sf(|S|)∑Xj≤Hj≤YjΠj∈S[⌈HjBj⌉≤i] ∑mi=1∑Sf(|S|)Πnj=1([j∉S](Yj−Xj+1)+[j∈S]∑YjHj=Xj[⌈HjBj⌉≤i]) ∑mi=1∑Sf(|S|)Πnj=1([j∉S](Yj−Xj+1)+[j∈S]max(min(Yj,Bji)−Xj+1),0)) 观察最后那个又有max又有min的东西，现在可以注意到对于每个j可以找到两个i作为分界点，在不同的i值域区间里最后那个式子会有不同的取值。 我们把所有关键点提取出来排序，然后现在关键点之间的区间里每个j影响固定了，我们考虑怎么计算答案。 我们设ff[i,j,k]表示考虑到第i个怪物，目前往S里选了j个怪物（这一维方便我们最后乘上容斥系数），然后有k个过了第一个分界点但是没过第二个分界点的（此时这样的怪物贡献是 Bji−Xj+1 ，可以把它给展开掉）,我们选择了k个 Bji (当然我们还没有乘上i)。这一维我们最后需要给它乘一个自然数幂和。
#include<cstdio> #include<algorithm> #define fo(i,a,b) for(i=a;i<=b;i++) #define fd(i,a,b) for(i=a;i>=b;i--) using namespace std; typedef long long ll; const int maxn=50+10,mo=1000000007; int f[maxn],ff[maxn][maxn][maxn],c[maxn][maxn],su[maxn][maxn],S[maxn],g[maxn],gg[maxn]; int x[maxn],y[maxn],d[maxn],b[maxn]; struct dong{ int x,y; friend bool operator <(dong a,dong b){ return a.x<b.x||a.x==b.x&&a.y>b.y; } } p[maxn*2]; bool bz[maxn],pd[maxn]; int i,j,k,l,t,n,m,q,a,tot,top,ans; void getS(int n,int m){ int i,j,k,t; if (n<0){ fo(i,0,m) S[i]=0; return; } S[0]=n; fo(i,1,m){ t=1; fd(j,n+1,n-i+1) if (j%(i+1)==0) t=(ll)t*(j/(i+1))%mo;else t=(ll)t*j%mo; fo(j,0,i-1){ k=(ll)su[i][j]*S[j]%mo; if ((i+j)%2==0) (t-=k)%=mo;else (t+=k)%=mo; } S[i]=t; } } void work(int l,int r){ int i,j,k,t; fo(i,0,n) fo(j,0,n) fo(k,0,n) ff[i][j][k]=0; ff[0][0][0]=1; fo(i,0,n-1) fo(j,0,i) fo(k,0,j) if (ff[i][j][k]){ (ff[i+1][j][k]+=(ll)ff[i][j][k]*(y[i+1]-x[i+1]+1)%mo)%=mo; if (pd[i+1]){ if (bz[i+1]) (ff[i+1][j+1][k]+=(ll)ff[i][j][k]*(y[i+1]-x[i+1]+1)%mo)%=mo; else{ (ff[i+1][j+1][k]+=(ll)ff[i][j][k]*(1-x[i+1])%mo)%=mo; (ff[i+1][j+1][k+1]+=(ll)ff[i][j][k]*b[i+1]%mo)%=mo; } } } getS(l-1,n); fo(i,0,n) g[i]=S[i]; getS(r,n); fo(i,0,n) gg[i]=S[i]; fo(j,0,n) fo(k,0,j) (ans+=(ll)ff[n][j][k]*f[j]%mo*(gg[k]-g[k])%mo)%=mo; } int main(){ scanf("%d%d%d%d",&n,&m,&q,&a); fo(i,1,n){ scanf("%d%d%d",&d[i],&x[i],&y[i]); b[i]=max(1,a-d[i]); } c[0][0]=1; su[0][0]=1; fo(i,1,n){ c[i][0]=1; fo(j,1,i){ c[i][j]=(c[i-1][j]+c[i-1][j-1])%mo; su[i][j]=(su[i-1][j-1]+(ll)su[i-1][j]*(i-1)%mo)%mo; } } fo(i,0,n){ t=0; fo(j,0,i-1) (t+=(ll)c[i][j]*f[j]%mo)%=mo; if (i>=q) f[i]=1-t;else f[i]=-t; } fo(i,1,n){ t=x[i]/b[i]; if (x[i]%b[i]) t++; p[++top].x=t; p[top].y=i; t=(y[i]+1)/b[i]; if ((y[i]+1)%b[i]) t++; p[++top].x=t; p[top].y=i; } p[++top].x=m; p[++top].x=1; p[top].y=-1; sort(p+1,p+top+1); fo(i,1,top){ if (p[i].y==0){ work(m,m); break; } if (p[i].y>0){ j=p[i].y; if (!pd[j]) pd[j]=1;else bz[j]=1; } if (p[i].x>=1&&p[i].x<p[i+1].x) work(p[i].x,p[i+1].x-1); } (ans+=mo)%=mo; printf("%d\n",ans); }