CSP 201409-5 拼图问题（给出一个n×m的方格图，现在要用如下L型的积木拼到这个图中......）

CSP 201409-5 拼图问题

一、题目信息

第一次写博客，有什么疏漏之处，欢迎各位大佬指出<(*￣▽￣*)/

题目要求

试题编号	201409-5
试题名称	拼图
时间限制	3.0s
试题名称	256.0MB

题目描述

给出一个n×m的方格图，现在要用如下L型的积木拼到这个图中，使得方格图正好被拼满，请问总共有多少种拼法。其中，方格图的每一个方格正好能放积木中的一块。积木可以任意旋转。给出一个n×m的方格图，现在要用如下L型的积木拼到这个图中，使得方格图正好被拼满，请问总共有多少种拼法。其中，方格图的每一个方格正好能放积木中的一块。积木可以任意旋转。

输入格式
　　输入的第一行包含两个整数n, m，表示方格图的大小。
输出格式
　　输出一行，表示可以放的方案数，由于方案数可能很多，所以请输出方案数除以1,000,000,007的余数。
样例输入
6 2
样例输出
4
样例说明
　　四种拼法如下图所示：

评测用例规模与约定
　　在评测时将使用10个评测用例对你的程序进行评测。
　　评测用例1和2满足：1<=n<=30，m=2。
　　评测用例3和4满足：1<=n, m<=6。
　　评测用例5满足：1<=n<=100，1<=m<=6。
　　评测用例6和7满足：1<=n<=1000，1<=m<=6。
　　评测用例8、9和10满足：1<=n<=10^15，1<=m<=7

二、题目分析

一般读完题目，可能第一反应试试暴力搜索，但从题目的位置（这是最后一道题）和给定的数据量级来看，是不可能用暴力搜索搜索解决该问题的，因为该题目中的方格图行数 row<=10¹⁵，每行共有的可能情况为2^m次方，所以方格图的总状态数目为 (2^m)ⁿ(n为行数，m为列数），这是一个非常恐怖的数字，当n比较大时，即使交给“神威太湖之光”来计算，也不可能是计算的出来。因此，还是老老实实解题吧(*￣︶￣)~

但在在做之前要考虑以下极为重要的两点问题：

用什么数据结构表示方格图（或者说每个格子）的状态？
如此量级的数值取值范围，如何提高效率？

做过此类题目的人应该比较敏感，其实这就是动态规划中经典的状态压缩模型。

状态表示
在这里用一个二进制数表示每一行格子的状态，如(1010_0000)₂，表示一行格子中第1，3个格子被填充，其余空白，如图：在这里用一个二进制数表示每一行格子的状态，如(1010_0000)₂，表示一行格子中第1，3个格子被填充，其余空白，如图：

所以，如果这是第i行，这里这一行的状态可以表示为d[i]=160,但是这里有n行，难道真的要用长度为n的数组保存状态，当然不是，既然是动态规划，必然是分析出该题目中对应的状态转移方程，用来表示放置积木的导致的状态转移过程。
状态转移矩阵
    其实很容易发现，为什么题目中只给出了一个"L"型的积木（可以加上高度为2的“I”型积木），而没有给高度为4的“I”型和高度为3的“L”型的积木(PS:我说的这些积木是什么样子，就不给出图片了，相信玩过俄罗斯方块的都知道)，是因为高度 <=2 ，上一行的放置操作只会影响到下一行的格子状态，而不会对下下行产生任何直接影响，因此很方便我们使用矩阵表示状态转移,因而，此处我们用一个二维数组status[][]表示状态转移矩阵，（划重点）status[i][j]表示的是在当前行初始状态为i（转换成二进制）的情况下，将当前行放满并且下一行最终状态为j（转换为二进制）的方案数（划重点）。例如status[0][105]（105二进制表示:01101001）表示下边的状态

    以为在i行初始状态状态为0(为空)的状态下，填满第i行并且让下一行状态为105的放置方案数量，再如status[9][41]（41的二进制:00101001）

这表示在第i行初始状态为9的情况下填满第i行并且让下一行最终状态为41的方案数量。
    分析了这么久，发现这道题渐渐有了眉目，如果我能计算出status[i][j] (0<=i<m,0<=j<m)下的所有可能组合方案数量，然后通过该该状态转移矩阵，由第1行初始状态为空，将其放满的所有可能情况作为第2行的初始状态值，计算第3行的可能最终状态值，由此递推，每完全放满一行，便进行一次状态转移，如果保证最后一行放满，岂不是能得到最终的题目所求的方案情况？下边我们进行求解状态转移矩阵

状态转移矩阵计算
状态转移矩阵的计算还是比较单纯的，直接利用DFS搜索即可，但是要注意一些问题，因为我们是用0~m(m<128=2⁷)表示的每一行格子的状态，那么改变格子的状态、状态检测等必然要涉及到位运算，如果不是很熟悉，要先去简单学习一下，此处不做讲解。
我们的计算思路是这样的，当前行格子的初始状态共有0~2^m-1的可能，也就是共有2^m种状态数目，所以每一种初始状态下将改行放满的方案数量直接循环2^m次即可。代码如下：

    #define FOR(i,start,end,step) for(int i=(start);i<(end);i+=(step))FOR(S,0,1<<m,1) search(S,S,0,0);void search(int S,int S1,int S2,int col){if(S1==(1<<m)-1){status[S][S2]++;return; }if(!check(S1,col)){// 注释1：如果col列不为空 search(S,S1,S2,col+1); }else{//注释:2if(check(S1,col+1)&&check(S2,col)) //"枪"search(S,put(S1,col,col+1),put(S2,col,-1),col+2);if(check(S1,col+1)&&check(S2,col+1)) //"7"search(S,put(S1,col,col+1),put(S2,col+1,-1),col+2);if(check(S2,col-1)&&check(S2,col)) //"J"search(S,put(S1,col,-1),put(S2,col-1,col),col+1);if(check(S2,col)&&check(S2,col+1)) //"L"search(S,put(S1,col,-1),put(S2,col,col+1),col+1);}}

简单解释一下代码：
（1）参数解释：S代表当前行的初始状态，S1代表当前行的此时状态，S2代表下一行的此时状态，col代表当前处在第几列（列数从1开始至m）。
（2）递归搜索的出口:当前行被放满，并做方案数加一或者单纯的列数遍历完回退。
（3）注释1：检查当前col列是否为空，如果不为空，直接本层递归，否则再判断是否可以放置某个积木。
（4）注释2：因为"L型"积木共有4中方向，因此需要尝试4中可能。
假设n=2,m=3,计算结果如下：

从结果可以看出（以status[1][2]==1为例,注意下标从0开始），在当前行行初始状态为1（001），保证下一行最终状态为2（010），将当前行放满的方案数量为1

4.方案数计算

说到此处了，如何进行放置n行后的方案数量的计算呢？实际上很简单，直接计算矩阵乘法即可，还是以上边的n=2,m=3的情况为例说明，如果求得result=status*status，则result[0][0]代表**进行1次状态转移**，保证第1行初始状态为0，第3行（m=2，3=m+1）刚好全空的（为0）的方案总数量（答案为result[0][0] =2），或者可以理解为**进行0次状态转移**，status[0][m]=2页即为所求，代表着第一行初始状态为0，第2行（m=2）**最终状态为m(全为1，最后一行全满)**。假设**n为行数，m为列数**，所以这里有两种理解方式：
（1）进行 **m-2 次状态转移**，也就是计算result=status^m-1^ ，最终让**第m行放满**（最后一行）的方案数量为result[0][m];
（2）进行**m-1状态转移**，也就是计算result=status^m^，最终让**第m+1行全空**（想象它存在，毕竟让它全空，不影响什么）的方案数量为result[0][0]; （PS：虽然在此题目中这两种都可以，但在下边会证明其实第2种才是最好的理解方式）

5. 快速幂与矩阵快速幂

走到这一步，基本上思路就十分清晰了，但还有一个问题，题目给出的行数上限为**10的15次方**，这是一个非常恐怖的数字，如果要计算这么多次128*128（根据你开辟的空间而定）的大矩阵相乘，肯定会超时，所以，这里要用上矩阵快速幂的知识，鉴于时间和篇幅的原因（(*^▽^*)其实是本人嫌麻烦<(*￣▽￣*)/），这里就不介绍快速幂和矩阵快速幂了，大家可以去(https://www.cnblogs.com/cmmdc/p/6936196.html)这里学习，我这里直接贴出矩阵快速幂的代码：

 long long t=n<<1;while(t>>=1>0){if((t&1)!=0) mat_mul(res,status);mat_mul(status,status);  }//矩阵相乘 void mat_mul(int A[1<<MAX_COL][1<<MAX_COL],int B[1<<MAX_COL][1<<MAX_COL]){long long matrix[1<<MAX_COL][1<<MAX_COL]={0};//临时矩阵FOR(i,0,1<<MAX_COL,1)FOR(j,0,1<<MAX_COL,1)FOR(k,0,1<<MAX_COL,1)matrix[i][j]=(matrix[i][j]+(long long)A[i][k]*(long long)B[k][j])%MOD;FOR(i,0,1<<MAX_COL,1) FOR(j,0,1<<MAX_COL,1) A[i][j]=(int)matrix[i][j];}

总结
到此为止，解决该题目的大问题基本上都被破解了，剩下就是写代码了♪(^∇*)，这里我给处C语言和java版的代码。补充一下，经过实测，当n取最大值时，运行时间约为0.5~0.6秒，完全在3秒的可控范围之内^_^：
C语言版：

#include <stdio.h>#define MOD 1000000007
#define MAX_COL 7#define FOR(i,start,end,step) for(int i=(start);i<(end);i+=(step))/*@author 阡陌红尘
*/
long long n;//方格图的行数
int m;//方格图的行列数
int status[1<<MAX_COL][1<<MAX_COL]; //状态转移矩阵
int res[1<<MAX_COL][1<<MAX_COL];//保存结果状态的矩阵
int put(int S,int col1,int col2){S|=1<<col1;//改变当前行/下一行，第col1列格子的状态if(col2!=-1) S|=1<<col2;//如果传值不为-1，则也改变col2列格子的状态return S;
}
int check(int S,int col){//判断当前行/下一行col列位置是否为空if(col>=0&&col<m&&(S&(1<<col))==0) return 1;return 0;
}
void search(int S,int S1,int S2,int col){if(S1==(1<<m)-1){//如果当前行被放满，进行状态转移status[S][S2]++;return; }if(!check(S1,col)){//如果col列不为空 search(S,S1,S2,col+1); }else{if(check(S1,col+1)&&check(S2,col)) //"枪"search(S,put(S1,col,col+1),put(S2,col,-1),col+2);if(check(S1,col+1)&&check(S2,col+1)) //"7"search(S,put(S1,col,col+1),put(S2,col+1,-1),col+2);if(check(S2,col-1)&&check(S2,col)) //"J"search(S,put(S1,col,-1),put(S2,col-1,col),col+1);if(check(S2,col)&&check(S2,col+1)) //"L"search(S,put(S1,col,-1),put(S2,col,col+1),col+1);}
}
//矩阵相乘
void mat_mul(int A[1<<MAX_COL][1<<MAX_COL],int B[1<<MAX_COL][1<<MAX_COL]){long long matrix[1<<MAX_COL][1<<MAX_COL]={0};//临时矩阵FOR(i,0,1<<MAX_COL,1)FOR(j,0,1<<MAX_COL,1)FOR(k,0,1<<MAX_COL,1)matrix[i][j]=(matrix[i][j]+(long long)A[i][k]*(long long)B[k][j])%MOD;FOR(i,0,1<<MAX_COL,1) FOR(j,0,1<<MAX_COL,1) A[i][j]=(int)matrix[i][j];
}
int main(void){scanf("%lld %d",&n,&m);FOR(S,0,1<<m,1) search(S,S,0,0);//构建单位矩阵FOR(i,0,1<<MAX_COL,1) FOR(j,0,1<<MAX_COL,1) res[i][j]=(i==j)?1:0;   //矩阵快速幂long long t=n<<1;while(t>>=1>0){//求status矩阵的n次方，n-1相乘，n-1次状态转移if((t&1)!=0) mat_mul(res,status);mat_mul(status,status);  }printf("%d",res[0][0]);
}

java版:

import java.util.Scanner;/*** @author 阡陌红尘**/
public class StackBlock {static long n;//行数 n <=Math.pow(10, 15);static int m;//列数 m <=7final int MOD=1000000007;//<=2147483647 = 2^31 - 1static final int MAX_COL=7;//status[i][j]指的是在当前行初始状态为i(转换为2进制)的情况下，把当前行填满并且下一行出现j状态(2进制)的方案数量static int[][] status=new int[1<<MAX_COL][1<<MAX_COL];//状态转移矩阵static int[][] result=new int[1<<MAX_COL][1<<MAX_COL];//存贮最终方案数量的矩阵public static void main(String[] args){StackBlock stackBlock=new StackBlock();Scanner scanner=new Scanner(System.in);n=scanner.nextLong();//1<=n<=10^15m=scanner.nextInt();//1<=m<=MAX_COLlong startTime=System.currentTimeMillis();//在m列中，共有2^m种状态,遍历搜索序号（状态）为0~(2^m)-1的可能for(int S=0;S<(1<<m);S++) stackBlock.search(S,S,0,0);for(int i=0;i<MAX_COL;i++)for(int j=0;j<MAX_COL;j++)result[i][j]=(i==j)?1:0;//矩阵快速幂 long b=n;//求status矩阵的n次方，n-1相乘，n-1次状态转移while(b>0){if((b&1)!=0) stackBlock.matMul(result,status);stackBlock.matMul(status,status);b>>=1;}//输出第1行原状态为全空0000_0000到放置后第n+1行为全空0000_0000并且保证1~n行全满的的方案数量System.out.println(result[0][0]);long endTime=System.currentTimeMillis();System.out.println("运行时间:"+(endTime-startTime)/1000.0+"秒");}//S:状态源 ，S1:当前行状态，S2:下一行状态public void search(int S,int S1,int S2,int col){//如果当前行被填满，即(11111111)的放满状态，递归结束，代表一种状态产生，进行状态转移if(S1==(1<<m)-1){status[S][S2]++;//方案数+1return;//退出递归}if(!check(S1,col)){//如果当前的位置不能放，向后遍历search(S,S1,S2,col+1);}else{/*if(check(S2,col)) //"I"search(S,put(S1,col,-1),put(S2,col,-1),col+1);if(check(S1,col+1)) //"--"search(S,put(S1,col,col+1),S2,col+2);*/if(check(S1,col+1)&&check(S2,col+1)) //"7"search(S,put(S1,col,col+1),put(S2,col+1,-1),col+2);if(check(S1,col+1)&&check(S2,col)) //"枪"search(S,put(S1,col,col+1),put(S2,col,-1),col+2);if(check(S2,col)&&check(S2,col-1)) //"J"search(S,put(S1,col,-1),put(S2,col-1,col),col+1);if(check(S2,col)&&check(S2,col+1)) //"L"search(S,put(S1,col,-1),put(S2,col,col+1),col+1);}}private int put(int S,int mx1,int mx2){//mx2=-1代表该位置不放格子S|=1<<mx1;//放置一个if(mx2>0) S|=1<<mx2;//放置两个return S;}private boolean check(int S,int mx){//如果在S状态下，第mx列(mx必须包含于0~m-1)的格子没有被占用if(mx>=0&&mx<m&&(S&(1<<mx))==0) return true;return false;}//矩阵相乘的方法 A=A*B,将相乘结果赋值给Aprivate void matMul(int[][] A,int[][] B){long[][] tmpMat=new long[1<<MAX_COL][1<<MAX_COL];for(int i=0;i<A.length;i++)for(int j=0;j<A[i].length;j++)for(int k=0;k<A[i].length;k++)//此处会发生int溢出，需要先用long保存取余后再转换为inttmpMat[i][j]=(tmpMat[i][j]+(long)A[i][k]*(long)B[k][j])%MOD;for(int i=0;i<A.length;i++)for(int j=0;j<A[i].length;j++)A[i][j]=(int)tmpMat[i][j];}
}

三、拓展延伸

记得之前说过的那两种理解方式吧，特意说过后者才是最标准的，为什么？如果我在原题目中增加一种积木类型，高度为2的“**I**”型积木，如图：

那么第一种理解方式会漏数方案数量，对于“I”型积木，它有一种比较特殊的摆放方式就是“平躺”放置，这样会导致它在当前行的放置操作不会直接对下一行产生影响，也就是在第一种方案中，会存在这样一种情况：1~ n-1 行完全放满，n-1行存在“平躺放置I型积木”，这样的话，只通过n-2次状态转移，导致第m行存在空缺的情况，所以再通过一次状态转移（第二种理解方式），可以只通过填充此类型的积木达到第m行也全满，图示如下：

如图所示，在第n-1行，m-3~m-2,m-1~m,分别放置了两个平躺着的I型积木，无法影响到第m行，如果使用第一种理解方式，会认为这是一种不可行的方式，所以还是得用第二种理解方式才是最正确的，要放满n行的方格图，就得有n-1次状态转移，也就是n-1次矩相乘，状态转移矩阵的n次方，代码如下：

#include <stdio.h>#define MOD 1000000007
#define MAX_COL 7#define FOR(i,start,end,step) for(int i=(start);i<(end);i+=(step))/*第10题 堆放积木@author 阡陌红尘
*/
long long n;//方格图的行数
int m;//方格图的行列数
int status[1<<MAX_COL][1<<MAX_COL]; //状态转移矩阵
int res[1<<MAX_COL][1<<MAX_COL];//保存结果状态的矩阵
int put(int S,int col1,int col2){S|=1<<col1;//改变当前行/下一行，第col1列格子的状态if(col2!=-1) S|=1<<col2;//如果传值不为-1，则也改变col2列格子的状态return S;
}
int check(int S,int col){if(col>=0&&col<m&&(S&(1<<col))==0) return 1;return 0;
}
void search(int S,int S1,int S2,int col){if(S1==(1<<m)-1){//如果当前行被放满，进行状态转移status[S][S2]++;return; }if(!check(S1,col)){//如果col列不为空 search(S,S1,S2,col+1); }else{if(check(S2,col)) //"I"search(S,put(S1,col,-1),put(S2,col,-1),col+1);if(check(S1,col+1)) //"--"search(S,put(S1,col,col+1),S2,col+2);if(check(S1,col+1)&&check(S2,col)) //"枪"search(S,put(S1,col,col+1),put(S2,col,-1),col+2);if(check(S1,col+1)&&check(S2,col+1)) //"7"search(S,put(S1,col,col+1),put(S2,col+1,-1),col+2);if(check(S2,col-1)&&check(S2,col)) //"J"search(S,put(S1,col,-1),put(S2,col-1,col),col+1);if(check(S2,col)&&check(S2,col+1)) //"L"search(S,put(S1,col,-1),put(S2,col,col+1),col+1);}
}
//矩阵相乘
void mat_mul(int A[1<<MAX_COL][1<<MAX_COL],int B[1<<MAX_COL][1<<MAX_COL]){long long matrix[1<<MAX_COL][1<<MAX_COL]={0};//临时矩阵FOR(i,0,1<<MAX_COL,1)FOR(j,0,1<<MAX_COL,1)FOR(k,0,1<<MAX_COL,1)matrix[i][j]=(matrix[i][j]+(long long)A[i][k]*(long long)B[k][j])%MOD;FOR(i,0,1<<MAX_COL,1) FOR(j,0,1<<MAX_COL,1) A[i][j]=(int)matrix[i][j];
}
int main(void){scanf("%lld %d",&n,&m);FOR(S,0,1<<m,1) search(S,S,0,0);//构建单位矩阵FOR(i,0,1<<MAX_COL,1) FOR(j,0,1<<MAX_COL,1) res[i][j]=(i==j)?1:0;   //矩阵快速幂long long t=n<<1;while(t>>=1>0){if((t&1)!=0) mat_mul(res,status);mat_mul(status,status);  }printf("%d",res[0][0]);
}

好累~,第一次写博客，难免有疏漏之处，如果有什么问题，欢迎评论区留言，嘻嘻(#^.^#)