NOI2013矩阵游戏

题目链接LUOGU1397矩阵游戏
题目大意
已知

$$\begin{cases} f_{1,1}=1 \\ f_{i,j}=a*f_{i,j-1}+b &\text[j!=1]\\ f_{i,1}=c*f_{i-1,m}+d&\text[i\neq1] \end{cases}$$
求 $f_{n,m}$

算法一 矩阵快速幂

我们可以构造矩阵

$$\begin{bmatrix} 1&f_{n,1} \end{bmatrix} \times \begin{bmatrix} 1&b\\ 0&a \end{bmatrix}^{m-1} = \begin{bmatrix} 1&f_{n,m} \end{bmatrix}$$
$\Rightarrow$
$$\begin{bmatrix} 1&c*f_{n-1,m}+d \end{bmatrix} \times { \begin{bmatrix} 1&b\\ 0&a \end{bmatrix} }^{m-1} = \begin{bmatrix} 1&f_{n,m} \end{bmatrix}$$
$\Rightarrow$
$$令T= \begin{bmatrix} 1&b\\ 0&a \end{bmatrix}^{m-1}$$
$$x=T_{2,2}$$
$$y=T_{1,2}$$
$$∴f_{n,m}=x*c*f_{n-1,m}+x*d+y$$
$\Rightarrow$
$$\begin{bmatrix} 1&f_{1,m} \end{bmatrix} \times { \begin{bmatrix} 1&x*d+y\\ 0&x*c \end{bmatrix} }^{n-1} = \begin{bmatrix} 1&f_{n,m} \end{bmatrix}$$
所以我们可以先快速幂求出 $f_{1,m}$然后再快速幂求出 $f_{n,m}$
但是由于 $n,m$很大 $n,m \le 10^{1000000}$我们用高精度显然超时,而且非常不好打
因为 $10^9+7$是一个质数,所以我们考虑费马小定理
$$a^{p-1} \equiv 1 \pmod p \Rightarrow a^{n} \equiv a^{n \pmod {p-1}}\pmod p$$
因为矩阵乘法也是满足费马小定理的
但是有一点问题(这里只举例说 $2\times2$的矩阵 $n\times n$的矩阵也是一样的)
$$\begin{bmatrix} a&b\\ c&d \end{bmatrix}^{p-1} \equiv \begin{bmatrix} 1&0\\ 0&1 \end{bmatrix} \text{if a$\neq$d} \pmod p \\[2ex] \begin{bmatrix} a&b\\ c&d \end{bmatrix}^{p} \equiv \begin{cases} \begin{bmatrix} a&0\\ 0&a \end{bmatrix}&\text{if a=d} \\[2ex] \begin{bmatrix} a&b\\ c&d \end{bmatrix}&\text{if a$\neq$d} \end{cases} \pmod p$$
又在这道题中两个矩阵的 $a$都为$1$
$$\Rightarrow \begin{bmatrix} 1&f_{n,1} \end{bmatrix} \times \begin{bmatrix} 1&b\\ 0&a \end{bmatrix}^{m-1} = \begin{bmatrix} 1&f_{n,1} \end{bmatrix} \times \begin{cases} \begin{bmatrix} 1&b\\ 0&a \end{bmatrix}^{m-1 \pmod p}&\text{if a=1}\\[2ex] \begin{bmatrix} 1&b\\ 0&a \end{bmatrix}^{m-1 \pmod {p-1}}&\text{if a$\neq$1}\\[2ex] \end{cases}$$
所以我们可以直接取模,对于第二个矩阵快速幂也要判断一下 $x*c$是否等于 $1$
否则会被UOJ上的数据$Hack$掉(当然你第二个不判断也可以过官方数据,第一个不判你也有85分)

这里贴代码（代码有一些奇奇怪怪的东西 当然你可以无视）

#include<cstdio>
#include<cstring>
#define re register int
#define fp(i,a,b) for(re i=a,I=b;i<=I;++i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
char ss[1<<17],*A=ss,*B=ss,ch;
inline char gc(){if(A==B){B=(A=ss)+fread(ss,1,1<<17,stdin);if(A==B)return EOF;}return*A++;}
template<class T>inline void sdf(T&x){
    while(ch=gc(),ch<48);x=ch^48;
    while(ch=gc(),48<=ch)x=10*x+ch-48;
}
inline void cis(char*s){
    while(ch=gc(),ch<48);*s++=ch;
    while(ch=gc(),48<=ch)*s++=ch;
}
char sr[1<<20],z[20];int C=-1,Z;
template<class T>inline void wer(T x){
    while(z[++Z]=x%10+'0',x/=10);
    while(sr[++C]=z[Z],--Z);
}
const int N=3,P=1e9+7;
typedef int arr[N];
typedef long long ll;
struct matrix{
    int a[N][N];
    inline int*operator[](re x){return a[x];}
    matrix(re x=0){memset(a,0,sizeof a);if(x==1)fp(i,0,N-1)a[i][i]=1;}
    inline matrix operator*(matrix b)const{
        matrix c;
        fp(k,1,2)fp(i,1,2)if(a[i][k])fp(j,1,2)
            c[i][j]=(c[i][j]+1ll*a[i][k]*b[k][j])%P;
        return c;
    }
    matrix operator^(ll b){
        matrix x(1),A=*this;
        for(;b;b>>=1,A=A*A)if(b&1)x=x*A;
        return x;
    }
}T;
ll ans=1,n,m,a,b,c,d,mod=P-1;char s1[1000010],s2[1000010];
int main(){
    file("oi");
    cis(s1);cis(s2);sdf(a);sdf(b);sdf(c);sdf(d);mod+=a==1;
    fp(i,0,strlen(s2)-1)m=(m*10+s2[i]-48)%mod;--m;if(m<0)m+=mod;
    T[2][2]=a;T[1][2]=b;T[1][1]=1;T=T^m;
    ans=(T[1][2]+T[2][2])%P;mod=P-1;
    T[1][2]=(1ll*T[2][2]*d+T[1][2])%P;
    T[2][2]=1ll*T[2][2]*c%P;mod+=T[2][2]==1;
    fp(i,0,strlen(s1)-1)n=(n*10+s1[i]-48)%mod;--n;if(n<0)n+=mod;
    T=T^n;ans=(T[1][2]+1ll*T[2][2]*ans)%P;wer(ans);
    fwrite(sr,1,C+1,stdout);
return 0;
}

算法二 数列通项公式

对于这种一阶递推式我们完全可以用等比数列的通项公式来求（反正上必修五老师会讲怎求的）注意接下来的运算都是带取模的,模数$P=10^9+7$

$$f_{n,m}=a*f_{n,m-1}+b\\ f_{n,m}=\begin{cases} f_{n,1}+b*(m-1)&\text{if a=1}\\ a^{m-1}*f_{n,1}+(a^{m-1}-1)*\frac{b}{a-1}&\text{if a$\neq$1} \end{cases}\\ 令x=\begin{cases} 1&\text{if a=1}\\ a^{m-1}&\text{if a$\neq$1} \end{cases}\\ y=\begin{cases} b*(m-1)&\text{if a=1}\\ (a^{m-1}-1)*\frac{b}{a-1}&\text{if a$\neq$1} \end{cases}\\ \Rightarrow f_{1,m}=x*f_{1,1}+y\\ \Rightarrow f_{n,m}=x*(c*f_{n-1,m}+d)+y=x*c*f_{n-1,m}+x*d+y\\ 令p=x*c\\q=x*d+y\\ \Rightarrow f_{n,m}=p*f_{n,m-1}+q= \begin{cases} f_{1,n}+q*(n-1)&\text{if p=1}\\ p^{n-1}*f_{1,m}+(p^{n-1}-1)*\frac{q}{p-1}&\text{if p$\neq$1} \end{cases}\\ $$
同样的道理我们要用费马小定理来缩小 $n,m$
处理方式和矩阵的一样要分情况讨论
因为当 $a=1$或 $p=1$是这是等差数列而 $\gcd (a,P)=1,\gcd (p,P)=1$所以循环长度是 P <script type="math/tex" id="MathJax-Element-34">P</script>
反之,我们就可以直接用快速幂运算了
这里贴代码

#include<cstdio>
#include<cstring>
#define re register int
#define rg register long long
#define fp(i,a,b) for(re i=a,I=b;i<=I;++i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
char ss[1<<17],*A=ss,*B=ss,ch;
inline char gc(){if(A==B){B=(A=ss)+fread(ss,1,1<<17,stdin);if(A==B)return EOF;}return*A++;}
template<class T>inline void sdf(T&x){
    while(ch=gc(),ch<48);x=ch^48;
    while(ch=gc(),48<=ch)x=10*x+ch-48;
}
inline void cis(char*s){
    while(ch=gc(),ch<48);*s++=ch;
    while(ch=gc(),48<=ch)*s++=ch;
}
char sr[1<<20],z[20];int C=-1,Z;
template<class T>inline void wer(T x){
    while(z[++Z]=x%10+'0',x/=10);
    while(sr[++C]=z[Z],--Z);
}
const int N=3,P=1e9+7;
typedef int arr[N];
typedef long long ll;
ll ans,n,m,a,b,c,d,mod;char s1[1000010],s2[1000010];
inline ll powMod(rg a,rg b){
    rg x=1;
    for(;b;b>>=1,a=(a*a)%P)if(b&1)x=(x*a)%P;
    return x;
}
int main(){
    cis(s1);cis(s2);sdf(a);sdf(b);sdf(c);sdf(d);mod=P-(a!=1);
    fp(i,0,strlen(s2)-1)m=(m*10+s2[i]-48)%mod;--m;if(m<0)m+=mod;
    rg x=a==1?1:powMod(a,m)%P,y=a==1?b*m%P:(powMod(a,m)-1+P)*b%P*powMod(a-1,P-2)%P;
    ans=x+y;y=(x*d+y)%P;x=x*c%P;mod=P-(x!=1);
    fp(i,0,strlen(s1)-1)n=(n*10+s1[i]-48)%mod;--n;if(n<0)n+=mod;
    if(x==1)ans=(ans+n*y%P)%P;
    else ans=(ans*powMod(x,n)+(powMod(x,n)-1+P)*y%P*powMod(x-1,P-2))%P;
    wer(ans);
    fwrite(sr,1,C+1,stdout);
return 0;
}

这两份代码的复杂度都是一样的，主要卡在读入上。读者可以考虑把n%P和n%(P-1)都存下来，不要先存字符串。这样的话会快一些（纯粹是为了刷个排名吧）代码这里就不贴了，留给读者自己思考。