HDU 4686 Arc of Dream （矩阵快速幂）

题意：An Arc of Dream is a curve defined by following function:

where
a 0 = A0
a i = a i-1 * AX+AY
b 0 = B0
b i = b i-1 * BX+BY
给出n，A0，AX，AY，B0，BX，BY，What is the value of AoD(N) modulo 1,000,000,007?

题解：矩阵快速幂
数据范围很恶心，先要对原数据取模，对于两个数的乘法也要马上取模。
题目给出的转移关系很明确，由于要求和，可得转移矩阵为
$$T=\left[\begin{array}{ccccc}AX& 0& AY& 0& 0\\ 0& BX& BY& 0& 0\\ 0& 0& 1& 0& 0\\ AX\ast BY& BX\ast AY& AY\ast BY& AX\ast BX& 0\\ AX\ast BY& BX\ast AY& AY\ast BY& AX\ast BX& 1\end{array}\right]$$
初始矩阵为
$$ma=\left[\begin{array}{c}A0(a_{i-1})\\ B0(b_{i-1})\\ 1\\ A0\ast B0((a_{i-1})\ast (b_{i-1}))\\ A0\ast B0(AoD)\end{array}\right]$$

#define _CRT_SECURE_NO_WARNINGS
#include<iostream>
#include<cstdio>
#include<string>
#include<cstring>
#include<algorithm>
#include<queue>
#include<stack>
#include<cmath>
#include<vector>
#include<fstream>
#include<set>
#include<map>
#include<sstream>
#include<iomanip>
#define ll long long
using namespace std;
//矩阵快速幂
const int mod = 1e9 + 7;
struct matrix {
    ll mat[5][5];
    matrix() { memset(mat, 0, sizeof(mat)); }
    matrix operator * (const matrix& b)const {
        matrix ans;
        for (int i = 0; i < 5; i++) {
            for (int j = 0; j < 5; j++) {
                ans.mat[i][j] = 0;
                for (int k = 0; k < 5; k++) {
                    ans.mat[i][j] = (ans.mat[i][j] + mat[i][k] * b.mat[k][j] % mod + mod) % mod;
                }
            }
        }
        return ans;
    }
};
matrix q_pow(matrix a, ll b) {
    matrix ans;
    memset(ans.mat, 0, sizeof(ans.mat));
    for (int i = 0; i < 5; i++) ans.mat[i][i] = 1;
    while (b) {
        if (b & 1) ans = ans * a;
        b >>= 1;
        a = a * a;
    }
    return ans;
}
ll n, a0, ax, ay, b0, bx, by;
int main() {
    while (~scanf("%lld%lld%lld%lld%lld%lld%lld", &n, &a0, &ax, &ay, &b0, &bx, &by)) {
        a0 %= mod, ax %= mod, ay %= mod, b0 %= mod, bx %= mod, by %= mod;
        if (n == 0) {
            puts("0");
        }
        else if (n == 1) {
            printf("%lld\n", a0 * b0 % mod);
        }
        else {
            matrix ma;
            ma.mat[0][0] = ax % mod, ma.mat[0][2] = ay % mod;
            ma.mat[1][1] = bx % mod, ma.mat[1][2] = by % mod;
            ma.mat[2][2] = 1;
            ma.mat[3][0] = ax * by % mod, ma.mat[3][1] = bx * ay % mod, ma.mat[3][2] = ay * by % mod, ma.mat[3][3] = ax * bx % mod;
            ma.mat[4][0] = ax * by % mod, ma.mat[4][1] = bx * ay % mod, ma.mat[4][2] = ay * by % mod, ma.mat[4][3] = ax * bx % mod, ma.mat[4][4] = 1;
            ma = q_pow(ma, n - 1);
            ll ans = ma.mat[4][0] * a0 % mod + ma.mat[4][1] * b0 % mod + ma.mat[4][2] + ma.mat[4][3] * (a0 * b0 % mod) % mod + a0 * b0 % mod;
            printf("%lld\n", ans % mod);
        }
    }
	return 0;
}