HDU4686——Arc of Dream矩阵快速幂

题目链接：

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

题目大意：

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

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

　　求∑ai*bi（i=0-->n-1）。

　　n不超过1018，A0，Ax，Ay，B0，Bx，BY不超过2*109。

题目分析：

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

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

矩阵1：

　　1　　ai　　bi　　ai*bi　　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　　ai+1　　bi+1　　ai+1*bi+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\n");
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;
}

HDU4686