题目传送门：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=; char c=getchar(),f=;
for(;c<''||''<c;c=getchar())if(c=='-')f=-;
for(;''<=c&&c<='';c=getchar())tmp=tmp*+c-'';
return tmp*f;
}
using namespace std;
struct mat{
ll num[][];
}mat1,mat2;
struct hp{
ll num[];
int len;
}n,m;
char s[],t[];
ll a,b,c,d;
mat times(mat a,mat b)
{
mat c;
c.num[][]=(a.num[][]*b.num[][]+a.num[][]*b.num[][])%mod;
c.num[][]=(a.num[][]*b.num[][]+a.num[][]*b.num[][])%mod;
c.num[][]=(a.num[][]*b.num[][]+a.num[][]*b.num[][])%mod;
c.num[][]=(a.num[][]*b.num[][]+a.num[][]*b.num[][])%mod;
return c;
}
mat power_num(mat a,ll b)
{
mat ans; ans.num[][]=ans.num[][]=; ans.num[][]=ans.num[][]=;
while(b){
if(b&)ans=times(ans,a);
a=times(a,a); b>>=;
}
return ans;
}
mat power(mat a,hp b)
{
mat ans; ans.num[][]=ans.num[][]=; ans.num[][]=ans.num[][]=;
for(register int i=;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[][]=a; mat1.num[][]=b; mat1.num[][]=; mat1.num[][]=;
mat2.num[][]=c; mat2.num[][]=d; mat2.num[][]=; mat2.num[][]=;
int len1=strlen(s),len2=strlen(t);
for(i=;i*<=len1;i++){
n.num[i]=;
for(register short j=;j;j--)
n.num[i]=n.num[i]*+s[len1-(i-)*-j]-'';
}
n.len=len1/;
if(len1%){
n.num[++n.len]=;
for(register short j=;j<len1%;j++)n.num[n.len]=n.num[n.len]*+s[j]-'';
}
for(i=;i*<=len2;i++){
m.num[i]=;
for(register short j=;j;j--)
m.num[i]=m.num[i]*+t[len2-(i-)*-j]-'';
}
m.len=len2/;
if(len1%){
m.num[++m.len]=;
for(register short j=;j<len2%;j++)m.num[m.len]=m.num[m.len]*+t[j]-'';
}
n.num[]-=;
for(i=;i<=n.len;i++)if(n.num[i]<)n.num[i]+=base,n.num[i+]-=;
if(n.num[n.len]==&&n.len>)--n.len;
m.num[]-=;
for(i=;i<=m.len;i++)if(m.num[i]<)m.num[i]+=base,m.num[i+]-=;
if(m.num[m.len]==&&m.len>)--m.len;
mat hang=power(mat1,m);
mat ans=times(power(times(hang,mat2),n),hang);
printf("%lld\n",(ans.num[][]+ans.num[][])%mod);
}
bzoj3240
解法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=; char c=getchar(),f=;
for(;c<''||''<c;c=getchar())if(c=='-')f=-;
for(;''<=c&&c<='';c=getchar())tmp=(tmp<<)+(tmp<<)+c-'';
return tmp*f;
}
using namespace std;
char s[],t[];
ll a,b,c,d,n,m;
ll power(ll a,ll b)
{
ll ans=;
while(b){
if(b&)ans=ans*a%mod;
a=a*a%mod; b>>=;
}
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=;
for(i=;i<len2;i++)m=(m*+t[i]-'')%(a>?mod-:mod);
ll p=power(a,m-)*c%mod,q=((a>?(power(a,m-)+mod-)*power(a+mod-,mod-)%mod:m-)*b%mod*c%mod+d)%mod;
for(i=;i<len1;i++)n=(n*+s[i]-'')%(p>?mod-:mod);
ll ans=(power(p,n)+(p>?(power(p,n)+mod-)*power(p+mod-,mod-)%mod:n)*q)%mod;
printf("%lld\n",(ans+mod-d)*power(c,mod-)%mod);
}
bzoj3240