hdu5015 233 Matrix 矩阵快速幂 矩阵构造方法

233 Matrix

题目的意思是有一个叫做233 Matrix的矩阵。

给定第一行元素a(0, 1) = 233, a(0, 2) = 2333, a(0, 3) = 23333, ..., a(0, n) = 10 * a(0, n - 1) + 3, (n >= 2)。

第一列元素a(1, 0), a(2, 0), a(3, 0), ..., a(n, 0)，和一个递推式子a(i ,j) = a(i - 1, j) + a(i, j - 1), (i,j ≠ 0)，让你求出a(n, m) mod 10000007的值。

其中，a(0, 0)应该等于0，n和m及第一列元素由输入给定。数据范围为n ≤ 10,m ≤ 10⁹

因为数据量很大，暴力求解肯定是不行的。必定超时，虽然给定的时间是5秒钟。

那接下来就是推公式了，而我最终也真的推出来个公式，但很遗憾，它只有当n和m相等的时候才成立，可是n最大才为10，这种思路也就被否定了。

既然给定的是一个矩阵的形式，难道可以用矩阵快速幂搞出来？可以试试。注意一点，如果能构造出来快速幂的话，也只有按列来推。

对于给定的递推式a(i ,j) = a(i - 1, j) + a(i, j - 1)，将它展开则a(i, j) = a(i - 1, j) + a(i, j - 1) = a(i - 2, j) + a(i - 1, j - 1) + a(i, j - 1) = a(i - 3, j) + a(i - 2, j - 1) + a(i - 1, j - 1) + a(i, j - 1) = ...

很快就可以发现规律了得到：a(i, j) = a(1, j - 1) + a(2, j - 1) + ... + a(i, j - 1) + a(0, j)。

如下图：

图片橙色部分等于蓝色部分的和，也就是上边的公式。

我们的目的是由当前蓝色的部分，推出与它紧邻的右边下一列蓝色的部分。

构造如下矩阵，我们希望找出一个n+1阶的方阵，使得如下矩阵式成立。

根据上述推导的公式，可得A矩阵的第一行为[1 0 ... 0 1]，第二行为[1 1 0 ... 0 1]，第三行为[1 1 1 0 ... 0 1]，。。。，第n行为[1 1 1 ... 1 1]。

遗憾的是第n+1行无论怎么取值也无法是矩阵式成立。

因为a(0, m + 1） = a（0，m）* 10 + 3，那么显然n+1行的矩阵不能满足要求，可以将矩阵变成n+2行最后一行为3。则矩阵就变成：

相应的A矩阵第一行为[1 0 ... 0 1 0]，第二行为[1 1 0 ... 0 1 0]，第三行为[1 1 1 0 ... 0 1 0]，。。。，第n行为[1 1 1 ... 1 1 0]，第n+1行为[0 0 0 0 .... 0 10 1]，

第n+2行为[0 0 0 0 ... 0 0 1]。则矩阵A的形式如下：

由于矩阵乘法的结合律，可得：

这样就可以对矩阵A^m使用矩阵快速幂了。

代码如下（代码中构造的矩阵可能和上述矩阵不符）：

/*************************************************************************    > File Name: matrix.cpp    > Author: gwq    > Mail: [e-mail]     > Created Time: 2014年09月18日 星期四 08时32分40秒 ************************************************************************/#include <cmath>#include <ctime>#include <cctype>#include <climits>#include <cstdio>#include <cstdlib>#include <cstring>#include <map>#include <set>#include <queue>#include <stack>#include <vector>#include <iostream>#include <algorithm>#define INF (INT_MAX / 10)#define SQR(x) ((x) * (x))#define rep(i, n) for (int i = 0; i < (n); ++i)#define repf(i, a, b) for (int i = (a); i <= (b); ++i)#define repd(i, a, b) for (int i = (a); i >= (b); --i)#define clr(arr, val) memset(arr, val, sizeof(arr))#define pb push_back#define sz(a) ((int)(a).size())#define mid(x, y) ((x + y) / 2)using namespace std;typedef set<int> si;typedef vector<int> vi;typedef map<int, int> mii;typedef long long ll;#define N 20#define MOD 10000007typedef struct Node {    int r, c;    int mat[N][N];    Node(int r, int c)    {        //构造相伴矩阵        this->r = r;        this->c = c;        clr(this->mat, 0);        for (int i = 0; i < r - 1; ++i) {            for (int j = 0; j < i + 1; ++j) {                this->mat[i][j] = 1;            }        }        this->mat[0][0] = 10;        this->mat[0][c - 1] = 1;        this->mat[r - 1][c - 1] = 1;    }    void show(void)    {        for (int i = 0; i < r; ++i) {            for (int j = 0; j < c; ++j) {                printf("%d%c", this->mat[i][j], (j == c - 1) ? '\n' : ' ');            }        }    }    Node operator *(const Node& u)    {        Node res(this->r, u.c);        clr(res.mat, 0);        for (int i = 0; i < res.r; ++i) {            for (int j = 0; j < res.c; ++j) {                for (int k = 0; k < this->c; ++k) {                    res.mat[i][j] = ((long long)this->mat[i][k] * u.mat[k][j] + (long long)res.mat[i][j]) % MOD;                }            }        }        return res;    }}Node;int main(int argc, char *argv[]){    int n, m;    while (scanf("%d%d", &n, &m) != EOF) {        Node res(n + 2, 1);        for (int i = 1; i <= n; ++i) {            scanf("%d", &res.mat[i][0]);        }        res.mat[0][0] = 233;        res.mat[n + 1][0] = 3;        Node A(n + 2, n + 2);        Node B(n + 2, n + 2);        clr(B.mat, 0);        for (int i = 0; i < n + 2; ++i) {            B.mat[i][i] = 1;        }        while (m != 0) {            if ((m & 1) == 1) {                B = B * A;            }            A = A * A;            m >>= 1;        }        res = B * res;        printf("%d\n", res.mat[n][0]);    }    return 0;}