Hdu 4686 Arc of Dream 矩阵快速幂

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

题目大意：已知a[0],ax,ay,  其中a[i] = ax * a[i-1] + ay;

         已知b[0],bx,by   其中b[i] = bx * b[i-1] + by;

         求 a[i]*b[i] (0 =< i <= (n-1)) 的和；

代码如下：

/* 由a[n-1]*b[n-1] = a[n-2]*b[n-2]*ax*bx + a[n-2]*ax*by + b[n-2]*ay*bx + ay * by; 
*  设前n项和为sum[n],则：
*     sum[n]=sum[n-1]+a[n-1]*b[n-1];        
*           =sum[n-1]+a[n-2]*b[n-2]*ax*bx + a[n-2]*ax*by + b[n-2]*ay*bx + ay * by;
                                                                                   {   1     0     0   0   0 } 
*                                                                                  { ax*bx  ax*bx  0   0   0 } 
*{sum[n],a[n-1]*b[n-1],a[n-1],b[n-1],1} = {sum[n-1],a[n-2]*b[n-2],a[n-1],b[n-2],1}*{ ax*by  ax*by  ax  0   0 } 
*                                                                                  { ay*bx  ay*bx  0   bx  0 } 
*                                                                                  { ay*by  ay*by  ay  by  1 } 
*/
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LL long long
#define mod 1000000007
LL n;
struct mat
{
    LL ma[5][5];
};
mat cal(mat a,mat b)
{
    mat c;
    memset(c.ma,0,sizeof(c.ma));
    for(int i=0;i<5;i++)
    for(int k=0;k<5;k++)
    for(int j=0;j<5;j++)
    c.ma[i][j]=(c.ma[i][j]+a.ma[i][k]*b.ma[k][j])%mod;
    return c;
}
mat pow(mat a,LL k)
{
    mat e;
    memset(e.ma,0,sizeof(e.ma));
    for(int i=0;i<5;i++)
        e.ma[i][i]=1;
    while(k>0)
    {
        if(k&1) e=cal(e,a);
        a=cal(a,a);
        k>>=1;
    }
    return e;
}

int main()
{
    LL a0,ax,ay,b0,bx,by;
    while(~scanf("%I64d",&n))
    {
        scanf("%I64d%I64d%I64d%I64d%I64d%I64d",&a0,&ax,&ay,&b0,&bx,&by);
        if(n==0) {printf("0\n");continue;}
        if(n==1) {printf("%I64d\n",a0*b0%mod);continue;}
        mat x;
        memset(x.ma,0,sizeof(x.ma));
        x.ma[0][0]=1;x.ma[0][1]=ax*bx%mod;x.ma[0][2]=ax*by%mod;x.ma[0][3]=bx*ay%mod;x.ma[0][4]=ay*by%mod;
        x.ma[1][0]=0;x.ma[1][1]=ax*bx%mod;x.ma[1][2]=ax*by%mod;x.ma[1][3]=bx*ay%mod;x.ma[1][4]=ay*by%mod;
        x.ma[2][0]=0;x.ma[2][1]=0;x.ma[2][2]=ax%mod;x.ma[2][3]=0;x.ma[2][4]=ay%mod;
        x.ma[3][0]=0;x.ma[3][1]=0;x.ma[3][2]=0;x.ma[3][3]=bx%mod;x.ma[3][4]=by%mod;
        x.ma[4][0]=0;x.ma[4][1]=0;x.ma[4][2]=0;x.ma[4][3]=0;x.ma[4][4]=1;
        x=pow(x,n-1);
        LL k=a0*b0%mod;
        LL ans=(x.ma[0][0]*k+k*x.ma[0][1]+a0*x.ma[0][2]+b0*x.ma[0][3]+x.ma[0][4])%mod;
        printf("%I64d\n",ans);
    }
    return 0;
}