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

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4686

#include<bits/stdc++.h>
using namespace std;

#define debug puts("YES");
#define rep(x,y,z) for(int (x)=(y);(x)<(z);(x)++)
#define ll long long

#define lrt int l,int r,int rt
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
#define root l,r,rt
#define mst(a,b) memset((a),(b),sizeof(a))
const int  maxn =1e3+5;
const int mod=1e9+7;
const int ub=5;
ll powmod(ll x,ll y){ll t; for(t=1;y;y>>=1,x=x*x%mod) if(y&1) t=t*x%mod; return t;}
ll gcd(ll x,ll y){return y?gcd(y,x%y):x;}
/*
把二维递推关系找好就可以了，
ai*bi可以用ai-1*bi-1和ai-1和bi-1和常数来表示，
那么我们列出求解向量为：sn，an*bn，an，bn，1
这个向量列出递推矩阵，然后直接快速幂求解即可。

*/

ll n,a0,ax,ay,b0,bx,by;

struct mat{
    ll a[ub][ub];
    mat()
    {
        memset(a,0,sizeof(a));
        for(int i=0;i<ub;i++) a[i][i]=1;
    }
    mat operator*(const mat& y) const
    {
        mat tp;
        for(int i=0;i<ub;i++) for(int j=0;j<ub;j++)
        {
            tp.a[i][j]=0;
            for(int k=0;k<ub;k++) (tp.a[i][j]+=a[i][k]*y.a[k][j])%=mod;
        }
        return tp;
    }
};

ll quick_pow(mat x,ll n)
{
    mat ret;
    for(;n;n>>=1,x=x*x) if(n&1) ret=x*ret;
    ll ans=0;
    ans=(ans+ret.a[0][1]*a0%mod*b0%mod)%mod;
    ans=(ans+a0*ret.a[0][2]%mod)%mod;
    ans=(ans+b0*ret.a[0][3]%mod)%mod;
    ans=(ans+ret.a[0][4])%mod;
    return ans;
}

int main()
{
    while(scanf("%lld",&n)!=EOF)
    {
        scanf("%lld%lld%lld",&a0,&ax,&ay);
        scanf("%lld%lld%lld",&b0,&bx,&by);
        mat x;memset(x.a,0,sizeof(x.a));
        x.a[0][0]=x.a[0][1]=1;
        x.a[1][1]=ax*bx%mod;x.a[1][2]=ax*by%mod;
        x.a[1][3]=ay*bx%mod;x.a[1][4]=ay*by%mod;
        x.a[2][2]=ax%mod;x.a[2][4]=ay%mod;
        x.a[3][3]=bx%mod;x.a[3][4]=by%mod;
        x.a[4][4]=1;
        ll ans=quick_pow(x,n);
        printf("%lld\n",ans);
    }
    return 0;
}