状态压缩递推(States Compressing Recursion,SCR)

 我们暂时避开状态压缩的定义，先来看一个小小的例题。

【引例】

    在 n*n(n≤20)的方格棋盘上放置n 个车(可以攻击所在行、列)，求使它们不能互相攻击的方案总数。

 【分析】

    这个题目之所以是作为引例而不是例题，是因为它实在是个非常简单的组合学问题：我们一行一行放置，则第一行有n 种选择，第二行n-1，……，最后一行只有 1 种选择，根据乘法原理，答案就是n! 。这里既然以它作为状态压缩的引例，当然不会是为了介绍组合数学。我们下面来看另外一种解法：状态压缩递推(States Compressing Recursion,SCR)。

    我们仍然一行一行放置。取棋子的放置情况作为状态，某一列如果已经放置棋子则为 1，否则为0。这样，一个状态就可以用一个最多20 位的二进制数2表示。例如n=5,第 1、3、4 列已经放置，则这个状态可以表示为01101(从右到左)。设f[s]为达到状态s 的方案数，则可以尝试建立f 的递推关系。

    考虑n=5,s=01101 。这个状态是怎么得到的呢？因为我们是一行一行放置的，所以当达到 s 时已经放到了第三行。又因为一行能且仅能放置一个车，所以我们知道状态s 一定来自：

①前两行在第 3、4  列放置了棋子(不考虑顺序，下同)，第三行在第 1 列放置；

②前两行在第 1、4 列放置了棋子，第三行在第3 列放置；

③前两行在第 1、3 列放置了棋子，第三行在第4 列放置。

这三种情况互不相交，且只可能有这三种情况。根据加法原理，f[s]应该等于这三种情况的和。写成递推式就是：

               f[01101]=f[01100]+f[01001]+f[00101]

根据上面的讨论思路推广之，得到引例的解决办法：                     

[cpp] view plaincopy

1. f[0]=1   
2. f[s]=∑f[s^2i]  

其中s∈[0…01,1…11]，s的右起第i+1位为1(3)。 

   反思这个算法，其正确性毋庸置疑(可以和n!对比验证)。但是算法的时间复杂度为O(n2^n )，空间复杂度O(2^n )，是个指数级的算法，比循环计算n!差了好多，它有什么优势？较大的推广空间。(4) 

Sample Problems 

【例1】(5)

    在n*n(n≤20)的方格棋盘上放置n 个车，某些格子不能放，求使它们不能互相攻击的方案总数。

=====================================================================================

1 本文所有例题的C++代码均可以在附件中找到。

2 方便、整齐起见，本文中不在数的后面标明进制。

3 考虑上节介绍的位运算的特殊应用，可以更精巧地实现。

4 还有一个很…的用处，即对新手说：“来看看我这个计算n!的程序，连这都看不懂就别OI 了~”

5 题目来源：经典组合学问题。

=====================================================================================

 【分析】

    对于这个题目，如果组合数学学得不够扎实，你是否还能一眼看出解法？应该很难。对于这个题目，确实存在数学方法(容斥原理)，但因为和引例同样的理由，这里不再赘述。

    联系引例的思路，发现我们并不需要对算法进行太大的改变。引例的算法是在枚举当前行(即s 中1 的个数，设为r)的放置位置(即枚举每个 1)，而对于例 1，第r 行可能存在无法放置的格子，怎么解决这个问题呢？枚举 1 的时候判断一下嘛！事实的确是这样，枚举 1 的时候判断一下是否是不允许放置的格子即可。

    但是对于n=20，O(n2^n )的复杂度已经不允许我们再进行多余的判断。所以实现这个算法时应该应用一些技巧。对于第r 行，我们用 a[r]表示不允许放置的情况，如果某一位不允许放置则为 1，否则为0，这可以在读入数据阶段完成。运算时，对于状态 s，用tmps=s^a[r]来代替 s 进行枚举，即不枚举 s 中的 1 转而枚举 tmps  中的 1。因为tmps 保证了无法放置的位为 0，这样就可以不用多余的判断来实现算法，代码中只增加了计算a 数组和r 的部分，而时间复杂度没有太大变化。

    这样，我们直接套用引例的算法就使得看上去更难的例 1 得到了解决。你可能会说，这题用容斥原理更快。没错，的确是这样。但是，容斥原理在这题上只有当棋盘为正方形、放入的棋子个数为n、且棋盘上禁止放置的格子较少时才有简单的形式和较快的速度。如果再对例 1 进行推广，要在 m*n  的棋盘上放置 k个车，那么容斥原理是无能为力的，而 SCR 算法只要进行很少的改变就可以解决问题(1) 。这也体现出了引例中给出的算法具有很大的扩展潜力。

 

    棋盘模型是状态压缩最好的展示舞台之一。下面再看几个和棋盘有关的题目。

 

 【例2】(2)

    给出一个n*m 的棋盘(n、m≤80,n*m≤80)，要在棋盘上放k(k≤20)个棋子，使得任意两个棋子不相邻。每次试验随机分配一种方案，求第一次出现合法方案时试验的期望次数，答案用既约分数表示。

 

 【分析】

    显然，本题中的期望次数应该为出现合法方案的概率的倒数，则问题转化为求出现合法方案的概率。而概率=合法方案数/方案总数 , 方案总数显然为C(n*m,k)，则问题转化为求合法方案数。整理一下，现在的问题是：在 n*m  的棋盘上放k 个棋子，求使得任意两个棋子不相邻的放置方案数。

 

    这个题目的状态压缩模型是比较隐蔽的。观察题目给出的规模，n、m≤80，这个规模要想用 SC 是困难的，若同样用上例的状态表示方法(放则为 1，不放为0)，2^80  无论在时间还是在空间上都无法承受。然而我们还看到 n*m≤80，这种给出数据规模的方法是不多见的，有什么玄机呢？能把状态数控制在可以承受的范围吗？稍微一思考，我们可以发现：9*9=81>80，即如果n,m 都大于等于 9，将不再满足n*m≤80 这一条件。所以，我们有n 或m 小于等于 8，而2^8 是可以承受的。我们假设 m≤n(否则交换，由对称性知结果不变)n  是行数 m  是列数，则每行的状态可以用 m 位的二进制数表示。但是本题和例 1 又有不同：例 1 每行每列都只能放置一个棋子，而本题却只限制每行每列的棋子不相邻。但是，上例中枚举当前行的放置方案的做法依然可行。我们用数组 s[1..num] 保存一行中所有的num 个放置方案，则s 数组可以在预处理过程中用DFS 求出，同时用c[i]保存第i 个状态中 1 的个数以避免重复计算。开始设计状态。如注释一所说，维数需要增加，原因在于并不是每一行只放一个棋子，也不是每一行都要求有棋子，原先的表示方法已经无法完整表达一个状态。我们用 f[i][j][k]表示第 i 行的状态为s[j]且前i 行已经放置了k 个棋子(2) 的方案数。沿用枚举当前行方案的做法，只要当前行的方案和上一行的方案不冲突即可，“微观”地讲，即s[snum[i]]和s[snum[i-1]]没有同为 1 的位，其中snum[x]表示第x 行的状态的编号。然而，虽然我们枚举了第 i 行的放置方案，但却不知道其上一行(i-1)的方案。为了解决这个问题，我们不得不连第i-1 的状态一起枚举，则可以写出递推式：

[cpp] view plaincopy

1. f[0][1][0]=1;  
2. f[i][j][k]=∑f[i-1][p][k-c[j]]  

其中s[1]=0，即在当前行不放置棋子；j和p是需要枚举的两个状态编号，且要求s[j]与s[p]不冲突，即s[j]&s[p]=0。(3) 当然，实现上仍有少许优化空间，例如第i行只和第i-1行有关，可以用滚动数组节省空间。 

    有了合法方案数，剩下的问题就不是很困难了，需要注意的就只有C(n*m,k)可能超出 64 位整数范围的问题，这可以通过边计算边用 GCD 约分来解决，具体可以参考附件中的代码。这个算法时间复杂度 O(n*pn*num)，空间复杂度(滚动数组)O(pn*num)，对于题目给定的规模是可以很快出解的。 

=====================================================================================

1 如果这样改编，则状态的维数要增加，如有疑问可以参考下面的几个例子，这里不再赘述。

2 题目来源：经典问题改编。

=====================================================================================

    通过上文的例题，读者应该已经对状态压缩有了一些感性的认识。下面这个题目可以作为练习。 

 

 【例3】(4)

    在 n*n(n≤10)的棋盘上放k 个国王(可攻击相邻的 8 个格子)，求使它们无法互相攻击的方案数。

 

 【分析】

    其实有了前面几个例子的分析，这个题目应该是可以独立解决的。不过既然确实有疑问，那我们就来分析一下。

=====================================================================================

1 运用简单的组合数学知识可以求出：在格数为m 的一行上放棋子且相邻两个棋子中间的空格不能少于d的num=g[m+d+1]，其中g[i]=1 (i=1..d+1);  g[j]=g[j-d-1]+g[j-1] (j>d)。对于本题，num=144。此式在下文也有应用。

2 为了避免歧义，下文中用pn 代替原题中的k 。

3 请读者停下来仔细揣摩这个递推式，否则可能影响对本文后面内容的理解！

4 题目来源：SGU.223  《Little Kings》

=====================================================================================

 首先，你应该能想到将一行的状态 DFS  出来(如果不能，请返回重新阅读，谢谢)，仍然设为s[1..num],同时仍然设有数组c[1..num]记录状态对应的 1 的个数。和例2 相同，仍然以f[i][j][k]表示第i 行状态为 s[j]，且前i 行已经放置了k 个棋子的方案数。递推式仍然可以写作：

[cpp] view plaincopy

1. f[0][1][0]=1;  
2. f[i][j][k]=∑f[i-1][p][k-c[j]]  

其中仍要求s[j]和s[p]不冲突。 

    可是问题出来了：这题不但要求不能行、列相邻，甚至不能对角线相邻！s[j]、s[p]不冲突怎么“微观地”表示呢？其实，稍微思考便可以得出方法：用 s[p]分别和 s[j]、s[j]*2、s[j]/2 进行冲突判断即可，原理很显然。解决掉这唯一的问题，接下来的工作就没有什么难度了。算法复杂度同例2。 

 

    下一个例题是状态压缩棋盘模型的经典题目，希望解决这个经典的题目能够增长你的自信。 

 

 【例4】(1)

    给出一个 n*m(n≤100,m≤10)的棋盘，一些格子不能放置棋子。求最多能在棋盘上放置多少个棋子，使得每一行每一列的任两个棋子间至少有两个空格。

 【分析】 (2)

    显然，你应该已经有 DFS 搜出一行可能状态的意识了(否则请重新阅读之前的内容3 遍，谢谢)，依然设为s[1..num],依旧有c[1..num]保存 s 中1 的个数，依照例 1 的预处理搞定不能放置棋子的格子。

    问题是，这个题目的状态怎么选？继续像例 2、3  那样似乎不行，原因在于棋子的攻击范围加大了。但是我们照葫芦画瓢：例 2、3  的攻击范围只有一格，所以我们的状态中只需要有当前行的状态即可进行递推，而本题攻击范围是两格，因此增加一维来表示上一行的状态。用f[i][j][k]表示第i 行状态为 s[j]、第i-1行状态为 s[k]时前i 行至多能放置的棋子数，则状态转移方程很容易写出：

[cpp] view plaincopy

1. f[i][j][k]=max{f[i-1][k][l]}+c[j]  

其中要求 s[j],s[k],s[l]互不冲突。

    因为棋子攻击范围为两格，可以直观地想象到num 不会很大。的确，由例2中得到的num 的计算式并代入d=2、m=10，得到num=60 。显然算法时间复杂度为O(n*num^3)，空间复杂度(滚动数组)O(num^2)。此算法还有优化空间。我们分别枚举了三行的状态，还需要对这三个状态进行是否冲突的判断，这势必会重复枚举到一些冲突的状态组合。我们可以在计算出 s[1..num]后算出哪些状态可以分别作为两行的状态，这样在DP 时就不需要进行盲目的枚举。这样修改后的算法理论上比上述算法更优，但因为 num  本身很小，所以这样修改没有显著地减少运，但由于位行时间。值得一提的是，本题笔者的算法虽然在理论上并不是最优 (3)运算的使用，截至2 月9  日，笔者的程序在PKU OJ 上长度最短，速度第二快。

    这个题目是国内比赛中较早出现的状态压缩题。它告诉我们状态压缩不仅可以像前几个例题那样求方案数，而且可以求最优方案，即状态压缩思想既可以应用到递推上(SCR)，又可以应用到DP 上(SCDP)，更说明其有广泛的应用空间。

=====================================================================================

1 题目来源：NOI2001   《炮兵阵地》；PKU.1185

2 读者应该独立思考一个小时后再看分析，它值得这样。

3 有种应用三进制的方法理论上可能更优，但因为手工转换进制效率远不如系统的位运算高，且无法使用位运算判断可行性，程序效率并不比文中的高，故不再赘述那种方法。下文的一些题目将用到多进制。

===================================================================================== 

    

    看了这么多棋盘模型应用状态压缩的实例，你可能会有疑问，难道状态压缩只在棋盘上放棋子的题目中有用？不是的。我们暂时转移视线，来看看状态压缩在其他地方的应用——覆盖模型。

 

 【例5】（1）

 

    给出n*m   (1≤n、m≤11)的方格棋盘，用1*2 的长方形骨牌不重叠地覆盖这个棋盘，求覆盖满的方案数。

 

 【分析】

    这也是个经典的组合数学问题：多米诺骨牌完美覆盖问题(或所谓二聚物问题)。                               

    显然，如果 n、m  都是奇数则无解( 由棋盘面积的奇偶性知)，否则必然有至少一个解(很容易构造出)，所以假设n、m 至少有一个偶数，且m≤n(否则交换)。我们依然像前面的例题一样把每行的放置方案 DFS出来，逐行计算。用 f[i][s]表示把前i-1 行覆盖满、第i 行覆盖状态为 s 的覆盖方案数。因为在第i 行上放置的骨牌最多也只能影响到第i-1 行，则容易得递推式：

[cpp] view plaincopy

1. f[0][1…11]=1  
2. f[i][s1]=∑f[i-1][s2]  

其中(s1,s2)整体作为一个放置方案，可以把所有方案 DFS 预处理出来。下面讨论一下本题的一些细节。 

    首先讨论DFS 的一些细节。对于当前行每一个位置，我们有3 种放置方法：①竖直覆盖，占据当前格和上一行同一列的格；②水平覆盖，占据当前格和该行下一格；③不放置骨牌，直接空格。如何根据这些枚举出每个(s1,s2)呢？下面介绍两种方法：

第一种：

        DFS 共 5 个参数，分别为：p(当前列号)，s1、s2(当前行和上一行的覆盖

    情况)，b1、b2(上一列的放置对当前列上下两行的影响，影响为 1 否则为0)。

    初始时 s1=s2=b1=b2=0。①p=p+1，s1=s1*2+1，s2=s2*2^2  ，b1=b2=0；②p=p+1，

    s1=s1*2+1，s2=s2*2+1，b1=1，b2=0；③p=p+1，s1=s1*2，s2=s2*2+1，b1=b2=0 。

    当p 移出边界且b1=b2=0 时记录此方案。

 

第二种：

        观察第一种方法，发现b2 始终为0，知这种方法有一定的冗余。换个更

    自然的方法，去掉参数b1、b2 。①p=p+1，s1=s1*2+1，s2=s2*2；②p=p+2，

    s1=s1*4+3，s2=s2*4+3；③p=p+1，s1=s1*2，s2=s2*2+1。当p 移出边界时记

    录此方案。这样，我们通过改变p 的移动距离成功简化了DFS 过程，而且这

    种方法更加自然。

 

=====================================================================================

1 题目来源：经典问题；PKU.2411      《Mondriaan's Dream》

2 注意！第i 行的放置方案用到第i-1 行的某格时，s2 中该格应为0 ！

=====================================================================================

    DFS 过程有了，实现方法却还有值得讨论的地方。前面的例题中，我们为什么总是把放置方案 DFS  预处理保存起来？是因为不合法的状态太多，每次都重新DFS 太浪费时间。然而回到这个题目，特别是当采用第二种时，我们的 DFS过程中甚至只有一个判断(递归边界)，说明根本没有多少不合法的方案，也就没有必要把所有方案保存下来，对于每行都重新 DFS  即可，这不会增加运行时间却可以节省一些内存。

 

    这个算法时间复杂度为多少呢？因为 DFS  时以两行为对象，每行 2^m ，共进行 n 次DFS，所以是O(n*4^m  )？根据“O”的上界意义来看并没有错，但这个界并不十分精确，也可能会使人误以为本算法无法通过 1≤n、m≤11 的测试数据,而实际上本算法可以瞬间给出 m=10,n=11时的解。为了计算精确的复杂度，必须先算出DFS 得到的方案数。

    考虑当前行的放置情况。如果每格只有①③两个选择，则应该有 2m种放置方案；如果每格有①②③这 3 个选择，且②中 p 只移动一格，则应该有 3^m 种放置方案。然而现在的事实是：每格有①②③这3 个选择，但②中p 移动2 格，所以可以知道方案数应该在2^m 和 3^m 之间。考虑第 i 列，则其必然是：第 i-1 列采用①③达到；第i-2 列采用②达到。设h[i]表示前i 列的方案数，则得到h[i]的递推式：

                             

[cpp] view plaincopy

1. h[0]=1,h[1]=2  
2. h[i]=2*h[i-1]+h[i-2]  

应用组合数学方法(1) 求得其通项公式h[m]=( (2 + sqrt(2))*(1 + sqrt(2)^m )/4 + (2 – sqrt(2))*(1 – sqrt(2))/4 )注意到式子的第二项是多个绝对值小于 1 的数的乘积，其对整个h[m]的影响甚小，故略去，得到方案数h[m]≈0.85*2.414^m ，符合2^m<h[m]<3^m  的预想。

 

    因为总共进行了 n 次 DFS，每次复杂度为 O(h[m])，所以算法总时间复杂度为O(n*h[m])=O(n*0.85*2.414m)，对m=10,n=11 不超时也就不足为奇了。应用滚动数组，空间复杂度为 O(2^m) 

 

    对于本题，我们已经有了公式和 SCR 两种算法。公式对于 m*n 不是很大的情况有效，SCR 算法在竞赛中记不住公式时对小的m、n 有效。如果棋盘规模为n*m(m≤10,n≤2^31-1)，则公式和 SCR 都会严重超时。有没有一个算法能在 1 分钟内解决问题呢(2) ？答案是肯定的，它仍然用到SC 思想。

 

=====================================================================================

1 特征方程等方法。这里不再赘述。

2 对于如此大的数据，高精度将成为瓶颈，故这里不考虑高精度的复杂度；我没有找到能 1s 出解的算法，

故延长时限到 1min。

=====================================================================================

    此算法中应用到一个结论：给出一个图的邻接矩阵 G(允许有自环，两点间允许有多条路径，此时G[i][j]表示i 到j  的边的条数)，则从某点a 走 k 步到某点b  的路径数为G^k  [a][b](1) 。本结论实际上是通过递推得到的，简单证明如下：从i走 k 步到j ，必然是从i 走 k-1 步到 t，然后从t 走 1 步到j ，根据加法原理，即G[k][i][j]=∑G[k-1][i][t]*G[t][j](2)  。是否感到这个式子很眼熟？没错，它和矩阵乘法 一模一样，即：G[k]=G[k-1]*G。因为矩阵乘法满足结合律，又由G[1]=G，所以我们得到结果：G[k]=G^k 。

    下面介绍这个算法。考虑一个有 2^m个顶点的图，每个顶点表示一行的覆盖状态，即 SCR 算法中的 s1 或 s2。如果(s1,s2)为一个放置方案，则在 s2 和 s1 之间连一条(有向)边，则我们通过DFS 一次可以得到一个邻接矩阵G。仍然按照逐行放置的思想来考虑，则要求我们每行选择一个覆盖状态，且相邻两行的覆盖状态(s1,s2)应为一个放置方案，一共有n 行，则要求选择n 个状态，在图中考虑，则要求我们从初始(第 0 行)顶点(1...111)n 步走到(1…111)，因为图的邻接矩阵是DFS 出来的，每条边都对应一个放置方案，所以可以保证走的每条边都合法。因此，我们要求的就是顶点(1…111)走 n步到达(1…111)的路径条数。由上面的结论知，本题的答案就是G^n  [1…111][1…111]。

    现在的问题是，如何计算 G  的n  次幂？连续 O(n)次矩阵乘法吗？不可取。矩阵的规模是 2^m *2^m  ，一次普通矩阵乘法要O((2^m )^3 )=O(8^m)，O(n)次就是O(n*8^m )， 比SCR 算法还要差得多。其实我们可以借用二分的思想。如果要计算3^8 的值，你会怎么算呢？直接累乘将需要进行7 次乘法。一种较简单的方法是：3*3=3  ，3^2 *3^2 =3^4  ，3^4 *3^4 =3  ，只进行了 3 次乘法，效率高了许多。因为矩阵乘法满足结合律，所以可以用同样的思路进行优化。这种思路用递归来实现是非常自然的，然而，本题的矩阵中可能有 2^10*2^10=2^20=1048576个元素，如果用(未经优化的)递归来实现，将可能出现堆栈溢出。不过庆幸的是我们可以非递归实现。用bin[]保存 n  的二进制的每一位，从最高位、矩阵G 开始，如果 bin[当前位]为 0，则把上一位得到的矩阵平方；如果为 1，则平方后再乘以G。这种方法的时间复杂度容易算出：O(logn)次矩阵乘法，每次O(8^m)，共O(8^m*logn)。

    这样对于m≤7就可以很快出解了。但对于m=n=8，上述算法都需要1s才能出解，无法令人满意。此算法还有优化空间。

    我们的矩阵规模高达 2^m*2^m  =4^m ，但是其中有用的(非0 的)有多少个呢？根据介绍 SCR 算法时得到的h[m]计算式，G 中有4^m-h[m]=4^m-0.85*2.414^m 个0，对于m=8，可以算出G 中98.5%的元素都是0，这是一个非常非常稀疏的矩阵，使用三次方的矩阵乘法有点大材小用。我们改变矩阵的存储结构，即第p 行第q 列的值为value 的元素可以用一个三元组(p,q,value)来表示，采用一个线性表依行列顺序来存储这些非0 元素。怎样对这样的矩阵进行乘法呢？观察矩阵乘法的计算式c[i][j] ∑a[i][k]*b[k][j] ，当a[i][k]或者b[k][j]为0 时，结果为0，对结果没有影响，完全可以略去这种没有意义的运算。则得到计算稀疏矩阵乘法的算法：枚举a 中的非0 元素，设为(p,q,v1)，在b  中寻找所有行号为q 的非0 元素(q,r,v2)，并把 v1*v2  的值累加到 c[p][r](4) 中。这个算法多次用到一个操作：找出所有行号为 q  的元素，则可以给矩阵附加一个数组hp[q]，表示线性表中第一个行号为 q的元素的位置，若不存在则hp[q]=0 。算出二维数组c 之后再对其进行压缩存储即可。此矩阵乘法的时间复杂度为 O( a.not0*b.not0/2^m + 4^m ) ，在最坏情况下，a.not0=b.not0=4^m ，算法的复杂度为O(8^m)，和经典算法相同。因为矩阵非常稀疏，算法复杂度近似为 O(4^m )(1)。考虑整个算法的时间复杂度：O(logn)次矩阵乘法， 每次 O(4^m)，则总时间复杂度O(logn*4^m)，对于m≤9 也可以很快出解了，对于m=10，n=2147483647，此算法在笔者机器上(Pm 1.6G，512M)运行时间少于20s。虽然仍然不够理想，但已经不再超时数小时。此算法空间复杂度为O(max_not0+4^m)，对于m=10，max_not0 小于 190000。

=====================================================================================

1 上标表示乘幂。

2 因为0 乘任何数都为0，所以G[k-1][i][t]或G[t][j]是否为0 对最后的和没有影响。

3 这里假设大整数间的乘法和小整数间的乘法耗费同样的时间，实际上并非这样。但对于大数据，我们只关心结果mod 某个常数的值，所以可以进行这样的假设。

4 此时c 是个2^m*2^m  的二维数组

=====================================================================================

 

    以上给出了公式、SCR、矩阵乘方这3 个算法，分别适用于不同的情况，本题基本解决。

 

    读者应该已经注意到，覆盖模型和棋盘模型有很多共同点，譬如都是在矩形某些位置放入棋子(或某种形状的骨牌)来求方案数(如上例)或最优解(下面将会给出几个例题)。但不难看出，覆盖模型和棋盘模型又有着很大的不同：棋盘模型中，只要棋子所在的位置不被别的棋子攻击到即可，而覆盖模型中，棋子的攻击范围也不可以重叠。所以简单来说，覆盖模型就是攻击范围也不能重叠的棋盘模型。下面再给出一个与上例类似的覆盖模型的例题以加深印象。

 

 【例6】(2)

    给出n*m   (1≤n、m≤9)的方格棋盘，用 1*2 的矩形的骨牌和L 形的(2*2 的去掉一个角)骨牌不重叠地覆盖，求覆盖满的方案数。

 

 【分析】

    观察题目条件，只不过是比例 5 多了一种L 形的骨牌，因此很自然地顺着例5  的思路走。本题中两种骨牌的最大长度和例5 一样，所以仍然用 f[i][s]表示把

前 i-1 行覆盖满、第 i 行覆盖状态为 s 的覆盖方案数，得到的递推式和例5 完全一样：

[cpp] view plaincopy

1. f[0][1…11]=1  
2. f[i][s1]=∑f[i-1][s2]  

其中(s1,s2)整体作为一个放置方案。例 5 中有两种 DFS 方案，其中第二种实现起来较第一种简单。但在本题中，新增的L形骨牌让第二种DFS难以实现，在例5中看起来有些笨拙的第一种DFS方案在本题却可以派上用场。回顾第一种DFS，我们有 5  个参数，分别为：p(当前列号)，s1、s2(当前行和对应的上一行的覆盖情况)，b1、b2(上一列的放置对当前列两行的影响，影响为 1 否则为0)。本题中，可选择的方案增多，故列表给出：

===================================================================================== 

1 其实，虽然G 为很稀疏的矩阵，但乘方后非0 元素增多，不一定还是稀疏矩阵。但对于m=n=8，计算结束后，矩阵中仍然有 80%的0，上述算法仍然高效。此处的复杂度是理想情况。

2 题目来源：SGU.131  《Hardwood Floor》

=====================================================================================

[cpp] view plaincopy

1. 覆盖情况      条件      参数s变化           参数b变化   
2.   0 0          无       s1=s1*2+b1          b1=0   
3.   0 0                   s2=s2*2+1-b2        b2=0   
4.   0 0         b1=0      s1=s1*2+1           b1=1   
5.   1 1                   s2=s2*2+1-b2        b2=0   
6.   1 0         b1=0      s1=s1*2+1           b1=0   
7.   1 0         b2=0      s2=s2*2             b2=0   
8.   1 0         b1=0      s1=s1*2+1           b1=1   
9.   1 1         b2=0      s2=s2*2             b2=0   
10.   0 1         b1=0      s1=s1*2+1           b1=1   
11.   1 1                   s2=s2*2+1-b2        b2=1   
12.   1 1         b2=0      s1=s1*2+b1          b1=1   
13.   0 1                   s2=s2*2             b2=1   
14.   1 1         b1=0      s1=s1*2+1           b1=0   
15.   1 0         b2=0      s2=s2*2             b2=1   

容易看出，在本题中此种DFS方式实现很简单。考虑其复杂度，因为L形骨牌不太规则，笔者没能找到一维的方案数的递推公式，因此无法给出复杂度的解析式。但当 m=9 时，算法共生成放置方案 79248 个，则对于 n=m=9,算法的复杂度为O(9*79248)，可以瞬间出解。和上例一样，本题也没有必要保存所有放置方案，也避免MLE。 

 

    那么，对于本题是否可以应用上题的矩阵算法呢？答案是肯定的，方法也类似，复杂度为O( 8m*logn)。然而，对于本题却不能通过上题的稀疏矩阵算法加速，原因在于刚开始时矩阵中只有1-79248/49=70%的 0，而运算结束后整个矩阵中只有2 个0，根本无法达到加速效果。

 

    由于有上题的铺垫，基本相同的本题也很快得到了解决。

 

 【例7】(1)

    给出n*m(n,m≤10)的方格棋盘，用 1*r 的长方形骨牌不重叠地覆盖这个棋盘，求覆盖满的方案数。

===================================================================================== 

1 题目来源：经典问题改编

=====================================================================================

 【分析】

    本题是例 5  的直接扩展。如果说例 5  中公式比 SCR 好，本题可以指出当公式未知时 SCR  依然是可行的算法。直接思考不容易发现方法，我们先考虑 r=3时的情况。首先，此问题有解当且仅当 m 或 n 能被3 整除。更一般的结论是：用 1*r 的骨牌覆盖满m*n 的棋盘，则问题有解当且仅当m 或n 能被r 整除。当r=2  时，则对应于例5 中m、n 至少有一个是偶数的条件。此结论的组合学证明从略。

    不同于例 5，1*3 骨牌的“攻击范围”已经达到了 3 行，可以想象例 5  中的表示方法已经无法正确表示所有状态，但其思路依然可以沿用。例 5 中用f[i][s]表示把前 i-1 行覆盖满、第 i 行覆盖状态为 s  的覆盖方案数，是因为当前行的放置方案至多能影响到上一行，状态中只要包含一行的覆盖状态即可消除后效性。本题中当前行的放置方案可以影响到上两行，故可以想到应保存两行的覆盖状态以消除后效性，即增加一维，用f[i][s1][s2]表示把前i-2 行覆盖满、第i-1 行覆盖状态为 s1、第i 行覆盖状态为 s2 的覆盖方案数。先不论上述表示方法是否可行(答案是肯定的)，r=2 时状态有2 维，r=3 时有3 维，推广后状态变量居然有r 维，这样的方法不具有推广价值，而且空间复杂度也太高。

    仔细分析上述方案，可以发现其失败之处。s1 的第p 位 s1p 为 1(覆盖)时，s2p是不可能为0 的(要求覆盖满)，则这两位(s1p , s2p )的(0,0),(0,1),(1,0),(1,1)四种组合中有一种不合法，而上述状态表示方法却冗余地保存了这个组合，造成空间复杂度过高，也进行了多余的计算1。通过上面的讨论可以知道，每一位只有 3种状态，引导我们使用三进制。我们用 f[i][s]表示把前 i-2 行覆盖满、第i-1 和第i 行覆盖状态为 s 的覆盖方案数，但这里状态s 不再是二进制，而是三进制：sp =0 表示 s1p =s2p =0；sp =1 表示 s1p =0，s2p =1；sp =2 表示 s1p =s2p =1。这样，我们就只保留了必要的状态，空间和时间上都有了改进。当r=4 时，可以类推，用四进制表示三行的状态，r=5 时用五进制……分别写出r=2,3,4,5 的程序，进行归纳，统一DFS   的形式, 可以把 DFS(p,s1,s2) 分为两部分：

[cpp] view plaincopy

1. ①for i=0 to r-1 do DFS(p+1,s1*r+i,s2*r+(i+1)mod r);  
2. ②DFS(p+r,s1*r^r +r^r-1,s2*r^r +r^r-1)   

问题解决。但DFS的这种分部方法是我们归纳猜想得到的，并没有什么道理，其正确性无法保证，我们能否通过某种途径证明它的正确性呢？仍以r=3为例。根据上面的讨论，sp  取值0 到2，表示两行第p 位的状态，但 sp并没有明确的定义。我们定义 sp  为这两行的第 p 位从上面一行开始向下连续的 1 的个数，这样的定义可以很容易地递推，递推式同上两例没有任何改变，却使得上述 DFS  方法变得很自然。

(2)

    分析算法的时间复杂度，同例 5 一样需要用到 DFS  出的方案个数h[m]，并且仿照例 5 中h[m]的递推式，我们可以得到 ： (3)

[cpp] view plaincopy

1. h[i]=r  (i=0~r-1)  
2. h[j]=r*h[j-1]+h[j-r] (j=r~m)  

理论上我们可以根据递推式得到例5中那样的精确的通项公式，但需要解高于三次的方程，且根多数为复数，无法得到例5那样简单优美的表达式，这里仅r=2..5,m=10时h[m]的值依次为：5741，77772，1077334，2609585，9784376。对于推广后的问题，例5的矩阵算法依然可行，但此时空间将是一个瓶颈。 

=====================================================================================

1 同样的，例4 中的状态表示方法也有冗余，但因为合法方案很少，又因为位运算的使用，算法仍然高效。

2 请尝试写出r=3 时用这种定义的递推式和相应的程序

3 其实从DFS 过程可以很直观地看出

=====================================================================================

【状态压缩DP部分参考代码】

以下代码纯属自己写的,如有更好的写法请赐教.谢谢.

文章已经写得很清楚了所以没有写注释.

[cpp] view plaincopy

1. // 引例  
2. #include<stdio.h>  
3. #include<string.h>  
4. const int MAXN = 1 << 20 ;  
5. __int64 f[MAXN + 23] ;  
6. int main() {  
7.     int N ;  
8.     while( scanf("%d" , &N) != EOF) {  
9.         memset( f , 0 , sizeof(f)) ;  
10.         int i , tmp = (1 << N) - 1 ;  
11.         f[0] = 1 ;  
12.         for( i = 1 ; i <= tmp ; i ++) {  
13.             int tt = i , x ;  
14.             while( tt != 0 ) {  
15.                 x = tt&(-tt) ;  
16.                 f[i] += f[i ^ x] ;  
17.                 tt ^= x ;  
18.             }  
19.         }  
20.         printf("%I64d\n" , f[tmp] ) ;  
21.     }  
22.     return 0 ;  
23. }  
24. ///////////////////////////////////////////////////////////////////////////////////////  
25. // 例 1  
26. #include<stdio.h>  
27. #include<stdlib.h>  
28. #include<string.h>  
29. #include<time.h>  
30. const int MAX_N = 1 << 8 ;  
31. __int64 f[MAX_N] ;  
32. int map[23][23] , a[23] ;  
33. int N ;  
34. void init() {  
35.     int i , j ;  
36.     srand( time ( NULL )) ;  
37.     for( i = 1 ; i <= N ; i ++) {  
38.         a[i] = 0 ;  
39.         for( j = 1 ; j <= N ; j ++) {  
40.             //scanf("%d" , &map[i][j] ) ;  
41.             map[i][j] = rand() % 2 ;  
42.             a[i] = a[i]*2 + map[i][j] ;  
43.             printf("%d " , map[i][j] ) ;  
44.         }  
45.         printf("\n") ;  
46.     }  
47.     for( i = 1 ; i <= N ; i ++)  
48.         printf("%d\n" , a[i] ) ;  
49. }  
50. void solve() {  
51.     int i , tmp = (1 << N) - 1 ;  
52.     memset( f , 0 , sizeof( f )) ;  
53.     f[0] = 1 ;  
54.     for( i = 1 ; i <= tmp ; i ++ ) {  
55.         int tt = i , j ;  
56.         int one[23] , top = 0 ;  
57.         while( tt != 0 ) {  
58.             one[++top] = tt&(-tt) ;  
59.             tt ^= one[top] ;  
60.         }  
61.         for( j = 1 ; j <= top ; j ++) {  
62.             if( ( one[j] & a[top] ) == 0 )  
63.                 f[i] += f[i^one[j]] ;  
64.         }  
65.     }  
66.     printf("%I64d\n" , f[tmp] ) ;  
67. }  
68. int main() {  
69.     while( scanf("%d" , &N) != EOF) {  
70.         init() ;  
71.         solve() ;  
72.     }  
73.     return 0 ;  
74. }  
75. //////////////////////////////////////////////////////////////////////////////////////////  
76. // 例 2  
77. #include<stdio.h>  
78. #include<string.h>  
79. int s[263] , c[263] , top  ;  
80. int f[10][36][10] ;  
81. int N , M , pn ;  
82. // 产生一行里所有的合法状态  
83. void DFS( int t , int state , int count , int* flag) {  
84.     if( t == M ) {  
85.         s[++top] = state ;  
86.         c[top] = count ;  
87.         return ;  
88.     }  
89.     flag[t + 1] = 0 ;  
90.     DFS( t + 1 , state*2 , count , flag ) ;  
91.     if( flag[t] == 0 ) {  
92.         flag[t + 1] = 1 ;  
93.         DFS( t + 1 , state*2 + 1 , count + 1 , flag) ;  
94.     }  
95. }  
96. void solve() {  
97.     memset( f , 0 , sizeof( f )) ;  
98.     f[0][1][0] = 1 ;  
99.     int i , j , k , p ;  
100.     for( i = 1 ; i <= N ; i ++) {  
101.         for( j = 1 ; j <= top ; j ++) {  
102.             for( p = 1 ; p <= top ; p ++) {  
103.                 for( k = c[j] ; k <= pn ; k ++) {  
104.                     if( (s[j]&s[p]) == 0 ) {  
105.                         f[i][j][k] += f[i - 1][p][k - c[j]] ;  
106.                     }  
107.                 }  
108.             }  
109.         }  
110.     }  
111.     int sum = 0 ;  
112.     for( i = 1 ; i <= top ; i ++) {  
113.         sum += f[N][i][pn] ;  
114.     }  
115.     printf("%d\n" , sum ) ;  
116. }  
117. int main() {  
118.     while( scanf("%d %d %d" , &N , &M , &pn ) != EOF) {  
119.         int  flag[23] = { 0 };  
120.         top = 0 ;  
121.         if( N < M ) {  
122.             N ^= M ;  
123.             M ^= N ;  
124.             N ^= M ;  
125.         }  
126.         printf("%d %d\n" , N , M ) ;  
127.         DFS( 0 , 0 , 0 , flag) ;  
128.         solve() ;  
129.     }  
130.     return 0 ;  
131. }  
132. ////////////////////////////////////////////////////////////////////////////////////////  
133. // 炮兵阵地 pku 1185  
134. #include<stdio.h>  
135. #include<string.h>  
136. int N , M ;  
137. int map[123][23] ;  
138. int s[123][100] ;  
139. int c[123][100] ;  
140. int a[123] = { 0 } ;  
141. int f[123][100][100] ;  
142. void init() {  
143.     getchar() ;  
144.     int i , j ;  
145.     char ch ;  
146.     for( i = 1 ; i <= N ; i ++) {  
147.         for( j = 1 ; j <= M ; j ++) {  
148.             scanf("%c" , &ch ) ;  
149.             map[i][j] = (ch == 'P'? 0 : 1 );  
150.         }  
151.         getchar() ;  
152.     }  
153. }  
154. // 对状态进行预处理  
155. void DFS(int id , int t , int state , int count , int* flag ) {  
156.     if( t == M ) {  
157.         s[id][++a[id]] = state ;  
158.         c[id][a[id]] = count ;  
159.         return ;  
160.     }  
161. // 判断地形 , 位置是否合法  
162.     if( map[id][t + 1] == 1 || flag[t] == 1 || ( t > 0 && flag[t - 1] == 1)) {  
163.         flag[t + 1] = 0 ;  
164.         DFS( id , t + 1 , state*2 , count , flag ) ;  
165.         return ;  
166.     }  
167.     flag[t + 1] = 0 ;  
168.     DFS( id , t + 1 , state * 2 , count , flag ) ;  
169.     flag[t + 1] = 1 ;  
170.     DFS( id , t + 1 , state * 2 + 1 , count + 1 , flag) ;  
171. }  
172. inline int Max( int x , int y ) {  
173.     return x > y ? x : y ;  
174. }  
175. int main() {  
176.     //freopen( "1185.txt" , "a+" , stdout ) ;  
177.     scanf("%d %d" , &N , &M ) ;  
178.     init() ;  
179.     int i , j , k , l , ans = 0 ;  
180.     for( i = 1 ; i <= N ; i ++) {  
181.         int flag[23] = { 0 } ;  
182.         DFS( i , 0 , 0 , 0 , flag) ;  
183.     }  
184.     a[0] = 62 ;  
185.     for( i = 1 ; i <= a[1] ; i ++ ) {  
186.         for( j = 1 ; j <= a[0] ; j ++) {  
187.             f[1][i][j] = c[1][i] ;  
188.             ans = Max( ans ,  f[1][i][j] ) ;  
189.         }  
190.     }  
191.     for( i = 2 ; i <= N ; i ++) {  
192.         for( j = 1 ; j <= a[i] ; j ++) {  
193.             for( k = 1 ; k <= a[i - 1] ; k ++) {  
194.                 int max = 0 ;  
195.                 for( l = 1 ; l <= a[i - 2] ; l ++) {  
196.                     // 判断是否合法  
197.                     if( (( s[i - 2][l]^s[i - 1][k])&s[i][j]) == 0 ) {  
198.                         max = Max( max , f[i - 1][k][l] ) ;  
199.                     }  
200.                 }  
201.                 f[i][j][k] = max + c[i][j] ;  
202.                 ans = Max( ans , f[i][j][k] ) ;  
203.             }  
204.         }  
205.     }  
206.     printf("%d\n" , ans ) ;  
207.     return 0 ;  
208. }  
209. ///////////////////////////////////////////////////  
210. // http://acm.pku.edu.cn/JudgeOnline/problem?id=2411  
211. #include<stdio.h>  
212. #include<string.h>  
213. #define MAX_N 1<<11 | 3  
214. __int64 f[2][MAX_N] ;  
215. int N , M ;  
216. // 把结果存放在数组 f[t][] 中 , s1 , s2 分别为前一行和当前行的状态 ,  
217. // m 为访问到了第几列  
218. void DFS( int t , int s1 , int s2 , int m) {  
219.     if( m == M ) {  
220.         f[t][s2] += f[1 - t][s1] ;  
221.         return ;  
222.     }  
223.     if( (m + 1)<= M ) {  
224.         // 当前行竖直放置  
225.         DFS( t , (s1 << 1) + 1 , s2 << 1 , m + 1 ) ;  
226.         // 当前行不放置  
227.         DFS( t , s1 << 1 , (s2 << 1) + 1 , m + 1 ) ;  
228.     }  
229.     if( (m + 2) <= M ) {  
230.         // 当前行水平放置  
231.         DFS( t , (s1 << 2) + 3 , (s2 << 2) + 3 , m + 2 ) ;  
232.     }  
233. }  
234. int main() {  
235.     while(scanf("%d %d" , &N , &M ) ) {  
236.         if( N == 0 && M == 0 )  
237.             break ;  
238.         //  都为奇数肯定无解  
239.         if( (( N&1 ) == 1) && (( M&1 ) == 1)) {  
240.             printf("0\n") ;  
241.             continue ;  
242.         }  
243.         if(N < M) { // 交换  
244.             int tmp = N ;  
245.             N = M , M = tmp ;  
246.         }  
247.         memset( f , 0 , sizeof( f )) ;  
248.         f[0][(1<<M) - 1] = 1 ;  
249.         int i , t = 0 ;  
250.         for( i = 1 ; i <= N ; i ++) {  
251.             DFS( 1 - t , 0 , 0 , 0 ) ;  
252.             memset( f[t] , 0 , sizeof( f[t] )) ;  
253.             t = 1 - t ;  
254.         }  
255.         printf("%I64d\n" , f[t][(1<<M) - 1]) ;  
256.     }  
257.     return 0 ;  
258. }  
259. //////////////////////////////////////////////////  
260. // 矩阵加速  
261. // http://acm.pku.edu.cn/JudgeOnline/problem?id=3420  
262. #include<stdio.h>  
263. #include<string.h>  
264. struct Matrix {  
265.     __int64 ma[16][16] ;  
266. };  
267. struct Matrix matrix ;  
268. int N ;  
269. __int64 M ;  
270. void DFS(int s1 , int s2 , int m) {  
271.     if( m == 4 ) {  
272.         matrix.ma[s1][s2] = 1 ;  
273.         return ;  
274.     }  
275.     if( (m + 1)<= 4 ) {  
276.         DFS( (s1 << 1) + 1 , s2 << 1 , m + 1 ) ;  
277.         DFS( s1 << 1 , (s2 << 1) + 1 , m + 1 ) ;  
278.     }  
279.     if( (m + 2) <= 4 ) {  
280.         DFS( (s1 << 2) + 3 , (s2 << 2) + 3 , m + 2 ) ;  
281.     }  
282. }  
283. void init() {  
284.     memset( matrix.ma , 0 , sizeof( matrix.ma ) ) ;  
285.     DFS( 0 , 0 , 0 ) ;  
286. }  
287. Matrix multiple( Matrix& a , Matrix& b ) {  
288.     int i , j , k ;  
289.     Matrix c ;  
290.     memset( c.ma , 0 , sizeof( c.ma ) ) ;  
291.     for( i = 0 ; i < 16 ; i ++) {  
292.         for( j = 0 ; j < 16 ; j ++) {  
293.             for( k = 0 ; k < 16 ; k ++) {  
294.                 c.ma[i][j] = ( c.ma[i][j] + a.ma[i][k] * b.ma[k][j] ) % M  ;  
295.             }  
296.         }  
297.     }  
298.     return c ;  
299. }  
300. Matrix solve(int n , Matrix& a) {  
301.     if( n == 1 )  
302.         return a ;  
303.     Matrix tmp = solve( n / 2 , a ) ;  
304.     tmp = multiple( tmp , tmp ) ;  
305.     if( (n % 2) == 1 )  
306.         tmp = multiple( tmp , a ) ;  
307.     return tmp ;  
308. }  
309. int main() {  
310.     init() ;  
311.     while( scanf("%I64d %I64d" , &N , &M ) ) {  
312.         if( N == 0 && M == 0 )  
313.             break ;  
314.         if( M == 1 ) {  
315.             printf("0\n") ;  
316.             continue ;  
317.         }  
318.         Matrix ans = solve(N , matrix ) ;  
319.         printf("%I64d\n" , ans.ma[15][15] ) ;  
320.     }  
321.     return 0 ;  
322. }