HDU - 5950 Recursive sequence（矩阵快速幂）

Recursive sequence
Farmer John likes to play mathematics games with his N cows. Recently, they are attracted by recursive sequences. In each turn, the cows would stand in a line, while John writes two positive numbers a and b on a blackboard. And then, the cows would say their identity number one by one. The first cow says the first number a and the second says the second number b. After that, the i-th cow says the sum of twice the (i-2)-th number, the (i-1)-th number, and i4. Now, you need to write a program to calculate the number of the N-th cow in order to check if John’s cows can make it right.
Input
The first line of input contains an integer t, the number of test cases. t test cases follow.
Each case contains only one line with three numbers N, a and b where N,a,b < 2³¹ as described above.
Output
For each test case, output the number of the N-th cow. This number might be very large, so you need to output it modulo 2147493647.
Sample Input
2
3 1 2
4 1 10
Sample Output
85
369

题意：题目给出了一个递推式f(n)=f(n-1)+2*f(n-2)+n⁴，然后给出了第一项a，和第二项b，要求第n项的值。

思路：看到递推式，然后范围较大的就想到用矩阵快速幂，这题的难点就是如何求出这个关系矩阵，因为这里有一个n⁴，要找到 n⁴ 和 (n-1)⁴ 的关系,以4⁴为例：

4⁴=C(4,0)*3⁴+C(4,1)*3³+C(4,2)*3²+C(4,3)*3¹+C(4,4)*3⁰

（注意这里的C(n,k)指的是组合数！）
知道了这个是怎么转换过来了以后，就可以推出关系矩阵:

| 1 2 1 4 6 4 1 |
| 1 0 0 0 0 0 0 |
| 0 0 1 4 6 4 1 |
| 0 0 0 1 3 3 1 |
| 0 0 0 0 1 2 1 |
| 0 0 0 0 0 1 1 |
| 0 0 0 0 0 0 1 |

代码：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <queue>
#include <stack>
#include <vector>
#include <set>
#include <map>
#include <list>
#include <deque>
#include <iomanip>
#include <iostream>
#include <algorithm>
using namespace std;

#define IOS                      \
    ios::sync_with_stdio(false); \
    cin.tie(0);                  \
    cout.tie(0)
#define lowbit(a) (a & (-a))
#define ll long long
const int inf = 0x3f3f3f3f;
const int maxn = 2e5 + 10;
const int minn = 1e3 + 10;
const double PI = acos(-1.0);
#define mod 2147493647
#define eps 1e-10
/*------------------------------------------------------*/

int N = 7;

struct Matrix
{
    ll mp[7][7];
    Matrix()
    {
        memset(mp, 0, sizeof(mp));
    }
    void init()
    {
        memset(mp, 0, sizeof(mp));
        for (int i = 0; i < N; i++)
            mp[i][i] = 1;
    }
};

Matrix mul(Matrix a, Matrix b)
{
    Matrix ans;
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            for (int k = 0; k < N; k++)
                ans.mp[i][j] = ((ans.mp[i][j] % mod) + ((a.mp[i][k] % mod) * (b.mp[k][j] % mod)) % mod) % mod;
    return ans;
}

Matrix mat_pow(Matrix a, ll b)
{
    Matrix ans;
    ans.init();
    while (b)
    {
        if (b & 1)
            ans = mul(ans, a);
        a = mul(a, a);
        b >>= 1;
    }
    return ans;
}

int main()
{
    int t;
    scanf("%d", &t);
    while (t--)
    {
        ll n, a, b;
        scanf("%lld %lld %lld", &n, &a, &b);
        if (n == 1)
            printf("%lld\n", a % mod);
        else if (n == 2)
            printf("%lld\n", b % mod);
        else
        {
            Matrix cnt;
            cnt.init();
            cnt.mp[0][0] = 1, cnt.mp[0][1] = 2, cnt.mp[0][2] = 1, cnt.mp[0][3] = 4, cnt.mp[0][4] = 6, cnt.mp[0][5] = 4, cnt.mp[0][6] = 1;
            cnt.mp[1][0] = 1, cnt.mp[1][1] = 0, cnt.mp[1][2] = 0, cnt.mp[1][3] = 0, cnt.mp[1][4] = 0, cnt.mp[1][5] = 0, cnt.mp[1][6] = 0;
            cnt.mp[2][0] = 0, cnt.mp[2][1] = 0, cnt.mp[2][2] = 1, cnt.mp[2][3] = 4, cnt.mp[2][4] = 6, cnt.mp[2][5] = 4, cnt.mp[2][6] = 1;
            cnt.mp[3][0] = 0, cnt.mp[3][1] = 0, cnt.mp[3][2] = 0, cnt.mp[3][3] = 1, cnt.mp[3][4] = 3, cnt.mp[3][5] = 3, cnt.mp[3][6] = 1;
            cnt.mp[4][0] = 0, cnt.mp[4][1] = 0, cnt.mp[4][2] = 0, cnt.mp[4][3] = 0, cnt.mp[4][4] = 1, cnt.mp[4][5] = 2, cnt.mp[4][6] = 1;
            cnt.mp[5][0] = 0, cnt.mp[5][1] = 0, cnt.mp[5][2] = 0, cnt.mp[5][3] = 0, cnt.mp[5][4] = 0, cnt.mp[5][5] = 1, cnt.mp[5][6] = 1;
            cnt.mp[6][0] = 0, cnt.mp[6][1] = 0, cnt.mp[6][2] = 0, cnt.mp[6][3] = 0, cnt.mp[6][4] = 0, cnt.mp[6][5] = 0, cnt.mp[6][6] = 1;
            cnt = mat_pow(cnt, n - 2);
            printf("%lld\n", (cnt.mp[0][0] * b % mod + cnt.mp[0][1] * a % mod + cnt.mp[0][2] * 16 % mod + cnt.mp[0][3] * 8 % mod + cnt.mp[0][4] * 4 % mod + cnt.mp[0][5] * 2 % mod + cnt.mp[0][6] % mod) % mod);
        }
    }
    return 0;
}