题目链接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 都为
⇒[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上的数据
这里贴代码(代码有一些奇奇怪怪的东西 当然你可以无视)
#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 <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)都存下来,不要先存字符串。这样的话会快一些(纯粹是为了刷个排名吧)代码这里就不贴了,留给读者自己思考。