HDUOJ 4686 Arc of Dream

Problem Description

An Arc of Dream is a curve defined by following function:

where
a0=A0a_0 = A_0a0=A0
ai=ai−1∗AX+AYa_i = a_{i-1}*A_X+A_Yai=ai−1∗AX+AY
b0=B0b_0 = B_0b0=B0
bi=bi−1∗BX+BYb_i = b_{i-1}*B_X+B_Ybi=bi−1∗BX+BY
What is the value of AoD(N) modulo 1,000,000,007?

Input

There are multiple test cases. Process to the End of File.
Each test case contains 7 nonnegative integers as follows:
N
A0 AX AY
B0 BX BY
N is no more than 1018, and all the other integers are no more than 2×109.

Output

For each test case, output AoD(N) modulo 1,000,000,007.

Sample Input

Sample Output

4
134
1902

题目其实不难，简单推理一下即可：
Ai−1=a0b0+a1b1+⋯+ai−2bi−2A_{i-1}=a_0b_0+a_1b_1+\cdots+a_{i-2}b_{i-2}Ai−1=a0b0+a1b1+⋯+ai−2bi−2
Ai=a0b0+a1b1+⋯+ai−2bi−2+ai−1bi−1A_i=a_0b_0+a_1b_1+\cdots+a_{i-2}b_{i-2}+a_{i-1}b_{i-1}Ai=a0b0+a1b1+⋯+ai−2bi−2+ai−1bi−1
＝Ai−1+ai−1bi−1\ \ \ \ \ ＝A_{i-1}+a_{i-1}b_{i-1}     ＝Ai−1+ai−1bi−1
＝Ai−1+(AXai−2+AY)(BXbi−2+BY)\ \ \ \ \ ＝A_{i-1}+(A_Xa_{i-2}+A_Y)(B_Xb_{i-2}+B_Y)     ＝Ai−1+(AXai−2+AY)(BXbi−2+BY)
＝Ai−1+AXBXai−2bi−2+AXBYai−2+AYBXbi−2+AYBY\ \ \ \ \ ＝A_{i-1}+A_XB_Xa_{i-2}b_{i-2}+A_XB_Ya_{i-2}+A_YB_Xb_{i-2}+A_YB_Y     ＝Ai−1+AXBXai−2bi−2+AXBYai−2+AYBXbi−2+AYBY
那么递推矩阵就出来了：
[Aiai−1bi−1ai−1bi−11]=[1AXBXAXBYAYBXAYBY0AXBXAXBYAYBXAYBY00AX0AY000BXBY00001][Ai−1ai−2bi−2ai−2bi−21]\left[ \begin{matrix} A_i\\ a_{i-1}b_{i-1}\\ a_{i-1}\\ b_{i-1}\\ 1\\ \end{matrix} \right]= \left[ \begin{matrix} 1&A_XB_X&A_XB_Y&A_YB_X&A_YB_Y\\ 0&A_XB_X&A_XB_Y&A_YB_X&A_YB_Y\\ 0&0&A_X&0&A_Y\\ 0&0&0&B_X&B_Y\\ 0&0&0&0&1\\ \end{matrix} \right]\left[ \begin{matrix} A_{i-1}\\ a_{i-2}b_{i-2}\\ a_{i-2}\\ b_{i-2}\\ 1\\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎡Aiai−1bi−1ai−1bi−11⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡10000AXBXAXBX000AXBYAXBYAX00AYBXAYBX0BX0AYBYAYBYAYBY1⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡Ai−1ai−2bi−2ai−2bi−21⎦⎥⎥⎥⎥⎤还有两个坑点：
1.n=0n=0n=0 时特判 A0=0A_0=0A0=0
2.中间过程多取余防爆，我个人是加了快速乘，不会的就多加几个取模就行了
AC代码如下：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
ll n,k;
struct mat
{ll m[5][5];
}a,ans;void init(mat &a){memset(a.m,0,sizeof(a.m));for(int i=0;i<5;i++) a.m[i][i]=1;
}mat mul(mat a,mat b)
{mat c;memset(c.m,0,sizeof(c.m));for(int i=0;i<5;i++){for(int k=0;k<5;k++){if(a.m[i][k]==0) continue;for(int j=0;j<5;j++){c.m[i][j]=(c.m[i][j]+a.m[i][k]%mod*b.m[k][j]%mod)%mod;if(c.m[i][j]<0) c.m[i][j]+=mod;}}}return c;
}mat mat_pow(mat a,ll k)
{mat ans;init(ans);while(k>0){if(k&1) ans=mul(ans,a);a=mul(a,a);k>>=1;}return ans;
}ll f(ll a,ll b)
{ll k=0;while(b){if(b&1) k=(k+a)%mod;a=(a+a)%mod;b>>=1;}return k;
}int main(){while(cin>>n){ll a0,ax,ay,b0,bx,by;cin>>a0>>ax>>ay>>b0>>bx>>by;if(n==0) {puts("0");continue;}memset(a.m,0,sizeof(a.m));memset(ans.m,0,sizeof(ans.m));a.m[0][0]=1,a.m[0][1]=f(ax,bx),a.m[0][2]=f(ax,by),a.m[0][3]=f(ay,bx),a.m[0][4]=f(ay,by);a.m[1][1]=f(ax,bx),a.m[1][2]=f(ax,by),a.m[1][3]=f(ay,bx),a.m[1][4]=f(ay,by);a.m[2][2]=ax%mod,a.m[2][4]=ay%mod;a.m[3][3]=bx%mod,a.m[3][4]=by%mod;a.m[4][4]=1;ans.m[0][0]=f(a0,b0),ans.m[1][0]=f(a0,b0),ans.m[2][0]=a0%mod,ans.m[3][0]=b0%mod,ans.m[4][0]=1;ans=mul(mat_pow(a,n-1),ans);cout<<ans.m[0][0]<<endl;}return 0;
}