HDU4686——Arc of Dream矩阵快速幂

题目链接：

　　http://acm.hdu.edu.cn/showproblem.php?pid=4686

题目大意：

　　已知a0=A0,　　a_i=Ax*a_i-1+Ay;

　　　　b0=B0,　　b_i=Bx*b_i-1+By;

　　求∑a_i*b_i（i=0-->n-1）。

　　n不超过10¹⁸，A0，Ax，Ay，B0，Bx，BY不超过2*10⁹。

题目分析：

　　因为n很大，不可能用递推来做，这个时候就想到了矩阵的方法。构造了好几个满足要求的，但都是仅仅满足ai或者bi的，最后才发现，把ai*bi按递推式展开，

a_i*b_i=Ax*By*a_i-1*b_i-1+Ax*By*a_i-1+Ay*Bx*b_i-1+By*Ay。将常数组合在一起构成一个矩阵，将变量组合在一起构成另一个矩阵，然后将ai*bi构造成矩阵递推式：

矩阵1：

　　1　　a_i　　b_i　　a_i*b_i　　si（求和）

　　0　　0　　0　　    0　　   0

　　0　　0　　0　　    0　　   0

　　0　　0　　0　　    0　　   0

　　0　　0　　0　　    0　　   0

矩阵2：

　　1　　Ay　　By　　Ay*By　　Ay*By

　　0　　Ax　　0　　  Ax*By　　Ax*By

　　0　　0　　  Bx　　Ay*Bx　　Ay*Bx

　　0　　0　　  0　　  Ax*By　　Ax*By

　　0　　0　　  0　　  0　　　　 1

矩阵3

　　1　　a_i+1　　b_i+1　　a_i+1*b_i+1　　si（求和）

　　0　　0　　0　　    0　　   0

　　0　　0　　0　　    0　　   0

　　0　　0　　0　　    0　　   0

　　0　　0　　0　　    0　　   0

显然　　矩阵1*矩阵2=矩阵3。根据递推关系呢，矩阵1（i=0）*（矩阵2）^n-1就能得到s（n-1）了。因而，用矩阵快速幂就能很快把问题解决了。

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
const long long mod=1000000007;
typedef struct
{
    long long m[5][5];
}mat;
mat X,Y;
mat multi(mat x,mat y)
{
    mat temp;
    for(int i=0;i<5;i++)
        for(int j=0;j<5;j++)
        {
            temp.m[i][j]=0;
            for(int k=0;k<5;k++)
                temp.m[i][j]+=x.m[i][k]*y.m[k][j]%mod;
            temp.m[i][j]%=mod;
        }
    return temp;
}

mat pow(long long k)//矩阵快速幂
{
    mat ans=X,p=Y;
    while(k)
    {
        if(k&1)
            ans=multi(ans,p);
        p=multi(p,p);
        k/=2;
    }
    return ans;
}

int main()
{
    long long n,a0,ax,ay,b0,bx,by;
    while(cin>>n>>a0>>ax>>ay>>b0>>bx>>by)
    {
        if(!n)//这边需要注意特判一下
        {
            printf("0
");
            continue;
        }
        memset(X.m,0,sizeof(X.m));
        memset(Y.m,0,sizeof(Y.m));
        X.m[0][0]=1;X.m[0][1]=a0;X.m[0][2]=b0;X.m[0][3]=a0*b0%mod;X.m[0][4]=a0*b0%mod;
        Y.m[0][0]=1;Y.m[0][1]=ay;Y.m[0][2]=by;Y.m[0][3]=ay*by%mod;Y.m[0][4]=ay*by%mod;
        Y.m[1][1]=ax;Y.m[1][3]=Y.m[1][4]=ax*by%mod;
        Y.m[2][2]=bx;Y.m[2][3]=Y.m[2][4]=ay*bx%mod;
        Y.m[3][3]=Y.m[3][4]=ax*bx%mod;
        Y.m[4][4]=1;
        mat ans=pow(n-1);
        long long s=ans.m[0][4]%mod;
        cout<<s<<endl;
    }
    return 0;
}