HDU-4686 Arc of Dream(推公式+矩阵快速幂)
function:
where
a 0 = A0
a i = a i-1*AX+AY
b 0 = B0
b i = b i-1*BX+BY
What is the value of AoD(N) modulo 1,000,000,007?
Input
N
A0 AX AY
B0 BX BY
N <= 1e18, other integers <= 2e9
思路:
将an*bn = (an-1 * AX + AY) * (bn-1 * BX + BY) 展开得
(an-1 * bn-1) * AX * BX +
an-1 * AX * BY +
bn-1 * BX * AY +
AY * BY
所以便得到了an*bn的递推公式,然后再新加一行来求前n项和即可。如图5*5矩阵,然后利用矩阵快速幂求解即可。
代码:
#include <algorithm>
#include <string.h>
#include <cstdio>
#define LL long long
using namespace std;
const LL mod = 1e9+7;
struct node
{
LL m[5][5];
node(){memset(m, 0, sizeof m);}
} base, ans;
LL N, A0, AX, AY, B0, BX, BY;
LL AB, A, B, S;
inline void init()
{
base.m[0][0] = AX*BX%mod, base.m[0][1] = AX*BY%mod;
base.m[0][2] = AY*BX%mod, base.m[0][3] = AY*BY%mod;
base.m[0][4] = 0, base.m[1][0] = 0;
base.m[1][1] = AX, base.m[1][2] = 0;
base.m[1][3] = AY, base.m[1][4] = 0;
base.m[2][0] = 0, base.m[2][1] = 0;
base.m[2][2] = BX, base.m[2][3] = BY;
base.m[2][4] = 0, base.m[3][0] = 0;
base.m[3][1] = 0, base.m[3][2] = 0;
base.m[3][3] = 1, base.m[3][4] = 0;
base.m[4][0] = AX*BX%mod, base.m[4][1] = AX*BY%mod;
base.m[4][2] = AY*BX%mod, base.m[4][3] = AY*BY%mod;
base.m[4][4] = 1;
ans.m[0][0] = 1, ans.m[0][1] = 0;
ans.m[0][2] = 0, ans.m[0][3] = 0;
ans.m[0][4] = 0, ans.m[1][0] = 0;
ans.m[1][1] = 1, ans.m[1][2] = 0;
ans.m[1][3] = 0, ans.m[1][4] = 0;
ans.m[2][0] = 0, ans.m[2][1] = 0;
ans.m[2][2] = 1, ans.m[2][3] = 0;
ans.m[2][4] = 0, ans.m[3][0] = 0;
ans.m[3][1] = 0, ans.m[3][2] = 0;
ans.m[3][3] = 1, ans.m[3][4] = 0;
ans.m[4][0] = 0, ans.m[4][1] = 0;
ans.m[4][2] = 0, ans.m[4][3] = 0;
ans.m[4][4] = 1;
}
node mul(node a, node b)
{
node ans1;
for(int i = 0; i < 5; ++i)
for(int j = 0; j < 5; ++j)
for(int k = 0; k < 5; ++k)
{
ans1.m[i][j] += a.m[i][k]*b.m[k][j]%mod;
ans1.m[i][j] %= mod;
}
return ans1;
}
node qpow(LL b)
{
init();
while(b)
{
if(b&1) ans = mul(ans, base);
base = mul(base, base);
b >>= 1;
}
return ans;
}
int main()
{
while(~scanf("%lld", &N))
{
scanf("%lld %lld %lld", &A0, &AX, &AY);
scanf("%lld %lld %lld", &B0, &BX, &BY);
if(!N){puts("0"); continue;}
A0 %= mod, AX %= mod, AY %= mod;
B0 %= mod, BX %= mod, BY %= mod;
AB = A0*B0%mod;
A = A0, B = B0, S = AB;
node res = qpow(N-1);
LL final = res.m[4][0]*AB%mod;
final += res.m[4][1]*A%mod; final %= mod;
final += res.m[4][2]*B%mod; final %= mod;
final += res.m[4][3]; final %= mod;
final += res.m[4][4]*S%mod; final %= mod;
printf("%lld\n", final);
}
return 0;
}
继续加油~