jzoj 3206.【HNOI模拟题】Axonometric Projection dp+组合计数

Description
让我们来考虑一个单位立方体建成的模型。这个建筑的底是一个n m的单位正方形网格。在每个正方形上面，堆着若干个（可能是0）个单位立方体。每个立方体属于其中一个立方体堆。

给出了一个建筑的左视图和正视图。请计算有多少种建筑，符合给出的左视图和正视图。答案可能很大，只要返回它除以10^9 + 9的余数即可。

Input
第一行是整数n。第二行描述了建筑的左视图。第i个数表示了由上往下看时第i行最高的立方体堆的高度。第三行是整数m。第四行描述了建筑的正视图。第i个数表示了由上往下看时第i列最高的立方体堆的高度。

Output
仅一个数，表示答案。

Sample Input
2
1 1
2
1 1

Sample Output
7

Data Constraint
对于10%的数据，1<= n,m <= 5。
对于30%的数据，1 <= n,m <=10。
对于60%的数据，1 <= n,m <= 25。
对于100%的数据，1 <= n,m <= 50，所有出现的数不超过10^4。

分析：
题目可以看成有一个$n\ast m$的网格图，第$i$行的最大值为$a_i$，第$i$列的最大值为$b_i$的方案数。
显然交换两个行的限制$a_i$和$a_j$，方案数不变；列同理。
我们对$a_i$和$b_i$排序。
网格图中$c_{i,j}$取值显然是$[0,min(a_i,b_j)]$。
我们举一个例子，设$a=2,2,3$，$b=1,2,3$，那么
$c=\left(\begin{array}{ccc}1& 2& 2\\ 1& 2& 2\\ 1& 2& 3\end{array}\right)$
我们发现，对于同一个数字，会组成一个$L$型，或者是退化的$L$型。
每一个$L$型之间相互独立，答案就是他们各自方案的积。

我们把这个$L$型补成一个大矩形减一个小矩形的方式。设小矩形长宽为$(n,m)$，$L$型靠近边界两边为$(x,y)$，大矩形表示为$(x+n,y+m)$，我们使用四元组$(x,y,n,m)$表示一个$L$型。

对于上面的例子，$2$的第$1$行、第二行和第二列需要至少有一个$2$（上限）。
统计答案可以这样算，用总的方案减去一行没有上限的方案，再减去两行没有上限的方案……

对于计算一个$(x,y,n,m)$的$L$型的答案，设上限为$h$，$f[i][j]$表示$(i,j,n,m)$的$L$型的答案，那么
$f[i][j]=(h+1{)}^{(i+n)\ast (j+m)-n\ast m}-{\sum }_{a=0}^i{\sum }_{b=0}^{j}f[a][b]\ast (\genfrac{}{}{0ex}{}{i}{a})\ast (\genfrac{}{}{0ex}{}{j}{b})\ast h^{(i+n)\ast (j+m)-(a+n)\ast (b+m)}$
后面表示把一部分强行钦定在$[0,h-1]$范围内。
可以预处理次方。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#define LL long long

const int maxp=1e4+7;
const int maxn=107;
const LL mod=1e9+9;

using namespace std;

int n,m,x,cn,cm,maxa,maxb;
int a[maxp],b[maxp];
LL f[maxn][maxn],C[maxn][maxn],p[maxn*maxn],ans;

LL power(LL x,LL y)
{
	if (y==0) return 1;
	if (y==1) return x;
	LL c=power(x,y/2);
	c=(c*c)%mod;
	if (y&1) c=(c*x)%mod;
	return c;
}

LL dp(int h,int x,int y,int n,int m)
{
	memset(f,0,sizeof(f));
	f[0][0]=1;
	for (int i=0;i<=x;i++)
	{
		for (int j=0;j<=y;j++)
		{
			LL d;
			d=(i+n)*(j+m)-n*m;
			f[i][j]=power(h,d);
			for (int a=0;a<=i;a++)
			{
				for (int b=0;b<=j;b++)
				{
					if ((a==i) && (b==j)) continue;
					d=(i+n)*(j+m)-(a+n)*(b+m);
					f[i][j]=(f[i][j]+mod-f[a][b]*C[i][a]%mod*C[j][b]%mod*p[d]%mod)%mod;
				}
			}
		}
	}
	return f[x][y];
}

int main()
{
	scanf("%d",&n);
	for (int i=1;i<=n;i++)
	{
		scanf("%d",&x);
		a[x]++;
		maxa=max(maxa,x);
	}
	scanf("%d",&m);
	for (int i=1;i<=m;i++)
	{
		scanf("%d",&x);
		b[x]++;
		maxb=max(maxb,x);
	}
	if (maxa!=maxb)
	{
		printf("0\n");
		return 0;
	}	
	C[0][0]=1;
	for (int i=1;i<=50;i++)
	{
		C[i][0]=1;
		for (int j=1;j<=i;j++) C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;
	}
	cn=cm=0;
	ans=1;
	for (int i=1e4;i>=0;i--)
	{
		if ((!a[i]) && (!b[i])) continue;
		p[0]=1;
		for (int j=1;j<=n*m;j++) p[j]=(p[j-1]*(LL)i)%mod;
		ans=(ans*dp(i+1,a[i],b[i],cn,cm))%mod;
		cn+=a[i],cm+=b[i];
	}
	printf("%lld\n",ans);
}