【bzoj3240 && 洛谷P1397】矩阵游戏[NOI2013]（矩阵乘法+卡常）
题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=3240
这道题其实有普通快速幂+费马小定理的解法……然而我太弱了，一开始只想到了矩阵乘法的方法。
首先定义两个矩阵：
$ A_{1} = \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix} $
$ A_{2} = \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} $
于是我们就可以得到这样的式子：
$ \begin{aligned} \begin{bmatrix} f_{n,m} \\ 1 \end{bmatrix} & = A_{1} \begin{bmatrix} f_{n,m-1} \\ 1 \end{bmatrix} \\ & = A_{1}^{m-1} \begin{bmatrix} f_{n,1} \\ 1 \end{bmatrix} \\ & = A_{1}^{m-1} A_{2}\begin{bmatrix} f_{n-1,m} \\ 1 \end{bmatrix} \\ & = ( A_{1}^{m-1} A_{2} )^{n-1} \begin{bmatrix} f_{1,m} \\ 1 \end{bmatrix} \\ & = ( A_{1}^{m-1} A_{2} )^{n-1} A_{1}^{m-1} \begin{bmatrix} f_{1,1} \\ 1 \end{bmatrix} \end{aligned} $
然后用一发10进制矩阵快速幂能解决这道题了。
然而……这种做法跑的极慢。在洛谷上还能以近2000msAC，放到bzoj的6元cpu上跑就有些力不从心了，，所以得卡常数。
经过了十几次提交，使用了奥义·卡常数：10^18进制进制快速幂+循环展开+register后终于卡进了时限。。。
代码：
#include<cstdio> #include<cstring> #include<cmath> #include<cstdlib> #include<ctime> #include<algorithm> #include<queue> #include<vector> #define ll long long #define max(a,b) (a>b?a:b) #define min(a,b) (a<b?a:b) #define inf 0x3f3f3f3f #define mod 1000000007 #define base 1000000000000000000ll #define eps 1e-18 inline ll read() { ll tmp=0; char c=getchar(),f=1; for(;c<'0'||'9'<c;c=getchar())if(c=='-')f=-1; for(;'0'<=c&&c<='9';c=getchar())tmp=tmp*10+c-'0'; return tmp*f; } using namespace std; struct mat{ ll num[2][2]; }mat1,mat2; struct hp{ ll num[1000010]; int len; }n,m; char s[1000010],t[1000010]; ll a,b,c,d; mat times(mat a,mat b) { mat c; c.num[0][0]=(a.num[0][0]*b.num[0][0]+a.num[0][1]*b.num[1][0])%mod; c.num[0][1]=(a.num[0][0]*b.num[0][1]+a.num[0][1]*b.num[1][1])%mod; c.num[1][0]=(a.num[1][0]*b.num[0][0]+a.num[1][1]*b.num[1][0])%mod; c.num[1][1]=(a.num[1][0]*b.num[0][1]+a.num[1][1]*b.num[1][1])%mod; return c; } mat power_num(mat a,ll b) { mat ans; ans.num[0][0]=ans.num[1][1]=1; ans.num[0][1]=ans.num[1][0]=0; while(b){ if(b&1)ans=times(ans,a); a=times(a,a); b>>=1; } return ans; } mat power(mat a,hp b) { mat ans; ans.num[0][0]=ans.num[1][1]=1; ans.num[0][1]=ans.num[1][0]=0; for(register int i=1;i<=b.len;i++){ ans=times(ans,power_num(a,b.num[i])); a=power_num(a,base); } return ans; } int main() { register int i; scanf("%s",s); scanf("%s",t); a=read(); b=read(); c=read(); d=read(); mat1.num[0][0]=a; mat1.num[0][1]=b; mat1.num[1][0]=0; mat1.num[1][1]=1; mat2.num[0][0]=c; mat2.num[0][1]=d; mat2.num[1][0]=0; mat2.num[1][1]=1; int len1=strlen(s),len2=strlen(t); for(i=1;i*18<=len1;i++){ n.num[i]=0; for(register short j=18;j;j--) n.num[i]=n.num[i]*10+s[len1-(i-1)*18-j]-'0'; } n.len=len1/18; if(len1%18){ n.num[++n.len]=0; for(register short j=0;j<len1%18;j++)n.num[n.len]=n.num[n.len]*10+s[j]-'0'; } for(i=1;i*18<=len2;i++){ m.num[i]=0; for(register short j=18;j;j--) m.num[i]=m.num[i]*10+t[len2-(i-1)*18-j]-'0'; } m.len=len2/18; if(len1%18){ m.num[++m.len]=0; for(register short j=0;j<len2%18;j++)m.num[m.len]=m.num[m.len]*10+t[j]-'0'; } n.num[1]-=1; for(i=1;i<=n.len;i++)if(n.num[i]<0)n.num[i]+=base,n.num[i+1]-=1; if(n.num[n.len]==0&&n.len>1)--n.len; m.num[1]-=1; for(i=1;i<=m.len;i++)if(m.num[i]<0)m.num[i]+=base,m.num[i+1]-=1; if(m.num[m.len]==0&&m.len>1)--m.len; mat hang=power(mat1,m); mat ans=times(power(times(hang,mat2),n),hang); printf("%lld\n",(ans.num[0][0]+ans.num[0][1])%mod); }
解法2：
我们可以发现把$f_{1,i}=af_{i-1}+b$不断地展开后就能得到$f_{1,i}=a^i f_{1,1}+b\sum_{k=0}^{i-1}a^{k}$，即$f_{2,1}=a^{i}cf_{1,1}+bc\sum_{k=0}^{m-1}a^{k}+d$，于是我们可以用等比数列求和公式化简这个式子，又因为1e9+7是质数，所以可以由费马小定理将指数对1e9+6取模（不过需特判$a=1$的情况）。
然后我们可以发现这个式子也是$f_{i,1}=af_{i-1,1}+b$的形式（其中$ a,b $是常数），于是用上面的方法求出$f_{n+1,1}$，然后$f_{n,m}$就好求了。
代码：
#include<cstdio> #include<cstring> #include<cmath> #include<cstdlib> #include<ctime> #include<algorithm> #include<queue> #include<vector> #define ll long long #define max(a,b) (a>b?a:b) #define min(a,b) (a<b?a:b) #define inf 0x3f3f3f3f #define mod 1000000007 #define eps 1e-18 inline ll read() { ll tmp=0; char c=getchar(),f=1; for(;c<'0'||'9'<c;c=getchar())if(c=='-')f=-1; for(;'0'<=c&&c<='9';c=getchar())tmp=(tmp<<3)+(tmp<<1)+c-'0'; return tmp*f; } using namespace std; char s[1000010],t[1000010]; ll a,b,c,d,n,m; ll power(ll a,ll b) { ll ans=1; while(b){ if(b&1)ans=ans*a%mod; a=a*a%mod; b>>=1; } return ans; } int main() { int i; scanf("%s",s); scanf("%s",t); a=read(); b=read(); c=read(); d=read(); int len1=strlen(s),len2=strlen(t); m=0; for(i=0;i<len2;i++)m=(m*10+t[i]-'0')%(a>1?mod-1:mod); ll p=power(a,m-1)*c%mod,q=((a>1?(power(a,m-1)+mod-1)*power(a+mod-1,mod-2)%mod:m-1)*b%mod*c%mod+d)%mod; for(i=0;i<len1;i++)n=(n*10+s[i]-'0')%(p>1?mod-1:mod); ll ans=(power(p,n)+(p>1?(power(p,n)+mod-1)*power(p+mod-1,mod-2)%mod:n)*q)%mod; printf("%lld\n",(ans+mod-d)*power(c,mod-2)%mod); }
