HDU 4686 Arc of Dream 矩阵快速幂,线性同余 难度:1

时间:2023-12-21 09:36:14

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

当看到n为小于64位整数的数字时,就应该有个感觉,acm范畴内这应该是道矩阵快速幂

Ai,Bi的递推式题目已经给出,

Ai*Bi=Ax*Bx*(Ai-1*Bi-1)+Ax*By*Ai-1+Bx*Ay*Bi-1+Ay*By

AoD(n)=AoD(n-1)+AiBi

构造向量I{AoD(i-1),Ai*Bi,Ai,Bi,1}

初始向量为I0={0,A0*B0,A0,B0,1}

构造矩阵A{

1,0,0,0,0,

1,Ax*Bx,0,0,0,
0,Ax*By,Ax,0,0,

0,Bx*Ay,0,Bx,0,

0,Ay*By,Ay,By,1,

}

则I0*(A)^n则为包含结果的向量,此时I[0]即为解答

注意:

1.据说有n=0,直接输出0

2.long long的使用

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
const ll mod = 1e9+7;
ll n,ao,bo,ax,bx,ay,by;
void init(ll** ans,ll** base){
fill(ans[0],ans[0]+5,0);
fill(base[0],base[0]+25,0);
ans[0][1]=ao*bo%mod;
ans[0][2]=ao;
ans[0][3]=bo;
ans[0][4]=1;
base[0][0]=1;
base[1][0]=1;base[1][1]=ax*bx%mod;
base[2][1]=ax*by%mod;base[2][2]=ax;
base[3][1]=bx*ay%mod;base[3][3]=bx;
base[4][1]=ay*by%mod;base[4][2]=ay;base[4][3]=by;base[4][4]=1;
}
void print(ll **t ,int m,int r){
for(int i=0;i<m;i++){
for(int j=0;j<r;j++){
printf("%I64d ",t[i][j]);
}
puts("");
}
} void multi(ll** t,ll** b,ll ** tmp,int m,int r,int s){
fill(tmp[0],tmp[0]+m*s,0);
for(int i=0;i<m;i++){
for(int j=0;j<s;j++){
for(int k=0;k<r;k++){
tmp[i][j]+=(t[i][k]*b[k][j])%mod;
tmp[i][j]%=mod;
}
}
}
} void qpow(ll** ans,ll ** base,ll**tmp,ll time,int r,int m){
while(time>0){
if(time&1){
multi(ans,base,tmp,r,m,m);
copy(tmp[0],tmp[0]+r*m,ans[0]);
time--;
}
multi(base,base,tmp,m,m,m);
copy(tmp[0],tmp[0]+m*m,base[0]);
time>>=1;
}
} int main(){
ll ** ans = new ll*[1],**tmp=new ll*[5],** base =new ll*[5];
ans[0]=new ll[5],tmp[0]=new ll[25],base[0]=new ll[25];
for(int i=1;i<5;i++){
tmp[i]=tmp[0]+5*i;
base[i]=base[0]+5*i;
}
while(scanf("%I64d",&n)==1){
scanf("%I64d%I64d%I64d",&ao,&ax,&ay);
ao%=mod;ax%=mod;ay%=mod;
scanf("%I64d%I64d%I64d",&bo,&bx,&by);
bo%=mod;bx%=mod;by%=mod;
if(n==0)puts("0");
else {
init(ans,base);
qpow(ans,base,tmp,n,1,5);
printf("%I64d\n",ans[0][0]);
} }
delete ans;delete base;delete tmp;
return 0;
}