【计算理论与算法分析设计】 6. Tri Tiling (HDU-1143) (POJ-2663)

Tri Tiling

Description

In how many ways can you tile a 3xn rectangle with 2x1 dominoes?
Here is a sample tiling of a 3x12 rectangle.

Input

Input consists of several test cases followed by a line containing -1. Each test case is a line containing an integer 0 <= n <= 30.

Output

For each test case, output one integer number giving the number of possible tilings.

Sample Input

2
8
12
-1

Sample Output

3
153
2131

Source

Waterloo local 2005.09.24

思路分析

这个题其实是在《组合数学》课程当中学到的，在第1章的习题和第7章“递推关系与生成函数”一节都有提到。再加上之前做过类似的题目，就想到这应该是有一个递推关系的。首先，n为奇数的时候，方案数显然为0；n为偶数时，设g_n为可能的摆放方法数，则g_n满足递推关系：$g_{n}=4g_{n-2}-g_{n-4}$gn​=4gn−2​−gn−4​。下面的这个图比较显明地展现了递推关系的由来：

本题数据范围为0 <= n <= 30，在int范围内可以完成递推求解。但是，在《组合数学》课程中，老师把本题的数据范围扩大到了10⁹，要求输出结果模9973的余数。将数据范围扩大后，一项一项递推将会TLE，这就需要运用矩阵快速幂的思想去解决问题了。

先来谈一谈什么是快速幂。所谓快速幂，即是快速地计算a^b的一种算法，其原理是利用了b的二进制表示形式。比如我们要计算a¹³，13的二进制是1101，我们只需要计算a⁸，a⁴和a即可。下面的代码实现了这一想法：

int quickpow(int a,int b) //快速幂模板
{
    ans=1;
    while(b) //如果b不是0
    {
        if(b&1) ans*=a; //b&1得到的是b的二进制最低位，如果最低位是1，ans就乘以a
        a*=a; //a变成原来的平方
        b>>=1; //将b右移1位，考虑二进制的下一位
    }
    return ans;
}

a¹³的例子计算过程如下：

1. ans=1，b=13，b不是0，进入while循环；
2. b&1为真，因为1101最后一位是1，ans=a，a=a²，b>>=1，b变成6（110）；
3. b&1为假，a=a⁴，b>>=1，b变成3（11）；
4. b&1为真，ans=a*a⁴=a⁵，a=a⁸，b>>=1，b变成1（1）；
5. b&1为真，ans=a⁵*a⁸=a¹³，a=a¹⁶，b>>=1，b变成0；
6. b为0，退出while循环，返回ans=a¹³。

与常规的乘法进行比较，快速幂的时间复杂度为O(logn)，大大提高了计算a^b的效率。

而矩阵快速幂，即是把快速幂模板中的ans和a都换成矩阵，将数字的乘法变成矩阵乘法。那么为什么本题可以用矩阵快速幂来求解呢？就是因为g_n有递推关系。考虑这一递推关系，我们首先构造得到：
$\begin{bmatrix} g_{4}\\ g_{2} \end{bmatrix}= \begin{bmatrix} 4 & -1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} g_{2}\\ g_{0} \end{bmatrix}$[g4​g2​​]=[41​−10​][g2​g0​​]
为什么要写成这样的形式呢？因为这样才能进一步递推下去。我们有：
$\begin{bmatrix} g_{6}\\ g_{4} \end{bmatrix}= \begin{bmatrix} 4 & -1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} g_{4}\\ g_{2} \end{bmatrix}= \begin{bmatrix} 4 & -1\\ 1 & 0 \end{bmatrix}^{2} \begin{bmatrix} g_{2}\\ g_{0} \end{bmatrix}$[g6​g4​​]=[41​−10​][g4​g2​​]=[41​−10​]2[g2​g0​​]
所以，
$\begin{bmatrix} g_{n}\\ g_{n-2} \end{bmatrix}= \begin{bmatrix} 4 & -1\\ 1 & 0 \end{bmatrix}^{\frac{n}{2}-1} \begin{bmatrix} g_{2}\\ g_{0} \end{bmatrix}$[gn​gn−2​​]=[41​−10​]2n​−1[g2​g0​​]
我们可以使用矩阵快速幂的方法来计算矩阵的$\frac{n}{2}-1$2n​−1次幂。初始化如下：
$a=\begin{bmatrix} 4 & -1\\ 1 & 0 \end{bmatrix},\quad ans=\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}$a=[41​−10​],ans=[10​01​]
在快速幂中，ans初始化为1，那么类比一下，在矩阵快速幂中，ans初始化为单位阵。计算完成后，再乘以由g₂=3，g₀=1组成的矩阵，便可求得g_n了。

AC代码

(1) 原始题目

#include <cstdio>
int a[35];
int solve(int n)
{
    if (n==0) return 1;
    else if (n==2) return 3; //递归边界
    else if (a[n]!=0) return a[n]; //如果已经计算过了，直接返回，避免重复计算浪费时间
    else
    {
        a[n]=4*solve(n-2)-solve(n-4); //递归计算
        return a[n];
    }
}
int main()
{
    int n;
    while (scanf ("%d",&n))
    {
        if (n==-1) break;
        else if (n%2==1) printf ("0\n");
        else printf ("%d\n",solve(n));
    }
    return 0;
}

(2) 增大数据范围后的题目，使用矩阵快速幂

#include <cstdio>
struct matrix //定义矩阵类，方便后边调用函数进行矩阵乘法
{
    long long x[2][2];
}a,b;
const int mod=9973;
void init() //初始化
{
    a.x[0][0]=4;
    a.x[0][1]=-1;
    a.x[1][0]=1;
    a.x[1][1]=0;
    b.x[0][0]=1;
    b.x[0][1]=0;
    b.x[1][0]=0;
    b.x[1][1]=1;
}
matrix multiply(matrix b,matrix a) //矩阵乘法函数
{
    matrix c;
    for (int i=0;i<2;i++)
    {
        for (int j=0;j<2;j++)
        {
            c.x[i][j]=0;
            for (int k=0;k<2;k++)
            {
                c.x[i][j]+=b.x[i][k]*a.x[k][j];
            }
        }
    }
    for (int i=0;i<2;i++) //一定要注意取模，但取模是很费时间的，尤其是大数不宜取模过多，否则很有可能会超时
    {
        for (int j=0;j<2;j++)
        {
            c.x[i][j]%=mod;
        }
    }
    return c;
}
int main()
{
    int n;
    while (scanf ("%d",&n))
    {
        if (n==-1) break;
        else if (n%2==1) printf ("0\n");
        else if (n==0) printf ("1\n");
        else if (n==2) printf ("3\n"); //特判直接输出
        else
        {
            init();
            int cnt=n/2-1;
            while (cnt) //矩阵快速幂
            {
                if (cnt&1) b=multiply(b,a);
                a=multiply(a,a);
                cnt>>=1;
            }
            long long ans=((b.x[0][0]*3)%mod+(b.x[0][1])%mod)%mod;
            if (ans<0) ans+=mod; //这里可能结果会是负数，要判断一下
            printf ("%lld\n",ans);
        }
    }
    return 0;
}

后来学过面向对象之后，发现矩阵快速幂其实可以写一个矩阵类，然后重载乘法运算符，这样就可以在矩阵对象之间正常地使用乘号进行矩阵乘法运算了。