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 longusing 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;}

继续加油~