[CSP-S模拟测试]:Silhouette（数学）

Posted 2022-12-25 wzc521

题目描述

　　有一个$n\\times n$的网格，在每个格子上堆叠了一些边长为$1$的立方体。
　　现在给出这个三维几何体的正视图和左视图，求有多少种与之符合的堆叠立方体的方案。两种方案被认为是不同的，当且仅当某个格子上立方体的数量不同。
　　输出答案对$10^9+7$取模的结果。

---

输入格式

　　从文件$silhouette.in$中读入数据。
　　第一行一个整数$n$。
　　第二行$n$个整数，第$i$个表示正视图中从左到右第$i$个位置的高度$A_i$。
　　第三行$n$个整数，第$i$个表示左视图中从左到右第$i$个位置的高度$B_i$。

---

输出格式

　　输出到文件$silhouette.out$中。
　　输出一行表示答案。

---

样例

样例输入1：

2
1 2
2 1

样例输出1：

5

样例输入2：

3
3 1 3
2 3 2

样例输出2：

175

样例输入3：

3
1 1 1
3 3 3

样例输出3：

0

---

数据范围与提示

样例$1$解释：

正视图和左视图：

    _    _
 _|_|  |_|_
|_|_|  |_|_|

如果用$2\\times 2$的矩阵来表示每个格子上堆叠的立方体个数，则五种方案可以表示为：

1 2  0 2  1 2  0 2  1 2
1 1  1 0  1 0  1 1  0 1

数据范围：

对于所有数据，有$1\\leqslant n\\leqslant 10^5,1\\leqslant A_i,B_i\\leqslant 10^9$。
$\\bullet Subtask1(3\\%)$，$n=1$。
$\\bullet Subtask2(14\\%)$，$n\\leqslant 3,A_i,B_i\\leqslant 4$。
$\\bullet Subtask3(17\\%)$，$n\\leqslant 16$。
$\\bullet Subtask4(24\\%)$，$n\\leqslant 100$。
$\\bullet SUbtask5(13\\%)$，$n\\leqslant 3,000$。
$\\bullet Subtask6(18\\%)$，$A_i,B_i$分别构成了一个$1$至$n$的排列。
$\\bullet Subtask7(11\\%)$，没有特殊的约束。

---

题解

显然对于这道题，我们要求以下方程的解的个数：

$$\\forall i\\in [1,n],\\max \\limitsj=1^n x_i,j=A_i,\\max \\limits_j=1^n x_j,i=B_i$$

我们可以现将$A,B$排序，因为这样对结果没有影响，简单证明一下：

先来考虑列，对于每一个$A_i$，无论哪一列在前，哪一列在后，最终所有列在第$i$行的最大值都需要是$A_i$；行同理。

如果最大的$A_i$不等于最大的$B_i$，那么无解，否则一定有解，不再证明（主要是没有人问过）。

然后我们从大到小枚举$A,B$中每一个值$S$，对于排好序后每一个$S$，它们会是这样的：

类似分成了一层一层的。

然后我们枚举每一层。

先说$S_1$这一层，对于$S_1$有一个特殊性质，$S_1$是所有$S$中最大的，先设$S_1$所涉及到的这个矩形的长宽分别为$a,b$，如下图：

我们在设$f[i]$表示$a$行中，至少有$i$行一定不合法的方案数，这里说的不合法是指高度没有到达$S$，也就是说不能对投影做贡献；之所以我们要设为至少，是因为这样每两行之间就可以不互相影响了。

给出状态转移方程：

$$f[i]=C_a^i\\times (S^i\\times ((S+1)^a-i-S^a-i))^b$$

来解释一下状态转移方程：

首先，组合数必不可少，因为我们要从$a$行里选出$i$行；再解释后面的$b$次方，因为有$b$列，这$b$列互不影响（前面有解释）；再来解释$S^i$，因为有$i$个位置一定不合法，这$i$个位置可以选$0\\sim S-1$这$S$个高度；最后来理解$(S+1)^a-i-S^a-i$，因为还剩下$a-i$个位置，每一个位置选多高都可以，所以有$(S+1)^a-i$，但是我们要保证这一列一定合法，所以还要减去$S^a-i$。

因为我们$f[i]$的定义为至少有$i$行不合法的方案数，所以我们还需要求得恰好有$i$行不合法的方案数，而我们只需要知道恰好有$0$行不合法的方案数（设为$res$），也就是都合法的方案数，考虑容斥，则有：

$$res=\\sum \\limits_i=0^a(-1)^i\\times f[i]$$

那么，我们在来考虑一般情况，也就是上图中$S_2$以后的所有$S$，因为每当我们处理完一个$S$之后，下一个区域的形状只有$L$行或矩形，如下图（红色为上一次处理的区域，紫色为这一次处理的区域）：

而无论对于哪种情况，我们都按如下划分方式将其划分为两部分：

对于矩形的两种情况，无非就是没有了那一部分，可以认为举行是一种特殊的$L$行。

为方便，我们不妨将$L$行区域按如下图方式标号$a,b,c,d$：

那么，现在我给出状态转移方程：

$$f[i]=C_a^i\\times (S^i\\times ((S+1)^a+c-i-S^a+c-i))^b\\times (S^i\\times (S+1)^a-i)^d$$

现在来解释这个更复杂的式子，先明确一下，因为上图中红色区域已经处理完毕，所以对于蓝色区域，其已经满足行，而没有满足列，对于绿色区域同理。

先来解释式子中的组合数$C_a^i$，你可能会存在疑问，为什么有$a+c$行，而不是$C_a+c^i$；那是因为上面我也明确过了，对于$c$行，行已经满足，我们不能让它不满足，而这$i$行的不满足我们只能从$a$行出；所以是$C_a^i$，而不是$C_a+c^i$。

再来解释$(S^i\\times ((S+1)^a+c-i-S^a+c-i))^b$，对于$b$次方和$S^i$，上面对于$S$是最大的值的情况已经做了解释，在此就不做过多的赘述；直接解释$(S+1)^a+c-i-S^a+c-i$，这个也很好理解，与上面情况类似，有$i$行不合法，那么还有$a+c-i$个位置可以合法，而我们要保证这一列一定合法，所以还要减去不合法的情况。

最后来看后面新填进来的这一部分，应该能看的出来，这部分是为了处理矩形$a\\times d$的，也就是区域$2$，式子与前面大体类似，不再赘述，只来思考这样一个问题：网上有另一篇题解（在我写之前可能也只有那一篇，而我就是按照那一篇的思路调出来的），他的后面这部分是$S^i\\times (S+1)^a-i-S^a-i)^d$，但是我的式子里却没有$-S^a-i$，他的式子的漏洞就在这里，而为什么不用这部分呢？因为我们减去$S^a-i$是为了让那一列合法，还是开始我明确的那个问题，因为我们在处理红色区域的时候已经保证了其合法，而红色区域的$S$要比矩形$a\\times d$（绿色区域）的大，也就是高，所以无论如何其正视图都已经不变化了，变化的只有左视图，而我们需要做的就是让左视图合法，不用考虑列合不合法的问题，所以不用减去$S^a-i$。

对于统计合法方案数的容斥，其式子亦是：

$$res=\\sum \\limits_i=0^a (-1)^i\\times f[i]$$

对于$c=0$和$d=0$的情况，我们将其带入上式会发现对结果并无影响，于是我们可以将特殊情况的式子和普通情况的式子合并成一个。

这样我们从达到小不断枚举每一个$S$即可求出所有合法方案数。

时间复杂度：$\\Theta(n\\log n)$。

期望得分：$100$分。

实际得分：$100$分。

---

代码时刻

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int n;
int A[100001],B[100001],S[200001];
long long jc[100001],inv[100001];
long long ans=1;
long long qpow(long long x,long long y)

	long long res=1;
	while(y)
	
		if(y&1)res=res*x%mod;
		x=x*x%mod;
		y>>=1;
	
	return res;

void pre_work()

	jc[0]=1;
	for(long long i=1;i<=100000;i++)jc[i]=jc[i-1]*i%mod;
	inv[100000]=qpow(jc[100000],mod-2);
	for(long long i=100000;i>0;i--)inv[i-1]=inv[i]*i%mod;

long long get_C(int x,int y)return jc[x]*inv[y]%mod*inv[x-y]%mod;
long long lucas(int x,int y)

	if(!y)return 1;
	return get_C(x%mod,y%mod)*lucas(x/mod,y/mod)%mod;

long long solve(int a,int b,int c,int d,int s)

	long long res=0;
	for(int i=0;i<=a;i++)
	
		long long now=lucas(a,i)*qpow(qpow(s,i)*((qpow(s+1,a+c-i)-qpow(s,a+c-i)+mod)%mod)%mod,b)%mod*qpow(qpow(s,i)*qpow(s+1,a-i)%mod,d)%mod;
		if(i&1)res=(res-now+mod)%mod;
		else res=(res+now)%mod;
	
	return res;

int main()

	pre_work();
	scanf("%d",&n);
	for(int i=1;i<=n;i++)scanf("%d",&A[i]);S[i]=A[i];
	for(int i=1;i<=n;i++)scanf("%d",&B[i]);S[n+i]=B[i];
	sort(A+1,A+n+1);
	sort(B+1,B+n+1);
	if(A[n]!=B[n])puts("0");return 0;
	sort(S+1,S+2*n+1);
	S[0]=unique(S+1,S+2*n+1)-S-1;
	int prea=n+1,preb=n+1,nowa=n,nowb=n;
	for(int i=S[0];i;i--)
	
		while(A[nowa-1]==S[i]&&nowa-1)nowa--;
		while(B[nowb-1]==S[i]&&nowb-1)nowb--;
		ans=ans*solve(prea-nowa,preb-nowb,n-prea+1,n-preb+1,S[i])%mod;
		prea=nowa;preb=nowb;
	
	printf("%lld",ans);
	return 0;

---

rp++