解题思路：
aibi=(ai-1*AX+AY)(bi-1*BX+BY)=AX*BX*ai-1bi-1+AX*BY*ai-1+BX*AY*bi-1+AY*BY
n-1 n-1 n-1 n-1
∑ aibi = AX*BX*∑ ai-1bi-1 + AX*BY*∑ ai-1 + BX*AY*∑ bi-1 + (n-1)*AY*BY +A0*B0
i=0 i=0 i=0 i=0
n-1 n-1
∑ ai = AX*∑ ai-1 + (n-1)*AY + A0
i=0 i=0
n-1 n-1
∑ bi = BX*∑ bi-1 + (n-1)*BY + B0
i=0 i=0
所以有：
注意1和0的时候要特判。
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #define mod 1000000007 5 #define ll long long 6 using namespace std; 7 struct rec{ 8 ll v[5][5]; 9 rec(){memset(v,0,sizeof v);} 10 }; 11 void debug(rec a){ 12 cout<<a.v[0][0]<<" "<<a.v[0][1]<<endl; 13 cout<<a.v[0][2]<<" "<<a.v[0][3]<<endl; 14 } 15 rec mul(rec a,rec b){ 16 rec ans; 17 ll sum; 18 //debug(ans); 19 //debug(a),debug(b); 20 for(int i=0;i<5;i++){ 21 for(int j=0;j<5;j++){ 22 for(int k=0;k<5;k++){ 23 //cout<<a.v[i][k]<<" "<<b.v[k][j]<<endl; 24 sum=a.v[i][k]%mod*b.v[k][j]%mod; 25 ans.v[i][j]=(ans.v[i][j]+sum%mod)%mod; 26 } 27 //cout<<ans.v[i][j]<<endl; 28 } 29 } 30 return ans; 31 } 32 33 rec pow(rec a,ll n){ 34 if(n==1) return a; 35 rec ans=pow(a,n/2); 36 ans=mul(ans,ans); 37 if(n&1) return mul(ans,a); 38 else return ans; 39 } 40 ll n,a1,b1; 41 ll a0,ax,ay,b0,bx,by; 42 ll f(rec a){ 43 ll x1=a.v[0][0],x2=a.v[0][1],x3=a.v[0][2],x4=a.v[0][3],x5=a.v[0][4]; 44 return ((((x1*a0%mod*b0%mod+x2*a0%mod)%mod+x3*b0%mod)%mod+x4)%mod+x5)%mod; 45 } 46 47 int main(){ 48 while(cin>>n){ 49 cin>>a0>>ax>>ay; 50 cin>>b0>>bx>>by; 51 if(n==0) {cout<<0<<endl;continue;} 52 if(n==1) {cout<<a0*b0%mod<<endl;continue;} 53 rec a; 54 a.v[0][0]=ax%mod*bx%mod,a.v[0][1]=ax%mod*by%mod,a.v[0][2]=bx%mod*ay%mod, 55 a.v[0][3]=ay%mod*by%mod,a.v[0][4]=a0%mod*b0%mod; 56 a.v[1][0]=0,a.v[1][1]=ax%mod,a.v[1][2]=0,a.v[1][3]=ay%mod,a.v[1][4]=a0%mod; 57 a.v[2][0]=0,a.v[2][1]=0,a.v[2][2]=bx%mod,a.v[2][3]=by%mod,a.v[2][4]=b0%mod; 58 a.v[3][0]=0,a.v[3][1]=0,a.v[3][2]=0,a.v[3][3]=1,a.v[3][4]=1; 59 a.v[4][0]=0,a.v[4][1]=0,a.v[4][2]=0,a.v[4][3]=0,a.v[4][4]=1; 60 //debug(a); 61 rec ans=pow(a,n-1); 62 cout<<f(ans)%mod<<endl; 63 } 64 return 0; 65 }