由式子的性质发现都是线性的，考虑构造矩阵，先有式子，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]即为所求。
贴代码：
#include<cstdio>
#include<cstring>
typedef long long int ll;
const int p = ;
ll ax,ay,bx,by,a0,b0;
struct matrix
{
ll m[][];
} A;
inline void init()
{
memset(A.m,,sizeof(A.m));
A.m[][] =ax;
A.m[][] = ay;
A.m[][] = bx;
A.m[][] = by;
A.m[][] = ax*by%p;
A.m[][] = ay*bx%p;
A.m[][] = ax*bx%p;
A.m[][] = ay*by%p;
A.m[][] = A.m[][] = A.m[][] = ;
}
inline matrix mul(ll a[][],ll b[][])
{
matrix ans;
memset(ans.m,,sizeof(ans.m));
for(int i=; i<=; ++i)
for(int j=; j<=; ++j)
for(int k=; k<=; ++k)
ans.m[i][j] = (ans.m[i][j] + a[i][k]*b[k][j]%p)%p;
return ans;
}
inline matrix qPow(ll x)
{
matrix ans;
memset(ans.m,,sizeof(ans.m));
for(int i=; i<=; ++i)
ans.m[i][i] =;
init();
while(x)
{
if(x&) ans = mul(ans.m,A.m);
A = mul(A.m,A.m);
x >>= ;
}
return ans;
}
int main()
{
// freopen("in.txt","r",stdin);
ll n;
while(~scanf("%I64d",&n))
{
scanf("%I64d%I64d%I64d%I64d%I64d%I64d",&a0,&ax,&ay,&b0,&bx,&by);
matrix ans = qPow(n);
ll res=;
res = (res + ans.m[][]*a0%p)%p;
res = (res + ans.m[][]*b0%p)%p;
res = (res + ans.m[][]*((a0*b0)%p)%p)%p;
res = (res + ans.m[][])%p;
printf("%I64d\n",res);
}
return ;
}