S - Arc of Dream 矩陣快速冪

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

where 
a  ₀ = A0 
a  _i = a  _i-1*AX+AY 
b  ₀ = B0 
b  _i = b  _i-1*BX+BY 
What is the value of AoD(N) modulo 1,000,000,007?

InputThere are multiple test cases. Process to the End of File. 
Each test case contains 7 nonnegative integers as follows: 
N 
A0 AX AY 
B0 BX BY 
N is no more than 10 ¹⁸, and all the other integers are no more than 2×10 ⁹.OutputFor each test case, output AoD(N) modulo 1,000,000,007.Sample Input

1
1 2 3
4 5 6
2
1 2 3
4 5 6
3
1 2 3
4 5 6

Sample Output

4
134
1902

                                                            |   AX   0   AXBY   AXBY  0  |

                                                                |   0   BX  AYBX    AYBX  0  |

{a[i-1]   b[i-1]   a[i-1]*b[i-1]  AoD[i-1]  1}* |   0   0   AXBX    AXBX   0  |  = {a[i]   b[i]   a[i]*b[i]  AoD[i]  1}

                                                               |    0   0     0          1     0    |

                                                               |  AY  BY   AYBY   AYBY   1   |

#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);
        }
    }
}