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

http://acm.hdu.edu.cn/showproblem.php?pid=5950

题意：给出 a，b，n，递推出 f(n) = f(n-1) + f(n-2) * 2 + n ^ 4. f(1) = a, f(2) = b.

思路：在比赛时候知道是矩阵快速幂，可是推不出矩阵.那个n^4不知道怎么解决。结束后问其他人才知道要构造一个7 * 7的矩阵，而不是3 * 3的..

转自：http://blog.csdn.net/spring371327/article/details/52973534

#include <cstdio>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <string>
#include <cmath>
#include <queue>
#include <vector>
using namespace std;
#define N 1010
#define INF 0x3f3f3f3f
#define MOD 2147493647
typedef long long LL;

struct matrix
{
    LL a[7][7];

    void init() {
        memset(a, 0, sizeof(a));
        for(int i = 0; i < 7; i++) a[i][i] = 1;
    }

    matrix operator * (matrix b) {
        matrix ans;
        LL tmp;
        for(int i = 0; i < 7; i++) {
            for(int j = 0; j < 7; j++) {
                ans.a[i][j] = 0;
                for(int k = 0; k < 7; k++) {
                    tmp = a[i][k] * b.a[k][j] % MOD;
                    ans.a[i][j] = (ans.a[i][j] + tmp % MOD) % MOD;
                }
            }
        }
        return ans;
    }
};

matrix q_pow(matrix a, LL b)
{
    matrix ans;
    ans.init();
    while(b) {
        if(b & 1) ans = ans * a;
        b >>= 1;
        a = a * a;
    }
    return ans;
}

int main()
{
    matrix mo;
    memset(mo.a, 0, sizeof(mo.a));
    mo.a[0][1] = 1;
    mo.a[1][0] = 2; mo.a[1][1] = 1, mo.a[1][2] = 1, mo.a[1][3] = 4, mo.a[1][4] = 6, mo.a[1][5] = 4, mo.a[1][6] = 1;
    mo.a[2][2] = 1, mo.a[2][3] = 4, mo.a[2][4] = 6, mo.a[2][5] = 4, mo.a[2][6] = 1;
    mo.a[3][3] = 1, mo.a[3][4] = 3, mo.a[3][5] = 3, mo.a[3][6] = 1;
    mo.a[4][4] = 1, mo.a[4][5] = 2, mo.a[4][6] = 1;
    mo.a[5][5] = 1, mo.a[5][6] = 1;
    mo.a[6][6] = 1;
    int t;
    scanf("%d", &t);
    while(t--) {
        long long n, a, b;
        scanf("%I64d%I64d%I64d", &n, &a, &b);
        if(n == 1) printf("%I64d
", a);
        else if(n == 2) printf("%I64d
", b);
        else {
            matrix ans = q_pow(mo, n - 2);
            LL sum = 0;
            sum = (sum + ans.a[1][0] * a) % MOD;
            sum = (sum + ans.a[1][1] * b) % MOD;
            sum = (sum + ans.a[1][2] * 16) % MOD;
            sum = (sum + ans.a[1][3] * 8) % MOD;
            sum = (sum + ans.a[1][4] * 4) % MOD;
            sum = (sum + ans.a[1][5] * 2) % MOD;
            sum = (sum + ans.a[1][6]) % MOD;
            printf("%I64d
", sum);
        }
    }
    return 0;
}