题解 LOJ2026 「JLOI / SHOI2016」成绩比较

稍微改动一下题面中定义的变量名。我们定义每门课的最高分为(h_i)（也就是题面中的(u_i)）。定义被B神碾压的人数为(K)（而不是题面中的小写(k)，这个小写(k)我们下面会有别的用处）。

题目要求恰好有(K)位同学被B神碾压。设答案为(f_K)。我们可以先求出至少有(k)位同学被碾压的情况数(g_{k})。那么根据二项式反演（也就是简单容斥）可知：

[f_{K}=sum_{k=K}^{n}(-1)^{k-K}{kchoose K}g_k ]

考虑求(g_k)。首先，当(k>min_{i=1}^{m}(n-r_i))时，显然不可能有(k)位同学被碾压，因此此时(g_k=0)。

其他的情况。我们不妨先写出一个朴素的式子，再考虑优化其复杂度。

首先钦定除B神外的(n-1)位同学中的(k)位是被碾压的。然后依次考虑每一门课。从没有被B神碾压的(n-k-1)位同学中，选出(n-r_i-k)位虽然没有全面被B神碾压，但是该门课得分小于等于B神的同学。枚举B神该门课的得分，这样就能根据排名情况求出其他人得分的方案数。写成式子就是：

[g_k={n-1choose k}prod_{i=1}^{m}{n-k-1choose n-r_i-k}sum_{j=1}^{h_i}j^{n-r_i}(h_i-j)^{r_i-1} ]

暴力按照这个式子模拟，求(g)的时间复杂度(O(mh))。总时间复杂度(O(nmh))。期望得(40)分。

考虑优化。不难发现瓶颈在于求(sum_{j=1}^{h_i}j^{n-r_i}(h_i-j)^{r_i-1})。而这一部分又和(k)无关，所以可与对每个(i)预处理出(s_i=sum_{j=1}^{h_i}j^{n-r_i}(h_i-j)^{r_i-1})。把式子中的((h_i-j)^{r_i-1})用二项式定理展开，再把(j)的次幂合并到一起。得到：

[egin{align} s_i=&sum_{j=1}^{h_i}j^{n-r_i}(h_i-j)^{r_i-1}\ =&sum_{j=1}^{h_i}j^{n-r_i}sum_{l=0}^{r_i-1}{r_i-1choose l}h^{r_i-1-l}(-j)^{l}\ =&sum_{j=1}^{h_i}j^{n-r_i}sum_{l=0}^{r_i-1}{r_i-1choose l}h^{r_i-1-l}(-1)^lj^l\ =&sum_{l=0}^{r_i-1}{r_i-1choose l}h^{r_i-1-l}(-1)^lsum_{j=1}^{h_i}j^{n-r_i+l} end{align} ]

于是问题进一步转化为求(sum_{j=1}^{h_i}j^{n-r_i+l})。发现(h_i)很大但(n-r_i+l)不大。所以这是一个典型的自然数幂求和问题。自然数幂求和的方法很多。例如，可以使用拉格朗日插值法。用这个方法，可以在(O(n))的时间里求出(sum_{j=1}^{h_i}j^{n-r_i+l})。于是求出单个(s_i)的时间复杂度优化为(O(n^2))，求出(s_{1dots m})的时间复杂度就是(O(mn^2))。后面用预处理好的(s_i)，套一个二项式反演，可以(O(nm))求出答案。故总时间复杂度为(O(mn^2))。

参考代码（在LOJ查看）：

//problem:LOJ2026
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

const int MAXN=110;
const int MOD=1e9+7;
inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
inline void add(int& x,int y){x=mod1(x+y);}
inline void sub(int& x,int y){x=mod2(x-y);}
inline int pow_mod(int x,int i){int y=1;while(i){if(i&1)y=(ll)y*x%MOD;x=(ll)x*x%MOD;i>>=1;}return y;}

int fac[MAXN+5],ifac[MAXN+5];
inline int comb(int n,int k){
	assert(k>=0);
	if(n<k)return 0;
	return (ll)fac[n]*ifac[k]%MOD*ifac[n-k]%MOD;
}
void facinit(int lim=MAXN){
	fac[0]=1;
	for(int i=1;i<=lim;++i)fac[i]=(ll)fac[i-1]*i%MOD;
	ifac[lim]=pow_mod(fac[lim],MOD-2);
	for(int i=lim-1;i>=0;--i)ifac[i]=(ll)ifac[i+1]*(i+1)%MOD;
}

int n,m,K,h[MAXN+5],rk[MAXN+5],s[MAXN+5];

int main() {
	facinit();
	cin>>n>>m>>K;
	for(int i=1;i<=m;++i)cin>>h[i];
	int mn=n;
	for(int i=1;i<=m;++i)cin>>rk[i],mn=min(mn,n-rk[i]);
	static int pw[MAXN+5][MAXN+5];
	for(int i=1;i<=MAXN;++i){
		for(int j=0;j<=MAXN;++j){
			pw[i][j]=pow_mod(i,j);
		}
	}
	for(int i=1;i<=m;++i){
		//预处理s[i]
		/*
		for(int j=1;j<=h[i];++j){
			add(s[i],(ll)pow_mod(j,n-rk[i])*pow_mod(h[i]-j,rk[i]-1)%MOD);
		}
		*/
		for(int l=0;l<=rk[i]-1;++l){
			static int pre[MAXN+5],suf[MAXN+5];
			int lim=n-rk[i]+l+2;
			pre[0]=1;
			suf[lim+1]=1;
			for(int j=1;j<=lim;++j)pre[j]=(ll)pre[j-1]*mod2(h[i]%MOD-j)%MOD;
			for(int j=lim;j>=1;--j)suf[j]=(ll)suf[j+1]*mod2(h[i]%MOD-j)%MOD;
			int sum=0;
			/*
			for(int j=1;j<=h[i];++j){
				add(sum,pow_mod(j,n-rk[i]+l));
			}
			*/
			for(int j=1,y=0;j<=lim;++j){
				//add(y,pow_mod(j,n-rk[i]+l));
				add(y,pw[j][n-rk[i]+l]);
				int tmp=(ll)y*pre[j-1]%MOD*suf[j+1]%MOD*ifac[j-1]%MOD*ifac[lim-j]%MOD;
				if((lim-j)&1)sub(sum,tmp);
				else add(sum,tmp);
			}
			int tmp=(ll)pow_mod(h[i],rk[i]-1-l)*comb(rk[i]-1,l)%MOD*sum%MOD;
			if(l&1)sub(s[i],tmp);
			else add(s[i],tmp);
		}
	}
	int ans=0;
	for(int k=K;k<=mn;++k){
		int f=comb(n-1,k);
		for(int i=1;i<=m;++i){
			f=(ll)f*s[i]%MOD*comb(n-k-1,n-rk[i]-k)%MOD;
		}
		if((k-K)%2==1)sub(ans,(ll)comb(k,K)*f%MOD);
		else add(ans,(ll)comb(k,K)*f%MOD);
	}
	cout<<ans<<endl;
	return 0;
}