HDU 4686 Arc of Dream（矩阵快速幂）

Description
已知，其中
a0 = A0
ai = ai-1*AX+AY
b0 = B0
bi = bi-1*BX+BY
求出AoD(n)（mod 1000000007）的值
Input
多组输入，每组用例第一行为一整数N，第二行为三个整数A0,AX,AY,第三行为三个整数B0,BX,BY,以文件尾结束输入
Output
对于每组用例，输出AoD(N)(mod 1000000007)
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
Solution
构造矩阵，其中si=AoD(i+1)，那么由
可以得到下列矩阵转移方程

剩下的套矩阵快速幂模板即可，注意两点，其一是n=0需特判，其二是构造原始矩阵和转移矩阵时就需要模
Code

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define maxn 6
typedef long long ll;
const ll mod=1000000007ll;
struct Mat
{
    ll mat[maxn][maxn];//矩阵 
    ll row,col;//矩阵行列数 
};
Mat mod_mul(Mat a,Mat b,ll p)//矩阵乘法 
{
    Mat ans;
    ans.row=a.row;
    ans.col=b.col;
    memset(ans.mat,0,sizeof(ans.mat));
    for(ll i=0;i<ans.row;i++)       
        for(ll k=0;k<a.col;k++)
            if(a.mat[i][k])
                for(ll j=0;j<ans.col;j++)
                {
                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                    ans.mat[i][j]%=p;
                }
    return ans;
}
Mat mod_pow(Mat a,ll k,ll p)//矩阵快速幂 
{
    Mat ans;
    ans.row=a.row;
    ans.col=a.col;
    for(ll i=0;i<a.row;i++)
        for(ll j=0;j<a.col;j++)
            ans.mat[i][j]=(i==j);
    while(k)
    {
        if(k&1)ans=mod_mul(ans,a,p);
        a=mod_mul(a,a,p);
        k>>=1;
    }
    return ans;
}
int main()
{
    ll N,A0,AX,AY,B0,BX,BY;
    while(scanf("%lld%lld%lld%lld%lld%lld%lld",&N,&A0,&AX,&AY,&B0,&BX,&BY)!=EOF)
    {
        if(!N)//N=0特判 
        {
            printf("0\n");
            continue;
        }
        Mat A,B,ans;
        //初始化 
        memset(A.mat,0,sizeof(A.mat));
        memset(B.mat,0,sizeof(B.mat));
        //构造原始矩阵 
        A.row=1,A.col=5;
        A.mat[0][0]=1,A.mat[0][1]=A0%mod,A.mat[0][2]=B0%mod,A.mat[0][3]=A.mat[0][4]=A0*B0%mod;
        //构造转移矩阵 
        B.row=B.col=5;
        B.mat[0][0]=B.mat[4][4]=1;
        B.mat[0][1]=AY%mod;B.mat[1][1]=AX%mod;
        B.mat[0][2]=BY%mod;B.mat[2][2]=BX%mod;
        B.mat[0][3]=B.mat[0][4]=AY*BY%mod;
        B.mat[1][3]=B.mat[1][4]=AX*BY%mod;
        B.mat[2][3]=B.mat[2][4]=AY*BX%mod;
        B.mat[3][3]=B.mat[3][4]=AX*BX%mod;
        ans=mod_pow(B,N-1,mod);//矩阵快速幂 
        ans=mod_mul(A,ans,mod);//A*B^(N-1)即为答案 
        printf("%lld\n",ans.mat[0][4]);
    } 
    return 0;
}