[题解]bzoj3240(NOI2013)矩阵游戏

Description

婷婷是个喜欢矩阵的小朋友,有一天她想用电脑生成一个巨大的n行m列的矩阵(你不用担心她如何存储)。她生成的这个矩阵满足一个神奇的性质：若用F[i][j]来表示矩阵中第i行第j列的元素，则F[i][j]满足下面的递推式:

F[1][1]=1

F[i,j]=a*F[i][j-1]+b (j!=1)

F[i,1]=c*F[i-1][m]+d (i!=1)

递推式中a,b,c,d都是给定的常数。

现在婷婷想知道F[n][m]的值是多少,请你帮助她。由于最终结果可能很大，你只需要输出F[n][m]除以1,000,000,007的余数。

Input

一行有六个整数n,m,a,b,c,d。意义如题所述

Output

包含一个整数，表示F[n][m]除以1,000,000,007的余数

Sample Input

3 4 1 3 2 6

Sample Output

85

HINT

样例中的矩阵为：

1 4 7 10

26 29 32 35

76 79 82 85

1<=N,M<=10^1000 000,1<=a,b,c,d<=10^9

Solution

其实就几个线性递推在一起。矩阵不好写，我们考虑推通项公式。

有：（标准形式）

$$f_i=a*f_{i-1}+b$$

那么：(可以自己推一下，很容易得到的结论。)

$$f_m=f_1*a^{m-1}+b*\sum_{i=0}^{m-2}a^i$$

所以如果我们知道了题目中每一行的第一个数字，那么我们就可以很快求出该行第m个。现在的问题转化为推出第n行第一个数字。

怎么做呢？我们考虑：

$$f_{i,1}=c*f_{i-1,m}+d$$

将之前的结论带入，得：

$$f_{i,1}=c*[f_{i-1,1}*a^{m-1}+b*\sum_{j=0}^{m-2}a^j]+d$$

$$f_{i,1}=c*a^{m-1}*f_{i-1,1}+b*c*\sum_{j=0}^{m-2}a^j+d$$

令$a'=c*a^{m-1}$，$b'=b*c*\sum_{j=0}^{m-2}a^j+d$

于是第一列也转化成了标准形式：

$$f_{i,1}=a'*f_{i-1,1}+b'$$

我们依旧可以很快求出我们想要的$f_{n,1}$，然后求出$f_{n,m}$。

现在还有一些细节上的处理问题。首先我们要能快速求$a^b$，其中b为高精度整数。我们考虑把b按十进制位分解，预处理出$a^{10^i}$,然后就可以扫一遍b，看看他的第i位上是几，是几就把答案乘上几个$a^{10^i}$。

同时$a^{m-1}$可以看作$a^m/a$，所以还需要求a在mod $10^9+7$下的逆元。由于模数是比a大的质数，所以a和模数互质，根据费马小定理可知a的逆元是$a^{10^9+5}$，直接快速幂。

还有一个问题就是怎么快速求$\sum_{i=0}^{m-2}a^i$。

我们令$A=\sum_{i=0}^{m-2}a^i$，那么：

$$a*A=\sum_{i=1}^{m-1}a^i$$

所以

$$a*A-A=a^{m-1}-1$$

即

$$A=(a^{m-1}-1)/(a-1)$$

用快速幂和逆元求即可。

这样，问题完美解决～

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;typedef long long LL;
const int mod=1000000007,maxn=1000010;
int a,b,c,d;
LL e,f,g,ans;
int n[maxn],m[maxn],powa[maxn],powe[maxn];
char s1[maxn],s2[maxn];
int pow(int x,int y){int re=1;for(;y;y>>=1,x=1ll*x*x%mod){if(y&1){re=1ll*re*x%mod;}}return re;
}
int pow(int x,int y[],int p[]){int re=1;for(int i=1;i<=y[0];i++){for(int j=1;j<=y[i];j++){re=1ll*re*p[i]%mod;}}return re;
}
int inv(int x){return pow(x,mod-2);
}
int getint(int x[]){int re=0;for(int i=x[0];i>=1;i--){re=(1ll*re*10%mod+x[i])%mod;}return re;
}int main(){scanf("%s%s%d%d%d%d",s1,s2,&a,&b,&c,&d);n[0]=strlen(s1);m[0]=strlen(s2);for(int i=1;i<=n[0];i++)n[i]=s1[n[0]-i]-'0';for(int i=1;i<=m[0];i++)m[i]=s2[m[0]-i]-'0';if(a!=1){powa[1]=a;for(int i=2;i<=m[0];i++){powa[i]=pow(powa[i-1],10);}int temp=1ll*pow(a,m,powa)*inv(a)%mod;e=1ll*temp*c%mod;f=(1ll*b*((temp-1+mod)%mod)%mod*inv(a-1)%mod*c%mod+d)%mod;}else{e=c;f=(1ll*b*c%mod*(getint(m)-1)%mod+d)%mod;}if(e!=1){powe[1]=e;for(int i=2;i<=n[0];i++){powe[i]=pow(powe[i-1],10);}int temp=1ll*pow(e,n,powe)*inv(e)%mod;g=(1ll*f*((temp-1+mod)%mod)%mod*inv(e-1)%mod+temp)%mod;}else{g=(1ll*f*(getint(n)-1)%mod+1)%mod;}if(a!=1){int temp=1ll*pow(a,m,powa)*inv(a)%mod;ans=(1ll*b*((temp-1+mod)%mod)%mod*inv(a-1)%mod+temp*g%mod)%mod;}else{ans=(1ll*b*(getint(m)-1)%mod+g)%mod;}printf("%lld\n",ans);return 0;
}