由式子的性质发现都是线性的，考虑构造矩阵，先有式子，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]即为所求。
贴代码：
1 #include<cstdio> 2 #include<cstring> 3 typedef long long int ll; 4 const int p = 1000000007; 5 ll ax,ay,bx,by,a0,b0; 6 struct matrix 7 { 8 ll m[7][7]; 9 } A; 10 inline void init() 11 { 12 memset(A.m,0,sizeof(A.m)); 13 A.m[1][1] =ax; 14 A.m[1][5] = ay; 15 A.m[2][2] = bx; 16 A.m[2][5] = by; 17 A.m[3][1] = ax*by%p; 18 A.m[3][2] = ay*bx%p; 19 A.m[3][3] = ax*bx%p; 20 A.m[3][5] = ay*by%p; 21 A.m[4][3] = A.m[4][4] = A.m[5][5] = 1; 22 } 23 inline matrix mul(ll a[7][7],ll b[7][7]) 24 { 25 matrix ans; 26 memset(ans.m,0,sizeof(ans.m)); 27 for(int i=1; i<=5; ++i) 28 for(int j=1; j<=5; ++j) 29 for(int k=1; k<=5; ++k) 30 ans.m[i][j] = (ans.m[i][j] + a[i][k]*b[k][j]%p)%p; 31 return ans; 32 } 33 inline matrix qPow(ll x) 34 { 35 matrix ans; 36 memset(ans.m,0,sizeof(ans.m)); 37 for(int i=1; i<=5; ++i) 38 ans.m[i][i] =1; 39 init(); 40 while(x) 41 { 42 if(x&1) ans = mul(ans.m,A.m); 43 A = mul(A.m,A.m); 44 x >>= 1; 45 } 46 return ans; 47 } 48 int main() 49 { 50 // freopen("in.txt","r",stdin); 51 ll n; 52 while(~scanf("%I64d",&n)) 53 { 54 scanf("%I64d%I64d%I64d%I64d%I64d%I64d",&a0,&ax,&ay,&b0,&bx,&by); 55 matrix ans = qPow(n); 56 ll res=0; 57 res = (res + ans.m[4][1]*a0%p)%p; 58 res = (res + ans.m[4][2]*b0%p)%p; 59 res = (res + ans.m[4][3]*((a0*b0)%p)%p)%p; 60 res = (res + ans.m[4][5])%p; 61 printf("%I64d ",res); 62 } 63 return 0; 64 }