HDU 4686 矩阵快速幂 Arc of Dream

由式子的性质发现都是线性的,考虑构造矩阵,先有式子,a[i] = ax * a[i-1] + ay; b[i] = bx*b[i-1] +by;

a[i]*b[i] = ax*bx*a[i-1]*b[i-1] + ax*by*a[i-1] + bx*ay*b[i-1]+ay*by;

s[i] = s[i-1] + a[i-1]*b[i-1];

由此得到递推式 :设矩阵A=

ax 0 0 0 ay
0 bx 0 0 by
ax*by bx*ay ax*bx 0 ay*by
0 0 1 1 0
0 0 0 0 1

矩阵B[i]=(a[i-1],b[i-1],a[i-1]*b[i-1],s[i-1],1)' (转置),B[i] =(a[i],b[i],a[i]*b[i],s[i],1)' (转置),则有B[i] = A*B[i-1]

令s0 = 0,则有B[0] = (a0,b0,a0*b0,s0,1)',B[n] = A^n*B[0],矩阵乘法是服从结合律的,所以先用矩阵快速幂算出A^n,再算出B[n],那么B[n][4]即为所求。

贴代码:

 #include<cstdio>
#include<cstring>
typedef long long int ll;
const int p = ;
ll ax,ay,bx,by,a0,b0;
struct matrix
{
ll m[][];
} A;
inline void init()
{
memset(A.m,,sizeof(A.m));
A.m[][] =ax;
A.m[][] = ay;
A.m[][] = bx;
A.m[][] = by;
A.m[][] = ax*by%p;
A.m[][] = ay*bx%p;
A.m[][] = ax*bx%p;
A.m[][] = ay*by%p;
A.m[][] = A.m[][] = A.m[][] = ;
}
inline matrix mul(ll a[][],ll b[][])
{
matrix ans;
memset(ans.m,,sizeof(ans.m));
for(int i=; i<=; ++i)
for(int j=; j<=; ++j)
for(int k=; k<=; ++k)
ans.m[i][j] = (ans.m[i][j] + a[i][k]*b[k][j]%p)%p;
return ans;
}
inline matrix qPow(ll x)
{
matrix ans;
memset(ans.m,,sizeof(ans.m));
for(int i=; i<=; ++i)
ans.m[i][i] =;
init();
while(x)
{
if(x&) ans = mul(ans.m,A.m);
A = mul(A.m,A.m);
x >>= ;
}
return ans;
}
int main()
{
// freopen("in.txt","r",stdin);
ll n;
while(~scanf("%I64d",&n))
{
scanf("%I64d%I64d%I64d%I64d%I64d%I64d",&a0,&ax,&ay,&b0,&bx,&by);
matrix ans = qPow(n);
ll res=;
res = (res + ans.m[][]*a0%p)%p;
res = (res + ans.m[][]*b0%p)%p;
res = (res + ans.m[][]*((a0*b0)%p)%p)%p;
res = (res + ans.m[][])%p;
printf("%I64d\n",res);
}
return ;
}
上一篇:Deep_learning


下一篇:Oracle之ORA-00972: identifier is too long