HDU5950 Recursive sequence (矩阵快速幂加速递推) (2016ACM/ICPC亚洲赛区沈阳站 Problem C)

题目链接：传送门

题目：

Recursive sequence
Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 4173    Accepted Submission(s): 1813


Problem Description
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 < 231 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
Hint
In the first case, the third number is 85 = 2*1十2十3^4.
 In the second case, the third number is 93 = 2*1十1*10十3^4 and the fourth number is 369 = 2 * 10 十 93 十 4^4.

题目大意：

　　已知a₁ = a，a₂ = b，a_i = a_i-1 + 2a_i-2 + n⁴，求a_N，答案对2147493647取模。

　　N，a，b ＜ 2³¹。

思路：

　　因为N超级大，考虑快速幂，但是通项怎么都搞不出来。

　　又题目给的是递推式，考虑用矩阵快速幂。。。。

　　令Fn = [an-1, an]，则Fn-1 = [an-2，an-1]，但是这样an+1就算不上n4了。那就把n4加上去试试看：

　　再令Fn = [an-1，an，(n+1)4]，这样就可以推出an+1了，但是(n+1)4又不能递推。。。展开(n+1)4发现：(n+1)4 = n4 + 4n3 + 6n2 + 4n + 1，可以由n4、n3、n2、n、1递推出来。同时(n+1)3、(n+1)2、n+1、1都可以用n4、n3、n2、n、1递推出来，所以Fn和系数矩阵base就出来了：

　　Fn = [an-1，an，(n+1)4，(n+1)3，(n+1)2，n+1，1]；

$egin{bmatrix}
0 &2  &0  &0  &0  &0  &0 \
1 &1  &0  &0  &0  &0  &0 \
0 &1  &1  &0  &0  &0  &0 \
0 &0  &4  &1  &0  &0  &0 \
0 &0  &6  &3  &1  &0  &0 \
0 &0  &4  &3  &2  &1  &0 \
0 &0  &1  &1  &1  &1  &1
end{bmatrix}$

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
#define MAXN 7
const ll mod = 2147493647;

struct Matrix
{
    ll mat[MAXN][MAXN];
    Matrix() {}
    Matrix operator*(Matrix const &b)const
    {
        Matrix res;
        memset(res.mat, 0, sizeof(res.mat));
        for (int i = 0 ;i < MAXN; i++)
            for (int j = 0; j < MAXN; j++)
                for (int k = 0; k < MAXN; k++)
                    res.mat[i][j] = (res.mat[i][j]+this->mat[i][k] * b.mat[k][j])%mod;
        return res;
    }
}base;
Matrix pow_mod(Matrix base, ll n)
{
    Matrix res;
    memset(res.mat, 0, sizeof(res.mat));
    for (int i = 0; i < MAXN; i++)
        res.mat[i][i] = 1;
    while (n > 0)
    {
        if (n & 1) res = res*base;
        base = base*base;
        n >>= 1;
    }
    return res;
}//输入基础矩阵返回n次幂的矩阵;
//res.mat[0][0] 就是最终的值F(N+1)
//注意修改MAXN和mod

void init() {
    base.mat[0][0] = 0; base.mat[0][1] = 2; base.mat[0][2] = 0; base.mat[0][3] = 0; base.mat[0][4] = 0; base.mat[0][5] = 0; base.mat[0][6] = 0;
    base.mat[1][0] = 1; base.mat[1][1] = 1; base.mat[1][2] = 0; base.mat[1][3] = 0; base.mat[1][4] = 0; base.mat[1][5] = 0; base.mat[1][6] = 0;
    base.mat[2][0] = 0; base.mat[2][1] = 1; base.mat[2][2] = 1; base.mat[2][3] = 0; base.mat[2][4] = 0; base.mat[2][5] = 0; base.mat[2][6] = 0;
    base.mat[3][0] = 0; base.mat[3][1] = 0; base.mat[3][2] = 4; base.mat[3][3] = 1; base.mat[3][4] = 0; base.mat[3][5] = 0; base.mat[3][6] = 0;
    base.mat[4][0] = 0; base.mat[4][1] = 0; base.mat[4][2] = 6; base.mat[4][3] = 3; base.mat[4][4] = 1; base.mat[4][5] = 0; base.mat[4][6] = 0;
    base.mat[5][0] = 0; base.mat[5][1] = 0; base.mat[5][2] = 4; base.mat[5][3] = 3; base.mat[5][4] = 2; base.mat[5][5] = 1; base.mat[5][6] = 0;
    base.mat[6][0] = 0; base.mat[6][1] = 0; base.mat[6][2] = 1; base.mat[6][3] = 1; base.mat[6][4] = 1; base.mat[6][5] = 1; base.mat[6][6] = 1;
}

int main()
{
    init();
    int T;
    cin >> T;
    while (T--) {
        ll n, a, b;
        scanf("%lld%lld%lld", &n, &a, &b);
        Matrix F2;
        memset(F2.mat, 0, sizeof F2.mat);
        F2.mat[0][0] = a; F2.mat[0][1] = b; F2.mat[0][2] = 3*3*3*3; F2.mat[0][3] = 3*3*3; F2.mat[0][4] = 3*3; F2.mat[0][5] = 3; F2.mat[0][6] = 1;
        Matrix FN = F2*pow_mod(base, n-2);
        cout << FN.mat[0][1] << endl;
    }
    return 0;
}