hdu5950 矩阵dp

题目链接：点击打开链接

题意：

给定n， 第一个数 a，第二个数 b；

公式：Ai = A(i - 2) * 2 + A(i - 1) + i ^ 4；

求第 n 个数；

其中 n 的范围是 int；

理解：

由公式来看，该题很像斐波那契；

但是 n 的范围太大；

偶然想到斐波那契也可以用矩阵dp来推；

所以就应该用这个方法；

从公式可以看出，在不考虑 i ^ 4 的情况下：

A(i - 1)   [0 1]   A(i - 2)
         = [   ] *
A(i)       [2 1]   A(i - 1)

然后在根据：

A ^ k = C(0, k) * (A - 1) ^ k + C(1, k) * (A - 1) ^ (k - 1) + ... + C(k - 1, k) * (A - 1) ^ 1 + C(k, k) * (A - 1) ^ 0

其中C(0, k) 表示组合数；

这样可推出：

A(i - 1)      [0    1    0       0       0       0       0      ]   A(i - 2)
A(i)          [2    1    1       0       0       0       0      ]   A(i - 1)
(i + 1) ^ 4   [0    0    C(0, 4) C(1, 4) C(2, 4) C(3, 4) C(4, 4)]   i ^ 4
(i + 1) ^ 3 = [0    0    0       C(0, 3) C(1, 3) C(2, 3) C(3, 3)] * i ^ 3
(i + 1) ^ 2   [0    0    0       0       C(0, 2) C(1, 2) C(2, 2)]   i ^ 2
(i + 1) ^ 1   [0    0    0       0       0       C(0, 1) C(1, 1)]   i ^ 1
(i + 1) ^ 0   [0    0    0       0       0       0       C(0, 0)]   i ^ 0

则递乘的矩阵矩阵 M 为：

[0 1 0 0 0 0]
[2 1 1 0 0 0]
[0 0 1 6 4 1]
[0 0 0 3 3 1]
[0 0 0 1 2 1]
[0 0 0 0 1 1]
[0 0 0 0 0 1]

因此就可以通过以下快速幂的方式算出结果了

A(n - 1)                    a
A(n)                        b
(n + 1) ^ 4                 3 ^ 4
(n + 1) ^ 3 = M ^ (n - 2) * 3 ^ 3
(n + 1) ^ 2                 3 ^ 2
(n + 1) ^ 1                 3 ^ 1
(n + 1) ^ 0                 3 ^ 0

代码如下：

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;

const int MAXN = 1e6 + 10;
LL MOD = 2147493647UL;

typedef pair<int, int> PII;

#define X first
#define Y second

typedef vector<LL> vec;
typedef vector<vec> mat;

mat mul(mat &A, mat &B) {
    mat C(A.size(), vec(B[0].size()));
    for (int i = 0; i < A.size(); ++i) {
        for (int k = 0; k < B.size(); ++k) {
            for (int j = 0; j < B[0].size(); ++j) {
                C[i][j] = (C[i][j] + A[i][k] * B[k][j] % MOD) % MOD;
            }
        }
    }
    return C;
}

mat poww(mat A, LL n) {
    mat B(A.size(), vec(A.size()));
    for (int i = 0; i < A.size(); ++i) {
        B[i][i] = 1;
    }
    while (n > 0) {
        if (n & 1) {
            B = mul(B, A);
        }
        A = mul(A, A);
        n /= 2;
    }
    return B;
}

int main() {
    int t;
    cin >> t;
    while (t--) {
        LL n, a, b;
        cin >> n >> a >> b;
        mat A(7, vec(7));
        A[0][0] = 0, A[0][1] = 1, A[0][2] = 0, A[0][3] = 0, A[0][4] = 0, A[0][5] = 0, A[0][6] = 0;
        A[1][0] = 2, A[1][1] = 1, A[1][2] = 1, A[1][3] = 0, A[1][4] = 0, A[1][5] = 0, A[1][6] = 0;
        A[2][0] = 0, A[2][1] = 0, A[2][2] = 1, A[2][3] = 4, A[2][4] = 6, A[2][5] = 4, A[2][6] = 1;
        A[3][0] = 0, A[3][1] = 0, A[3][2] = 0, A[3][3] = 1, A[3][4] = 3, A[3][5] = 3, A[3][6] = 1;
        A[4][0] = 0, A[4][1] = 0, A[4][2] = 0, A[4][3] = 0, A[4][4] = 1, A[4][5] = 2, A[4][6] = 1;
        A[5][0] = 0, A[5][1] = 0, A[5][2] = 0, A[5][3] = 0, A[5][4] = 0, A[5][5] = 1, A[5][6] = 1;
        A[6][0] = 0, A[6][1] = 0, A[6][2] = 0, A[6][3] = 0, A[6][4] = 0, A[6][5] = 0, A[6][6] = 1;
        A = poww(A, n - 2);
        mat ans(7, vec(1));
        ans[0][0] = a, ans[1][0] = b, ans[2][0] = 3 * 3 * 3 * 3, ans[3][0] = 3 * 3 * 3;
        ans[4][0] = 3 * 3, ans[5][0] = 3, ans[6][0] = 1;
        ans = mul(A, ans);
        cout << ans[1][0] << endl;
    }

    return 0;
}