BZOJ3240 NOI2013 矩阵游戏 快速幂+矩阵乘法+逆元

题意：给定a b c d n m，f_i,j满足f_1,1=1  f_i,j=af_i,j-1+b  f_i,1=cf_i-1,m+d，求f_n,m%P

题解：

a=1且c=1时\[{f_{n,m}} = 1 + b(m - 1)n + d(n - 1)\]否则\[\left( {\begin{array}{*{20}{c}}
{{f_{n,m}}}&{bd}
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{f_{1,1}}}&{bd}
\end{array}} \right){\left( {\begin{array}{*{20}{c}}
a&0\\
{\frac{1}{d}}&1
\end{array}} \right)^{nm - n}}{\left( {\begin{array}{*{20}{c}}
c&0\\
{\frac{1}{b}}&1
\end{array}} \right)^{n - 1}}\]根据费马小定理，直接把指数压倒≤P

#include <cmath>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <climits>
#include <iostream>
#include <algorithm>
using namespace std;
#define P 1000000007
#define ll long long

struct MATRIX;

void Output(MATRIX a);

ll Quick_Mul(ll a,ll b){
    a%=P,b%=P;
    ll t=(b&1?a:0);
    while(b>>=1){
        a=(a+a)%P;
        if(b&1) t=(t+a)%P;
    }
    return t;
}

struct MATRIX{
    ll m[6][6];
    int l,r;
    MATRIX(){}
    MATRIX(int _l,int _r):l(_l),r(_r){ memset(m,0,sizeof(m));}
    MATRIX operator*(MATRIX a){
        MATRIX t=MATRIX(l,a.r);
        for(int i=1;i<=l;i++)
        for(int j=1;j<=a.r;j++)
        for(int k=1;k<=r;k++)
            t.m[i][j]=(t.m[i][j]+Quick_Mul(m[i][k],a.m[k][j]))%P;
        return t;
    }
    MATRIX operator^(ll x){
        MATRIX t=MATRIX(l,r),a=MATRIX(l,r);
        for(int i=1;i<=l;i++)
            for(int j=1;j<=r;j++)
                a.m[i][j]=m[i][j];

        if(x&1)
            for(int i=1;i<=l;i++)
                for(int j=1;j<=r;j++)
                    t.m[i][j]=m[i][j];
        else
            for(int i=1;i<=l;i++) t.m[i][i]=1;

        while(x>>=1){
            a=a*a;
            if(x&1) t=t*a;
        }
        return t;
    }
}a=MATRIX(1,2),b=MATRIX(2,2),c=MATRIX(2,2);

const int MAXN=1000000+2;
char S1[MAXN],S2[MAXN];
ll N,M,A,B,C,D;

void Output(MATRIX a){
    for(int i=1;i<=a.l;i++){
        for(int j=1;j<=a.r;j++)
            cout << a.m[i][j] << " ";
        cout << endl;
    }
    cout << endl;
}

void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,x,y);

    ll t=x;
    x=y,y=t-a/b*y;
}

ll Calc_Inv(ll a){
    ll x,y;
    exgcd(P,a,x,y);
    return (y+P)%P;
}

ll Calc(char *S,int p){
    int l=strlen(S);
    ll ret=0;
    for(int i=0;i<l;i++)
        ret=((ll)ret*10+S[i]-'0')%p;
    return ret;
}

int main(){
    scanf("%s %s %lld %lld %lld %lld",S1,S2,&A,&B,&C,&D);

    if(A==1){
        N=Calc(S1,P),M=Calc(S2,P);
        printf("%lld\n",((1+Quick_Mul(B,Quick_Mul(M-1,N)))+Quick_Mul(D,N-1))%P);
        return 0;
    }

    N=Calc(S1,P-1),M=Calc(S2,P-1);
    a.m[1][1]=1,a.m[1][2]=B*D%P;
    b.m[1][1]=A,b.m[1][2]=0,b.m[2][1]=Calc_Inv(D),b.m[2][2]=1;
    c.m[1][1]=C,c.m[1][2]=0,c.m[2][1]=Calc_Inv(B),c.m[2][2]=1;

    b=b^(M-1);
    a=a*(((b*c)^(N-1))*b);

    printf("%lld\n",a.m[1][1]);

    return 0;
}