POJ3420 递推/状态压缩DP +矩阵幂加速处理

最新推荐文章于 2022-03-02 10:02:44 发布

focus_best

最新推荐文章于 2022-03-02 10:02:44 发布

阅读量2.3k

点赞数 2

分类专栏： ACM--题解汇总 ACM--递推 ACM--状态压缩DP 文章标签： ACM

本文链接：https://blog.csdn.net/u013480600/article/details/19644847

版权

ACM--题解汇总同时被 3 个专栏收录

685 篇文章 11 订阅

订阅专栏

★★★

112 篇文章 0 订阅

订阅专栏

ACM--状态压缩DP

17 篇文章 1 订阅

订阅专栏

POJ3420 Quad Tiling

现在我们有一个4*n（1<=n<=10^9）矩阵，要求用1*2的骨牌填充，求填充的方法总数。由于结果较大，所以输出结果%m即可。

输入：包含多组实例，每组为n和m（0<m<=10^5）。当输入两个0的时候表示输入结束。

输出：输出摆放总数对m求余的结果。

分析：做本题首先要知道它的基础题型。见：http://blog.csdn.net/u013480600/article/details/19612171

这道题由于数据规模很大，且是4*n类型，所以不能用轮廓线DP来解。

解法一：递推+矩阵加速

首先尝试用递推解法来解。我们首先用手模拟画出4*2,4*3,4*4,4*5以及4*6的不可分割的矩阵。令a[n]表4*n规模的矩阵的摆放骨牌总数，b[n]表示4*n规模的不可分割的矩阵个数

我们会发现a[1]=b[1]=1，b[2]=4，且n>=3时有以下：

当n为奇数时：b[n]=2

当n为偶数时：b[n]=3

其中

（注意第4个矩形也是不可分割的，不可分割只是说4*n的n行分割，列永不分割的）

分别为4个不可分割的4*2矩阵构成。

所以按照这样a[n]表4*n规模的矩阵的摆放骨牌总数。

a[n]=a[n-1]*b[1]+a[n-2]*b[2]+…a[1]*b[n-1]+a[0]*b[n]

由上公式可得：

a[n]=a[n-1]+ a[n-2]*4+a[n-3]*2+a[n-4]*3+…a[0]*b[n]

a[n-1]= a[n-2]+a[n-3]*4+a[n-4]*2+a[n-5]*3+…a[0]*b[n-1]

以上两式子相加得：

a[n] = 5a[n-2]+6a[n-3]+5*（a[n-4]+a[n-5]+…a[0]）

且可得a[n-1] = 5a[n-3]+6a[n-4]+5*（a[n-5]+a[n-6]+…a[0]）

以上两式子相减得到：

a[n]=a[n-1]+5a[n-2]+a[n-3]-a[n-4] （n>=4且a[0]=a[1]=1，a[2]=5，a[3]=11）用轮廓线DP计算该题验证了数据，可知这个公式是对的。

由于n的规模为10^9所以如果直接递推速度也很慢。所以由上述红字公式可知如下矩阵递推关系：

* =

令矩阵A = 则要求a[n]的值只需要用A^(n-3) * 然后取最后一行的值即可。所以在此引入矩阵类模板，并且自己利用矩阵模板的*运算符写一个方阵的幂运算函数进行矩阵运算即可。

AC代码:0ms

#include<cstdio>
#include<cstring>
using namespace std;
const int MAXN=10;
const int MAXM=10;

int r,mod;
//矩阵类模板开始**********************************************
struct Matrix
{
    int n,m;
    int a[MAXN][MAXM];
    void clear()
    {
        n=m=0;
        memset(a,0,sizeof(a));
    }
    Matrix operator +(const Matrix &b)const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=m;
        for(int i=0; i<n; ++i)
            for(int j=0; j<m; j++)
                tmp.a[i][j] = a[i][j]+b.a[i][j];
        return tmp;
    }
    Matrix operator -(const Matrix &b)const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=m;
        for(int i=0; i<n; ++i)
            for(int j=0; j<m; j++)
                tmp.a[i][j] = a[i][j]-b.a[i][j];
        return tmp;
    }
    Matrix operator *(const Matrix &b)const
    {
        Matrix tmp;
        tmp.clear();
        tmp.n=n;
        tmp.m=b.m;
        for(int i=0; i<n; ++i)
            for(int j=0; j<b.m; j++)
                for(int k=0; k<m; k++)
                {
                    tmp.a[i][j] += a[i][k]*b.a[k][j];
                    tmp.a[i][j] %=mod;
                }
        return tmp;
    }
    Matrix get(int x)//方阵幂运算，*this必须为方阵且x>=0且为整数，得到（*this）的x次幂方阵
    {
        Matrix E;
        E.clear();
        E.n=E.m=n;
        for(int i=0; i<n; i++)
        {
            E.a[i][i]=1;
        }
        if(x==0)return E;
        else if(x==1)return *this;
        Matrix tmp = get(x/2);
        tmp = tmp*tmp;
        if(x%2)tmp = tmp*(*this);
        return tmp;
    }
};
//矩阵类模板结束**********************************************
int main()
{
    while(scanf("%d %d",&r,&mod)==2)
    {
        if(r==0&&mod==0)break;
        int a[]= {1,1,5,11};
        if(r<=3)
        {
            printf("%d\n",a[r]%mod);
            continue;
        }
        Matrix A;
        A.clear();
        A.n=A.m=4;
        A.a[0][1]=A.a[1][2]=A.a[2][3]=1;
        A.a[3][0]=-1;
        A.a[3][1]=1;
        A.a[3][2]=5;
        A.a[3][3]=1;
        A = A.get(r-3);//得到A的r-3次幂
        Matrix x;
        x.clear();
        x.n=4;
        x.m=1;
        x.a[0][0]=1;
        x.a[1][0]=1;
        x.a[2][0]=5;
        x.a[3][0]=11;//x此时是那个列向量（用于右乘A方阵的）
        x = A*x;//乘法顺序不要搞错
        printf("%d\n",(x.a[3][0]+mod)%mod);//这里可能为负值，所以需要+mod再求余
    }
    return 0;
}

解法二：状态压缩+矩阵加速

现在我们看下能否用基于状态压缩的DP来算这题呢？用pre表示一个兼容对的前一行的二进制形式的十进制值，用now表示一个兼容对的后一行的二进制形式的十进制值。由于该题是列固定为4，所以pre和now的值最多有16种，设A[pre][now]=1表示兼容对前行pre和后行now的二进制形式兼容有1种情况（该值不是1就是0）。则令B=A*A，则

B[2][3]=A[2][0]*A[0][3]

+A[2][1]*A[1][3]

+A[2][2]*A[2][3]

+A[2][3]*A[3][3]

即（A*A）[2][3]表示pre=2与now=3之间还间隔了一行，并且pre与now远程兼容有多少种方式。（pre与now的远程距离为1）.由于目标矩阵有n行，行号从1到n，则我们假想目标矩阵的上面还有1行是第0行，并且二进制形式为全1序列。所以我们要求的是第0行和第n行距离为n-1的远程兼容有多少种方法。

即求得就是（A^n）[0][15]的值，该值还要对mod求余。

求矩阵幂的时候注意，本来矩阵需要16*16的大小，结果你申请了100*100的大小，由于矩阵幂递归求，直接导致你如果求的矩阵幂太大比如A的100000次方时，程序内存不够，直接崩溃。但是如果用的是迭代求矩阵幂的方式可以随便开矩阵大小为100*100等。

递归方式时间为47ms,迭代方式时间为32ms。

AC代码：

//状态压缩DP+矩阵加速解法
#include<cstdio>
#include<cstring>
using namespace std;
const int MAXN=100;
const int MAXM=100;

int r,mod;
//矩阵类模板开始**********************************************
struct Matrix
{
    int n,m;
    int a[MAXN][MAXM];
    void clear()
    {
        n=m=0;
        memset(a,0,sizeof(a));
    }
    Matrix operator +(const Matrix &b)const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=m;
        for(int i=0; i<n; ++i)
            for(int j=0; j<m; j++)
                tmp.a[i][j] = a[i][j]+b.a[i][j];
        return tmp;
    }
    Matrix operator -(const Matrix &b)const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=m;
        for(int i=0; i<n; ++i)
            for(int j=0; j<m; j++)
                tmp.a[i][j] = a[i][j]-b.a[i][j];
        return tmp;
    }
    Matrix operator *(const Matrix &b)const
    {
        Matrix tmp;
        tmp.clear();
        tmp.n=n;
        tmp.m=b.m;
        for(int i=0; i<n; ++i)
            for(int j=0; j<b.m; j++)
                for(int k=0; k<m; k++)
                {
                    tmp.a[i][j] += a[i][k]*b.a[k][j];
                    tmp.a[i][j] %=mod;
                }
        return tmp;
    }
    Matrix get(int x)//方阵幂运算，*this必须为方阵且x>=0且为整数，得到（*this）的x次幂方阵
    {
        Matrix E;
        E.clear();
        E.n=E.m=n;
        for(int i=0; i<n; i++)
        {
            E.a[i][i]=1;
        }
        if(x==0)return E;
        else if(x==1)return *this;
        Matrix tmp = E,b=*this;
        while(x)
        {
            if(x&1)tmp = tmp*(b);
            b = b*b;
            x>>=1;
        }
        return tmp;
    }
};
//矩阵类模板结束**********************************************
Matrix path;
void dfs(int c,int pre,int now)
{
    if(c>4)return ;
    else if(c==4)
    {
        path.a[pre][now]++;
        return ;
    }
    dfs(c+1,(pre<<1)|1,(now<<1) );//不放
    dfs(c+1,(pre<<1),(now<<1)|1 );//上方
    dfs(c+2,(pre<<2)|3,(now<<2)|3 );//右方
}
int main()
{
    path.clear();
    path.n=16;
    path.m=16;
    dfs(0,0,0);
    while(scanf("%d %d",&r,&mod)==2)
    {
        Matrix A;
        A.clear();
        A.n=A.m=16;
        A=path;
        if(r==0&&mod==0)break;
        A=A.get(r);
        printf("%d\n",(A.a[15][15]+mod)%mod);//这里可能为负值，所以需要+mod再求余
    }
    return 0;
}