# 木塊砌墻 作者：July、caopengcs、紅色標記。致謝：fuwutu、demo。 時間：二零一三年八月十二日 __題目__：用 1×1×1, 1×2×1以及2×1×1的三種木塊（橫綠豎藍，且綠藍長度均為2），  搭建高長寬分別為K × 2^N × 1的墻，不能翻轉、旋轉（其中，0<=N<=1024，1<=K<=4）  有多少種方案，輸出結果 對1000000007取模。 舉個例子如給定高度和長度：N=1 K=2，則答案是7，即有7種搭法，如下圖所示：  __詳解__：此題很有意思，涉及的知識點也比較多，包括動態規劃，快速矩陣冪，狀態壓縮，排列組合等等都一一考察了個遍。而且跟一個比較經典的矩陣乘法問題類似：即用1 x 2的多米諾骨牌填滿M x N的矩形有多少種方案，M<=5，N<2^31，輸出答案mod p的結果  OK，回到正題。下文使用的圖示說明(所有看到的都是橫切面)：  首先說明“？方塊”的作用  “？方塊”，表示這個位置是空位置，可以任意擺放。 上圖的意思就是，當右上角被綠色木塊占用，此位置固定不變，其他位置任意擺放，在這種情況下的堆放方案數。 ###解法一、窮舉遍歷 初看此題，你可能最先想到的思路便是窮舉：用二維數組模擬墻，從左下角開始擺放，從左往右，從下往上，最后一個格子是右上角那個位置；每個格子把每種可以擺放木塊都擺放一次，每堆滿一次算一種用擺放方法。為了便于描述，為木塊的每個格子進行編號：  下面演示當n=1,k=2的算法過程（7種情況）：  窮舉遍歷在數據規模比較小的情況下還撐得住，但在0<=N<=1024這樣的數據規模下，此方法則立刻變得有心無力，因此我們得尋找更優化的解法。 ###解法二、遞歸分解 遞歸求解就是把一個大問題，分解成小問題，逐個求解，然后再解決大問題。 ####2.1、算法演示 假如有墻規模為(n,k)，如果從中間切開，被分為規模問(n-1,k)的兩堵墻，那么被分開的墻和原墻有什么關系呢？我們首先來看一下幾組演示。 #####2.1.1、n=1，k=2的情況 首先演示，__n=1,k=2__時的情況，如下圖2-1：  圖 2-1 上圖2-1中：  表示，左邊墻的所有堆放方案數 * 右邊墻所有堆放方案數 = 2 * 2 = 4  表示，當切開處有一個橫條的時候，空位置存在的堆放方案數。左邊*右邊 = 1*1 = 2；剩余兩組以此類推。 這個是排列組合的知識。 #####2.1.2、n=2，k=3的情況 其次，我們再來演示下面更具一般性的計算分解，即當__n=2,k=3__的情況，如下圖2-2：  圖 2-2 再從分解的結果中，挑選一組進行分解演示：  圖 2-3 通過圖2-2和圖2-3的分解演示，可以說明，最終都是分解成一列求解。在逐級向上匯總。 #####2.1.3、n=4，k=3的情況 我們再假設一堵墻n=4，k=3，也就是說，寬度是16，高度是3時，會有以下分解：  圖2-4 根據上面的分解的一個中間結果，再進行分解，如下：  圖2-5 通過上面圖2-1~圖2-5的演示可以明確如下幾點： 1.假設f(n)用于計算問題，那么f(n)依賴于f(n-1)的多種情況。 2.切開處有什么特殊的地方呢？通過上面的演示，我們得知被切開的兩堵墻從沒有互相嵌入的木塊（綠色木塊）到全是互相連接的木塊，相當于切口綠色木塊的全排列（即有綠色或者沒有綠色的所有排列），即有2^k種狀態（比如k=2，且有綠色用1表示，沒有綠色用0表示，那么就有00、01、10、11這4種狀態）。根據排列組合的性質，把每一種狀態下左右木墻堆放方案數相乘，再把所有乘積求和，就得到木墻的堆放結果數。以此類推，將問題逐步往下分解即可。 3.此外，從圖2-5中可以看出，除了需要考慮切口綠色木塊的狀態，還需要考慮最左邊一列和最右邊一列的綠色木塊狀態。我們把這兩種邊界狀態稱為左邊界狀態和右邊界狀態，分別用leftState和rightState表示。 且在觀察圖2-5被切分后，所有左邊的墻，他們的左邊界ls狀態始終保持不變，右邊界rs狀態從0~maxState, maxState = 2^k-1（有綠色方塊表示1，沒有表示0；ls表示左邊界狀態，rs表示右邊界狀態）：  圖2-6 同樣可以看出右邊的墻的右邊界狀態保持不變，而左邊界狀態從0~maxState。要堆砌的木墻可以看做是左邊界狀態=0，和右邊界狀態=0的一堵墻。 有一點可能要特別說明下，即上文中說，有綠色方塊的狀態表示標為1，無綠色方塊的狀態表示標為0，特意又拿上圖2-6標記了一些數字，以讓絕大部分讀者能看得一目了然，如下所示：  圖2-7 這下，你應該很清楚的看到，在上圖中，左邊木塊的狀態表示一律為010，右邊木塊的狀態表示則是000~111（即從下至上開始計數，右邊木塊rs的狀態用二進制表示為：000 001 010 011 100 101 110 111，它們各自分別對應整數則是：0 1 2 3 4 5 6 7）。 ####2.2、計算公式 通過圖2-4、圖2-5、圖2-6的分解過程，我們可以總結出下面公式（leftState=最左邊邊界狀態，rightState=最右邊邊界狀態）：  即：  接下來，分3點解釋下上述公式： __1__、上述函數返回結果是當左邊狀態為=leftState，右邊狀態=rightState時木墻的堆砌方案數，相當于直接分解的左右狀態都為0的情況，即直接分解f(n,k,0,0)即可。看到這，讀者可能便有疑問了，既然直接分解f(n,k,0,0)即可，為何還要加leftstate和leftstate兩個變量呢？回顧下2.1.3節中n=4，k=3的演示例子，即當n=4，k=3時，其分解過程即如下圖（上文2.1.3節中的圖2-4）  也就是說，剛開始直接分解f(4,3,0,0)，即n=4，k=3，leftstate=0，rightstate=0，但分解過程中leftstate和rightstate皆從0變化到了maxstate，故才讓函數的第3和第4個參數采用leftstate和rightstate這兩個變量的形式，公式也就理所當然的寫成了f(n,k,leftstate,rightstate)。 __2__、然后我們再看下當n=4，k=3分解的一個中間結果，即給定如上圖最下面部分中紅色框框所框住的木塊時：  它用方程表示即為 f(2,3,2,5)，怎么得來的呢？其實還是又回到了上文2.1.3節中，當n=2，k=3 時（下圖即為上文2.1.3節中的圖2-5和圖2-6）   左邊界ls狀態始終保持不變時，右邊界rs狀態從0~maxState；右邊界狀態保持不變時，而左邊界狀態從0~maxState。 故上述分解過程用方程式可表示為： __f(2,3,2,5) = f(1,3,2,0) * f(1,3,0,5)__ __+ f(1,3,2,1) * f(1,3,1,5)__ __+ f(1,3,2,2) * f(1,3,2,5)__ __+ f(1,3,2,3) * f(1,3,3,5)__ __+ f(1,3,2,4) * f(1,3,4,5)__ __+ f(1,3,2,5) * f(1,3,5,5)__ __+ f(1,3,2,6) * f(1,3,6,5)__ __+ f(1,3,2,7) * f(1,3,7,5)__ 說白了，我們曾在2.1節中從圖2-2到圖2-6正推推導出了公式，然上述過程中，則又再倒推推了一遍公式進行了說明。 __3__、最后，作者是怎么想到引入 leftstate 和rightstate 這兩個變量的呢？如紅色標記所說："因為切開后，發現綠色條，在分開處不斷的變化，當時也進入了死胡同，我就在想，藍色的怎么辦。后來才想明白，與藍色無關。每一種變化就是一種狀態，所以就想到了引入leftstate 和rightstate這兩個變量。" ####2.3、參考代碼 下面代碼就是根據上面函數原理編寫的。最終執行效率，n=1024,k=4 時，用時0.2800160秒（之前代碼用的是字典作為緩存，用時在1.3秒左右，后來改為數組結果，性能大增）。"" ```cs //copyright@紅色標記 12/8/2013 //updated@July 13/8/2013 using System; using System.Collections.Generic; using System.Text; using System.Collections; namespace HeapBlock { public class WoolWall { private int n; private int height; private int maxState; private int[, ,] resultCache; //結果緩存數組 public WoolWall(int n, int height) { this.n = n; this.height = height; maxState = (1 << height) - 1; resultCache = new int[n + 1, maxState + 1, maxState + 1]; //構建緩存數組，每個值默認為0； } /// <summary> /// 靜態入口。計算堆放方案數。 /// </summary> /// <param name="n"></param> /// <param name="k"></param> /// <returns></returns> public static int Heap(int n, int k) { return new WoolWall(n, k).Heap(); } /// <summary> /// 計算堆放方案數。 /// </summary> /// <returns></returns> public int Heap() { return (int)Heap(n, 0, 0); } private long Heap(int n, int lState, int rState) { //如果緩存數組中的值不為0，則表示該結果已經存在緩存中。 //直接返回緩存結果。 if (resultCache[n, lState, rState] != 0) { return resultCache[n, lState, rState]; } //在只有一列的情況，無法再進行切分 //根據列狀態計算一列的堆放方案 if (n == 0) { return CalcOneColumnHeapCount(lState); } long result = 0; for (int state = 0; state <= maxState; state++) { if (n == 1) { //在只有兩列的情況，判斷當前狀態在切分之后是否有效 if (!StateIsAvailable(n, lState, rState, state)) { continue; } result += Heap(n - 1, state | lState, state | lState) //合并狀態。因為只有一列，所以lState和rState相同。 * Heap(n - 1, state | rState, state | rState); } else { result += Heap(n - 1, lState, state) * Heap(n - 1, state, rState); } result %= 1000000007;//為了防止結果溢出，根據題目要求求模。 } resultCache[n, lState, rState] = (int)result; //將結果寫入緩存數組中 resultCache[n, rState, lState] = (int)result; //對稱的墻結果相同，所以直接寫入緩存。 return result; } /// <summary> /// 根據一列的狀態，計算列的堆放方案數。 /// </summary> /// <param name="state">狀態</param> /// <returns></returns> private int CalcOneColumnHeapCount(int state) { int sn = 0; //連續計數 int result = 1; for (int i = 0; i < height; i++) { if ((state & 1) == 0) { sn++; } else { if (sn > 0) { result *= CalcAllState(sn); } sn = 0; } state >>= 1; } if (sn > 0) { result *= CalcAllState(sn); } return result; } /// <summary> /// 類似于斐波那契序列。 /// f(1)=1 /// f(2)=2 /// f(n) = f(n-1)*f(n-2); /// 只是初始值不同。 /// </summary> /// <param name="k"></param> /// <returns></returns> private static int CalcAllState(int k) { return k <= 2 ? k : CalcAllState(k - 1) + CalcAllState(k - 2); } /// <summary> /// 判斷狀態是否可用。 /// 當n=1時，分割之后，左墻和右邊墻只有一列。 /// 所以state的狀態碼可能會覆蓋原來的邊緣狀態。 /// 如果有覆蓋，則該狀態不可用；沒有覆蓋則可用。 /// 當n>1時，不存在這種情況，都返回狀態可用。 /// </summary> /// <param name="n"></param> /// <param name="lState">左邊界狀態</param> /// <param name="rState">右邊界狀態</param> /// <param name="state">切開位置的當前狀態</param> /// <returns>狀態有效返回 true,狀態不可用返回 false</returns> private bool StateIsAvailable(int n, int lState, int rState, int state) { return (n > 1) || ((lState | state) == lState + state && (rState | state) == rState + state); } } } ``` 上述程序中， * WoolWall.Heap(1024,4); //直接通過靜態方法獲得結果 * new WoolWall(n, k).Heap();//通過構造對象獲得結果 #####2.3.1、核心算法講解 因為它最終都是分解成一列的情況進行處理，這就會導致很慢。為了提高速度，本文使用了緩存機制來提高性能。緩存原理就是，n,k,leftState,rightState相同的墻，返回的結果肯定相同。利用這個特性，每計算一種結果就放入到緩存中，如果下次計算直接從緩存取出。剛開始緩存用字典類實現，有網友給出了更好的緩存方法——數組。這樣性能好了很多，也更加簡單。程序結構如下圖所示：  上圖反應了Heap調用的主要方法調用，在循環中，result 累加 lResult 和 rResult。 ①在實際代碼中，首先是從緩存中讀取結果，如果沒有從緩存中讀取結果再進行計算。 分解到一列時，不再分解，直接計算結果 ```cs if (n == 0) { return CalcOneColumnHeap(lState); } ``` ②下面是整個程序的核心代碼，通過for循環，求和state=0到state=2^k-1的兩邊木墻乘積： ```cs for (int state = 0; state <= maxState; state++) { if (n == 1) { if (!StateIsAvailable(n, lState, rState, state)) { continue; } result += Heap(n - 1, state | lState, state | lState) * Heap(n - 1, state | rState, state | rState); } else { result += Heap(n - 1, lState, state) * Heap(n - 1, state, rState); } result %= 1000000007; } ``` 當n=1切分時，需要特殊考慮。如下圖：  圖2-8 看上圖中，因為左邊墻中間被綠色方塊占用，所以在（1,0）-（1,1）這個位置（位置的標記方法同解法一）不能再放綠色方塊。所以一些狀態需要排除，如state=2需要排除。同時在還需要合并狀態，如state=1時，左邊墻的狀態=3。 特別說明下：依據我們上文2.2節中的公式，如果第i行有這種木塊，state對應2^(i-1)，加上所有行的貢獻就得到state（0就是沒有這種橫跨木塊，2^k-1就是所有行都是橫跨木塊），然后遍歷state，還記得上文中的圖2-7么？  當第i行被這樣的木塊或這樣的木塊占據時，其各自對應的state值分別為： 1.當第1行被占據，state=1； 2.當第2行被占據，state=2； 3.當第1和第2行都被占據，state=3； 4.當第3行被占據，state=4； 5.當第1和第3行被占據，state=5； 6.當第2和第3行被占據，state=6； 7.當第1、2、3行全部都被占據，state=7。 至于原因，即如2.1.3節節末所說：二進制表示為：000 001 010 011 100 101 110 111，它們各自分別對應整數則是：0 1 2 3 4 5 6 7。 具體來說，下面圖中所有框出來的位置，不能有綠色的：  ③CalcOneColumnHeap(int state)函數用于計算一列時擺放方案數。 計算方法是， 求和被綠色木塊分割開的每一段連續方格的擺放方案數。每一段連續的方格的擺放方案通過CalcAllState方法求得。經過分析，可以得知CalcAllState是類似斐波那契序列的函數。 舉個例子如下（分步驟講述）： 1.令state = 4546（state=2^k-1，k最大為4，故本題中state最大在15，而這里取state=4546只是為了演示如何計算），二進制是：1000111000010。位置上為1，表示被綠色木塊占用，0表示空著，可以自由擺放。 2.1000111000010 被分割后 1 000 111 0000 1 0, 那么就有 000=3個連續位置， 0000=4個連續位置 ， 0=1個連續位置。 3.堆放結果=CalcAllState(3) + CalcAllState(4) + CalcAllState(1) = 3 + 5 + 1 = 9。 ####2.4、再次優化 上面程序因為調用性能的樹形結構，形成了大量的函數調用和緩存查找，所以其性能不是很高。 為了得到更高的性能，可以讓所有的運算直接依賴于上一次運算的結果，以防止更多的調用。即如果每次運算都算出所有邊界狀態的結果，那么就能為下一次運算提供足夠的信息。后續優化請[查閱此文第3節](http://blog.csdn.net/dw14132124/article/details/9038417#t2)。 ###解法三、動態規劃 相信讀到上文，不少讀者都已經意識到這個問題其實就是一個動態規劃問題，接下來咱們換一個角度來分析此問題。 ####3.1、暴力搜索不可行 首先，因為木塊的寬度都是1，我們可以想成2維的問題。也就是說三種木板的規格分別為1* 1, 1 * 2, 2 * 1。 通過上文的解法一，我們已經知道這個問題最直接的想法就是暴力搜索，即對每個空格嘗試放置哪種木板。但是看看數據規模就知道，這種思路是不可行的。因為有一條邊范圍長度高達2^1024，普通的電腦，2^30左右就到極限了。于是我們得想想別的方法。 ####3.2、另辟蹊徑 為了方便，我們把墻看做有2^n行，k列的矩形。這是因為雖然矩形木塊不能翻轉，但是我們同時擁有1*2和2*1的兩種木塊。 假設我們從上到下，從左到右考慮每個1*1的格子是如何被覆蓋的。顯然，我們每個格子都要被覆蓋住。木塊的特點決定了我們覆蓋一個格子最多只會影響到下一行的格子。這就可以讓我們暫時只考慮兩行。 假設現我們已經完全覆蓋了前(i–1)行。那么由于覆蓋前(i-1)行導致第i行也不“完整”了。如下圖： xxxxxxxxx ooxooxoxo 我們用x表示已經覆蓋的格子，o表示沒覆蓋的格子。為了方便，我們使用9列。 我們考慮第i行的狀態，上圖中，第1列我們可以用1*1的覆蓋掉，也可以用1*2的覆蓋前兩列。第4、5列的覆蓋方式和第1、2列是同樣的情況。第7列需要覆蓋也有兩種方式，即用1*1的覆蓋或者用2*1的覆蓋，但是這樣會導致第（i+1)行第7列也被覆蓋。第9列和第7列的情況是一樣的。這樣把第i行覆蓋滿了之后，我們再根據第(i+1)行被影響的狀態對下一行進行覆蓋。 那么每行有多少種狀態呢？顯然有2^k，由于k很小，我們只有大約16種狀態。如果我們對于這些狀態之間的轉換制作一個矩陣，矩陣的第i行第j列的數表示的是我們第m行是狀態i，我們把它完整覆蓋掉，并且使得第(m + 1)行變成狀態j的可能的方法數，這個矩陣我們可以暴力搜索出來，搜索的方式就是枚舉第m行的狀態，然后嘗試放木板，用所有的方法把第m行覆蓋掉之后，下一行的狀態。當然，我們也可以認為只有兩行，并且第一行是2k種狀態的一種，第二行起初是空白的，求使得第一行完全覆蓋掉，第二行的狀態有多少種類型以及每種出現多少次。 ####3.3、動態規劃 這個矩陣作用很大，其實我們覆蓋的過程可以認為是這樣：第一行是空的，我們看看把它覆蓋了，第2行是什么樣子的。根據第二行的狀態，我們把它覆蓋掉，看看第3行是什么樣子的。 如果我們知道第i行的狀態為s,怎么考慮第i行完全覆蓋后，第(i+1)行的狀態？那只要看那個矩陣的狀態s對應的行就可以了。我們可以考慮一下，把兩個這樣的方陣相乘得到得結果是什么。這個方陣的第i行第j個元素是這樣得到的，是第i行第k個元素與第k行第j個元素的對k的疊加。它的意義是上一行是第m行是狀態i，把第m行和第(m+1)行同時覆蓋住，第(m+2)行的狀態是j的方法數。這是因為中間第(m+1)行的所有狀態k，我們已經完全遍歷了。 于是我們發現，每做一次方陣的乘法，我們相當于把狀態推動了一行。那么我們要坐多少次方陣乘法呢？就是題目中墻的長度2ⁿ,這個數太大了。但是事實上，我們可以不斷地平方n次。也就是說我們可以算出A²,A⁴, A⁸, A¹⁶……方法就是不斷用結果和自己相乘，這樣乘n次就可以了。 因此，我們最關鍵的問題就是建立矩陣A。我們可以這樣表示一行的狀態，從左到右分別叫做第0列，第1列……覆蓋了我們認為是1，沒覆蓋我們認為是0，這樣一行的狀態可以表示為一個整數。某一列的狀態我們可以用位運算來表示。例如，狀態x第i列是否被覆蓋，我們只需要判斷x & (1 << i) 是否非0即可，或者判斷(x >> i) & 1， 用右移位的目的是防止溢出，但是本題不需要考慮溢出，因為k很小。 接下來的任務就是遞歸嘗試放置方案了 ####3.4、參考代碼 最終結果，我們最初的行是空得，要求最后一行之后也不能被覆蓋，所以最終結果是矩陣的第[0][0]位置的元素。另外，本題在乘法過程中會超出32位int的表示范圍，需要臨時用C/C++的long long，或者java的long。 參考代碼如下： ```c //copyright@caopengcs 12/08/2013 #ifdef WIN32 #define ll __int64 #else #define ll long long #endif // 1 covered 0 uncovered void cal(int a[6][32][32],int n,int col,int laststate,int nowstate) { if (col >= n) { ++a[n][laststate][nowstate]; return; } //不填 或者用1*1的填 cal(a,n, col + 1, laststate, nowstate); if (((laststate >> col) & 1) == 0) { cal(a,n, col + 1, laststate, nowstate | (1 << col)); if ((col + 1 < n) && (((laststate >> (col + 1)) & 1) == 0)) { cal(a,n, col + 2, laststate, nowstate); } } } inline int mul(ll x, ll y) { return x * y % 1000000007; } void multiply(int n,int a[][32],int b[][32]) { // b = a * a int i,j, k; for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) { for (k = b[i][j] = 0; k < n; ++k) { if ((b[i][j] += mul(a[i][k],a[k][j])) >= 1000000007) { b[i][j] -= 1000000007; } } } } } int calculate(int n,int k) { int i, j; int a[6][32][32],mat[2][32][32]; memset(a,0,sizeof(a)); for (i = 1; i <= 5; ++i) { for (j = (1 << i) - 1; j >= 0; --j) { cal(a,i, 0, j, 0); } } memcpy(mat[0], a[k],sizeof(mat[0])); k = (1 << k); for (i = 0; n; --n) { multiply(k, mat[i], mat[i ^ 1]); i ^= 1; } return mat[i][0][0]; } ``` ##參考鏈接及推薦閱讀 1. caopengcs，[木塊砌墻](http://blog.csdn.net/caopengcs/article/details/9928061) 2. 紅色標記，[木塊砌墻](http://blog.csdn.net/dw14132124/article/details/9038417) 3. LoveHarvy，[木塊砌墻](http://blog.csdn.net/wangyan_boy/article/details/9131501) 4. [在線編譯測試木塊砌墻問題](http://hero.pongo.cn/Question/Details?ID=36&ExamID=36) 5. hero上[木塊砌墻一題](http://hero.pongo.cn/Question/Details?ExamID=36&ID=36&bsh_bid=273040296)