S - Arc of Dream 矩阵快速幂

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<cmath>
#include<map>
#include<stack>
#include<set>
#include<string>
using namespace std;
typedef unsigned long long LL;
#define MAXN 30
#define L 100006
#define MOD 1000000007
#define INF 1000000009
const double eps = 1e-9;
/*
这个题目跟之前做的有一定区别，之前做的大多是表示一个序列的转移矩阵（用一列 列向量表示一个序列的几项）
而这个题给出了 an 的转移关系 bn的转移关系 要求 an*bn的前n项和
应当扩充列的范围（矩阵表达的含义增加）
【An Bn An*Bn Sum(n)】  转移关系是很好列的，剩下的就是矩阵快速幂
*/
LL n, a0, ax, ay, b0, bx, by;
struct Mat
{
    LL a[5][5];
    Mat()
    {
        memset(a, 0, sizeof(a));
    }
    Mat operator*(const Mat& rhs)
    {
        Mat ret;
        for (int i = 0; i < 5; i++)
        {
            for (int j = 0; j < 5; j++)
            {
                if(a[i][j])
                {
                    for (int k = 0; k < 5; k++)
                    {
                        ret.a[i][k] = (ret.a[i][k] + a[i][j] * rhs.a[j][k]) % MOD;
                    }
                }
            }
        }
        return ret;
    }
};
Mat fpow(Mat m, LL b)
{
    Mat ans;
    for (int i = 0; i < 5; i++)
        ans.a[i][i] = 1;
    while (b != 0)
    {
        if (b & 1)
            ans = m * ans;
        m = m * m;
        b = b / 2;
    }
    return ans;
}
int main()
{
    while (scanf("%lld", &n) != EOF)
    {
        scanf("%lld%lld%lld%lld%lld%lld", &a0, &ax, &ay, &b0, &bx, &by);
        if (n == 0)
        {
            printf("0\n");
            continue;
        }
        Mat M;
        M.a[0][0] = ax%MOD, M.a[0][2] = ax%MOD*by%MOD, M.a[0][3] = ax%MOD*by%MOD;
        M.a[1][1] = bx%MOD, M.a[1][2] = M.a[1][3] = bx%MOD*ay%MOD;
        M.a[2][2] = M.a[2][3] = ax%MOD*bx%MOD;
        M.a[3][3] = 1;
        M.a[4][0] = ay%MOD, M.a[4][1] = by%MOD, M.a[4][2] = M.a[4][3] = ay%MOD*by%MOD, M.a[4][4] = 1;
        if (n == 1)
            printf("%lld\n", a0*b0%MOD);
        else
        {
            M = fpow(M, n - 1);
            LL ans = 0;
            ans = (ans + a0*M.a[0][3] % MOD) % MOD;
            ans = (ans + b0*M.a[1][3] % MOD) % MOD;
            ans = (ans + +a0%MOD*b0%MOD*M.a[2][3] % MOD) % MOD;
            ans = (ans + a0%MOD*b0%MOD*M.a[3][3] % MOD) % MOD;
            ans = (ans +M.a[4][3] % MOD) % MOD;
            printf("%lld\n", ans);
        }
    }
}