【矩阵+十进制快速幂】[NOI2013]矩阵游戏

题目描述

婷婷是个喜欢矩阵的小朋友,有一天她想用电脑生成一个巨大的 n 行 m 列的矩阵（你不用担心她如何存储）。她生成的这个矩阵满足一个神奇的性质：
若用 Fi,j 来表示矩阵中第 i 行第 j 列的元素，则 Fi,j 满足下面的递推式：

⎧⎩⎨⎪⎪F1,1Fi,jFi,1===1a⋅Fi,j−1+bc⋅Fi−1,m+dj≠1i≠1

递推式中 a,b,c,d 都是给定的常数。 现在婷婷想知道 Fn,m 的值是多少，请你帮助她。由于最终结果可能很大，你只需要输出 Fn,m 除以 1000000007 的余数。

输入格式

包含一行有六个整数 n,m,a,b,c,d。意义如题所述。

输出格式

表示 Fn,m 除以 1000000007 的余数。

样例一

input

3 4 1 3 2 6

output

85

explanation

样例中的矩阵为：

⎛⎝⎜126764297973282103585⎞⎠⎟

样例二

见样例数据下载。

限制与约定

测试点编号	数据范围
1	1≤n,m≤10；1≤a,b,c,d≤1000
2	1≤n,m≤100；1≤a,b,c,d≤1000
3	1≤n,m≤103；1≤a,b,c,d≤109
4	1≤n,m≤103；1≤a,b,c,d≤109
5	1≤n,m≤109；1≤a=c≤109；1≤b=d≤109
6	1≤n,m≤109；a=c=1；1≤b,d≤109
7	1≤n,m,a,b,c,d≤109
8	1≤n,m,a,b,c,d≤109
9	1≤n,m,a,b,c,d≤109
10	1≤n,m,a,b,c,d≤109
11	1≤n,m≤101000；a=c=1；1≤b,d≤109
12	1≤n,m≤101000；1≤a=c≤109；1≤b=d≤109
13	1≤n,m≤101000；1≤a,b,c,d≤109
14	1≤n,m≤101000；1≤a,b,c,d≤109
15	1≤n,m≤1020000；1≤a,b,c,d≤109
16	1≤n,m≤1020000；1≤a,b,c,d≤109
17	1≤n,m≤101000000；a=c=1；1≤b,d≤109
18	1≤n,m≤101000000；1≤a=c≤109；1≤b=d≤109
19	1≤n,m≤101000000；1≤a,b,c,d≤109
20	1≤n,m≤101000000；1≤a,b,c,d≤109

时间限制：1s

空间限制：256MB

分析

可以很容易地推出矩阵

每一行

(Fi,j+11)=(Fi,j1)×(ab01)

行之间

(Fi+1,11)=(Fi,m1)×(cd01)

然后可以快速幂算出Fi,1和Fi,m，然后算出和Fi+1,1之间的矩阵。然后继续快速幂就好了。

由于指数是高精度数，所以使用十进制快速幂算比较方便。

代码

ps:高精度粘的版，很长。。。

#include<cstdio>
#include<algorithm>
#include<cstring>
#define MOD 1000000007
#define MAXLEN 1000000
using namespace std;
int a,b,c,d;
char s[MAXLEN+10];
void Read(int &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
void Read_mod(int &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=(1ll*x*10+c-'0')%(MOD-1);
            ungetc(c,stdin);
            return;
        }
}
struct Matrix{
    int a[2][2]; // f 1
    inline Matrix(){
        memset(a,0,sizeof a);
    }
    inline Matrix(int){
        memset(a,0,sizeof a);
        for(int i=0;i<2;i++)
            a[i][i]=1;
    }
    inline Matrix operator*(const Matrix &b)const{
        Matrix c;
        int i,j,k;
        for(i=0;i<2;i++)
            for(j=0;j<2;j++)
                if(a[i][j])
                    for(k=0;k<2;k++)
                        c.a[i][k]=(c.a[i][k]+1ll*a[i][j]*b.a[j][k])%MOD;
        return c;
    }
    inline Matrix &operator*=(const Matrix &b){
        return *this=*this*b;
    }
}A,B,C;
struct hp{
    int a[MAXLEN+10];
    hp(){
        memset(a,0,sizeof a);
        a[0]=1;
    }
    hp(int n){
        memset(a,0,sizeof a);
        a[0]=0;
        while(n){
            a[++a[0]]=n%10;
            n/=10;
        }
        if(!a[0])
            a[0]=1;
    }
    hp(char *s){
        memset(a,0,sizeof a);
        int len=strlen(s);
        for(int i=1;i<=len;i++)
            a[i]=s[len-i]-'0';
        a[0]=len;
    }
    hp operator*(const hp &b)const{
        hp c;
        int i,j,len=a[0]+b.a[0];
        for(i=1;i<=a[0];i++)
            for(j=1;j<=b.a[0];j++)
                c.a[i+j-1]+=a[i]*b.a[j];
        for(i=1;i<len;i++){
            c.a[i+1]+=c.a[i]/10;
            c.a[i]%=10;
        }
        while(len>1&&!c.a[len])
            len--;
        c.a[0]=len;
        return c;
    }
    hp operator/(int b)const{
        hp c;
        int d=0,i,len=a[0];
        for(i=a[0];i;i--){
            d=d*10+a[i];
            c.a[i]=d/b;
            d%=b;
        }
        while(len>1&&!c.a[len])
            len--;
        c.a[0]=len;
        return c;
    }
    hp operator+(const hp &b)const{
        hp c;
        int len=max(a[0],b.a[0]),i;
        for(i=1;i<=len;i++){
            c.a[i]+=a[i]+b.a[i];
            c.a[i+1]=c.a[i]/10;
            c.a[i]%=10;
        }
        len++;
        while(len>1&&!c.a[len])
            len--;
        c.a[0]=len;
        return c;
    }
    hp operator-(const hp &b)const{
        hp c;
        int i,len=a[0];
        for(i=1;i<=len;i++){
            c.a[i]+=a[i]-b.a[i];
            if(c.a[i]<0)
                c.a[i]+=10,c.a[i+1]--;
        }
        while(len>1&&!c.a[len])
            len--;
        c.a[0]=len;
        return c;
    }
    void operator*=(const hp &x){
        *this=*this*x;
    }
    void operator/=(const int &x){
        *this=*this/x;
    }
    void operator+=(const hp &x){
        *this=*this+x;
    }
    void operator-=(const hp &x){
        *this=*this-x;
    }
    void print(){
        for(int i=a[0];i;i--)
            printf("%d",a[i]);
    }
    bool operator>(const hp&b)const{
        if(a[0]>b.a[0])
            return 1;
        if(a[0]<b.a[0])
            return 0;
        for(int i=a[0];i;i--)
            if(a[i]>b.a[i])
                return 1;
            else if(a[i]<b.a[i])
                return 0;
        return 0;
    }
    bool operator<(const hp&b)const{
        if(a[0]<b.a[0])
            return 1;
        if(a[0]>b.a[0])
            return 0;
        for(int i=a[0];i;i--)
            if(a[i]<b.a[i])
                return 1;
            else if(a[i]>b.a[i])
                return 0;
        return 0;
    }
    bool operator<=(const hp&b)const{
        return !(*this>b);
    }
    hp operator/(const hp&b)const{
        hp l(0),r(*this),mid;
        while(l<r){
            mid=(l+r+1)/2;
            if(mid*b<=*this)
                l=mid;
            else
                r=mid-1;
        }
        return l;
    }
    void operator/=(const hp&b){
        *this=*this/b;
    }
}n,m;
template<class T>
T quick_pow(T a,int b){
    T ret(1);
    while(b){
        if(b&1)
            ret*=a;
        a*=a;
        b>>=1;
    }
    return ret;
}
void read(){
    scanf("%s",s);
    n=s;
    scanf("%s",s);
    m=s;
    Read(a),Read(b),Read(c),Read(d);
}
template<class T>
T quick_pow10(T a,hp b){
    int i;
    T ret(1),c[10];
    c[0]=ret;
    for(i=1;i<10;i++)
        c[i]=c[i-1]*a;
    for(i=b.a[0];i;i--){
        ret=quick_pow(ret,10);
        ret*=c[b.a[i]];
    }
    return ret;
}
void solve(){
    A.a[0][0]=a,A.a[1][0]=b,A.a[1][1]=1;
    B=quick_pow10(A,m-1);
    C.a[0][0]=c,C.a[1][0]=d,C.a[1][1]=1;
    A=B*C;
    A=quick_pow10(A,n-1);
    A*=B;
}
int main()
{
    read();
    solve();
    printf("%d\n",(A.a[0][0]+A.a[1][0])%MOD);
}