HDU4686 Arc of Dream 矩阵快速幂

#include<bits/stdc++.h>
#define N 5
#define mes(x) memset(x, 0, sizeof(x));
#define ll long long
const ll mod = 1e9+7;
const int MAX = 0x7ffffff;
using namespace std;
struct mat {
    ll a[N][N];  
    mat() {  
        memset(a, 0, sizeof(a));  
    } 
    mat operator * (mat b) {  
        mat c;  
        for (int i = 0; i < N; i++)  
            for (int j = 0; j < N; j++)  
                for (int k = 0; k < N; k++)  
                    c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j]) % mod;  
        return c;  
    }
};
mat f(mat b, ll m) {
    mat c;
    for (int i = 0; i < N; i++)  
        c.a[i][i] = 1;
    while (m) {
        if (m & 1)
            c = c * b;  
        b = b * b;  
        m >>= 1;  
    }
    return c;  
}
int main()
{
    ll n, A0,Ax, Ay, Bx,By,B0;
    mat A, s;
    while(~scanf("%lld", &n)){
        scanf("%lld%lld%lld%lld%lld%lld", &A0, &Ax, &Ay, &B0, &Bx, &By);
        if(n == 0){
            printf("0\n");
            continue;
        }
        mes(A.a);
        mes(s.a);
        s.a[0][0] = s.a[0][1] = (A0%mod*B0%mod)%mod;
        s.a[0][2] = B0%mod;
        s.a[0][3] = A0%mod;
        s.a[0][4] = 1;
        A.a[0][0] = 1; 
        A.a[1][1] = A.a[1][0] = (Ax%mod*Bx%mod)%mod;A.a[2][2] = Bx%mod;
        A.a[2][1] = A.a[2][0] = (Ay%mod*Bx%mod)%mod;A.a[3][3] = Ax%mod;
        A.a[3][1] = A.a[3][0] = (Ax%mod*By%mod)%mod;
        A.a[4][0] = A.a[4][1] = (Ay%mod*By%mod)%mod;
        A.a[4][2] = By;
        A.a[4][3] = Ay;
        A.a[4][4] = 1;
        A = f(A,n-1);
        s = s*A;
        printf("%lld\n", (mod+s.a[0][0])%mod);
    } 
}