NOI2013矩阵游戏

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

⎧⎩⎨⎪⎪f1,1=1fi,j=a∗fi,j−1+bfi,1=c∗fi−1,m+d[j!=1][i≠1]

求fn,m

算法一 矩阵快速幂

我们可以构造矩阵

[1fn,1]×[10ba]m−1=[1fn,m]

⇒

[1c∗fn−1,m+d]×[10ba]m−1=[1fn,m]

⇒

令T=[10ba]m−1

x=T2,2

y=T1,2

∴fn,m=x∗c∗fn−1,m+x∗d+y

⇒

[1f1,m]×[10x∗d+yx∗c]n−1=[1fn,m]

所以我们可以先快速幂求出f1,m然后再快速幂求出fn,m
但是由于n,m很大n,m≤101000000我们用高精度显然超时,而且非常不好打
因为109+7是一个质数,所以我们考虑费马小定理

ap−1≡1(modp)⇒an≡an(modp−1)(modp)

因为矩阵乘法也是满足费马小定理的
但是有一点问题(这里只举例说2×2的矩阵 n×n的矩阵也是一样的)

[acbd]p−1≡[1001]if a≠d(modp)[acbd]p≡⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪[a00a][acbd]if a=dif a≠d(modp)

又在这道题中两个矩阵的a都为1

⇒[1fn,1]×[10ba]m−1=[1fn,1]×⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪[10ba]m−1(modp)[10ba]m−1(modp−1)if a=1if a≠1

所以我们可以直接取模,对于第二个矩阵快速幂也要判断一下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=109+7

fn,m=a∗fn,m−1+bfn,m={fn,1+b∗(m−1)am−1∗fn,1+(am−1−1)∗ba−1if a=1if a≠1令x={1am−1if a=1if a≠1y={b∗(m−1)(am−1−1)∗ba−1if a=1if a≠1⇒f1,m=x∗f1,1+y⇒fn,m=x∗(c∗fn−1,m+d)+y=x∗c∗fn−1,m+x∗d+y令p=x∗cq=x∗d+y⇒fn,m=p∗fn,m−1+q={f1,n+q∗(n−1)pn−1∗f1,m+(pn−1−1)∗qp−1if p=1if p≠1

同样的道理我们要用费马小定理来缩小n,m
处理方式和矩阵的一样要分情况讨论
因为当a=1或p=1是这是等差数列而gcd(a,P)=1,gcd(p,P)=1所以循环长度是P
反之,我们就可以直接用快速幂运算了
这里贴代码

#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)都存下来，不要先存字符串。这样的话会快一些（纯粹是为了刷个排名吧）代码这里就不贴了，留给读者自己思考。