hdu 4686 Arc of Dream 自己推 矩阵快速幂

A.mat[0][0] = 1, A.mat[0][1] = 1, A.mat[0][2] = 0, A.mat[0][3] = 0, A.mat[0][4] = 0;
A.mat[1][0] = 0, A.mat[1][1] = AX*BX%Mod, A.mat[1][2] = AX*BY%Mod, A.mat[1][3] = AY*BX%Mod, A.mat[1][4] = AY*BY%Mod;
A.mat[2][0] = 0, A.mat[2][1] = 0, A.mat[2][2] = AX, A.mat[2][3] = 0, A.mat[2][4] = AY;
A.mat[3][0] = 0, A.mat[3][1] = 0, A.mat[3][2] = 0, A.mat[3][3] = BX, A.mat[3][4] = BY;
A.mat[4][0] = 0, A.mat[4][1] = 0, A.mat[4][2] = 0, A.mat[4][3] = 0, A.mat[4][4] = 1;

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long ll;
const ll Mod = 1e9 + 7;
const int N = 6;
int msize;

struct Mat
{
    ll mat[N][N];
};

Mat operator *(Mat a, Mat b)
{
    Mat c;
    memset(c.mat, 0, sizeof(c.mat));
    for(int k = 0; k < msize; ++k)
        for(int i = 0; i < msize; ++i)
            if(a.mat[i][k])
                for(int j = 0; j < msize; ++j)
                    if(b.mat[k][j])
                        c.mat[i][j] = (a.mat[i][k] * b.mat[k][j] + c.mat[i][j])%Mod;
    return c;
}

Mat operator ^(Mat a, ll k)
{
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    for(int i = 0; i < msize; ++i)
        c.mat[i][i]=1;
    for(; k; k >>= 1)
    {
        if(k&1) c = c*a;
        a = a*a;
    }
    return c;
}

int main()
{
    ll n,A0,B0,AX,AY,BX,BY;
    msize = 5;
    while(~scanf("%I64d",&n))
    {
        scanf("%I64d%I64d%I64d", &A0, &AX, &AY);
        scanf("%I64d%I64d%I64d", &B0, &BX, &BY);
        if(n==0)
        {
            puts("0");
            continue;
        }
        Mat A;
        A.mat[0][0] = 1, A.mat[0][1] = 1,         A.mat[0][2] = 0,         A.mat[0][3] = 0,         A.mat[0][4] = 0;
        A.mat[1][0] = 0, A.mat[1][1] = AX*BX%Mod, A.mat[1][2] = AX*BY%Mod, A.mat[1][3] = AY*BX%Mod, A.mat[1][4] = AY*BY%Mod;
        A.mat[2][0] = 0, A.mat[2][1] = 0,         A.mat[2][2] = AX,        A.mat[2][3] = 0,         A.mat[2][4] = AY;
        A.mat[3][0] = 0, A.mat[3][1] = 0,         A.mat[3][2] = 0,         A.mat[3][3] = BX,        A.mat[3][4] = BY;
        A.mat[4][0] = 0, A.mat[4][1] = 0,         A.mat[4][2] = 0,         A.mat[4][3] = 0,         A.mat[4][4] = 1;
        A = A^(n);
        ll AoD1 = A0*B0%Mod;
        printf("%I64d
",(A.mat[0][1]*AoD1%Mod + A.mat[0][2]*A0%Mod + A.mat[0][3]*B0%Mod + A.mat[0][4])%Mod);
    }
    return 0;
}