HDU 4686 Arc of Dream(快速幂矩阵)

#include <cstring>
#include <cstdio>
#include <string>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
using namespace std;
#define MOD 1000000007
#define LL long long
LL mat[5][5],p[5][5];
LL a1,b1,a0,b0;
int qmod(LL n)
{
    LL c[5][5];
    int i,j,k;
    LL ans = 0;
    while(n)
    {
        if(n&1)
        {
            memset(c,0,sizeof(c));
            for(i = 0; i < 5; i ++)
            {
                for(j = 0; j < 5; j ++)
                {
                    for(k = 0; k < 5; k ++)
                    {
                        c[i][j] += mat[i][k]*p[k][j];
                        c[i][j] %= MOD;
                    }
                }
            }
            memcpy(mat,c,sizeof(mat));
        }
        memset(c,0,sizeof(c));
        for(i = 0; i < 5; i ++)
        {
            for(j = 0; j < 5; j ++)
            {
                for(k = 0; k < 5; k ++)
                {
                    c[i][j] += p[i][k]*p[k][j];
                    c[i][j] %= MOD;
                }
            }
        }
        memcpy(p,c,sizeof(p));
        n >>= 1;
    }
    ans = ((mat[0][0]*a0)%MOD*b0)%MOD;
    ans = (ans + (mat[0][1]*a1)%MOD*b1)%MOD;
    ans = (ans +  (a1*mat[0][2])%MOD)%MOD;
    ans = (ans + (b1*mat[0][3])%MOD)%MOD;
    ans = (ans + mat[0][4])%MOD;
    return ans;
}
int main()
{
    LL n,ax,ay,bx,by;
    int i,j;
    while(cin>>n)
    {
        cin>>a0>>ax>>ay;
        cin>>b0>>bx>>by;
        memset(p,0,sizeof(p));
        a0 %= MOD;
        ax %= MOD;
        ay %= MOD;
        b0 %= MOD;
        bx %= MOD;
        by %= MOD;
        p[0][0] = p[0][1] = 1;
        p[1][1] = (ax*bx)%MOD;
        p[1][2] = (ax*by)%MOD;
        p[1][3] = (ay*bx)%MOD;
        p[1][4] = (ay*by)%MOD;
        p[2][2] = ax;
        p[2][4] = ay;
        p[3][3] = bx;
        p[3][4] = by;
        p[4][4] = 1;
        for(i = 0; i < 5; i ++)
        {
            for(j = 0; j < 5; j ++)
                mat[i][j] = i == j ? 1:0;
        }
        a1 = (a0*ax + ay)%MOD;
        b1 = (b0*bx + by)%MOD;
        if(n == 0)
            printf("0\n");
        else
            printf("%d\n",qmod(n-1));
    }
    return 0;
}