矩阵快速幂

矩阵运算

矩阵加法

矩阵加法是由两个行列数相等的矩阵进行的运算，结果也为同样的行列数。其具体运算方式为两矩阵对应位置的数相加，放入答案矩阵的相同位置里。

例如：

依据这个运算方式可以写出代码：

struct s{
    int a[105][105];
    int r, c;
    s(int _n=0,int _m= 0) {
        n=_n,m=_m;
        memset(a,0,sizeof a);
        if(!m)m=n;
    }
    s operator +(const s t) const {
        s ans(n,m);
        for (int i=1;i<=n;i++)
            for (int j=1;j<=m;j++)ans.a[i][j]=a[i][j]+t.a[i][j];return ans;
    }
};

矩阵数乘

矩阵数乘是由一个数字乘上一个矩阵。运算方式为矩阵中每个数分别乘上数字。

例如：

矩阵乘法

矩阵乘法相对较为复杂。首先假设有两个矩阵，行列数分别为n1行m1列和n2行m2列，只有n1=m2时，两矩阵才能做乘法运算。结果为一个n2行m1列的矩阵。

运算方法如下：

首先将第一个矩阵第一行的每一个元素乘上第二个矩阵第一列的每一个元素再求和。将结果作为结果矩阵的第一行第一列。

例如：

按照此步骤结果为1*1+2*3=7

将 7放入结果矩阵。

然后再用第一行乘上第二列，以此类推，直到将m2列全部乘完。此时结果矩阵有一行n2个元素。

接着用第二行继续从第二个矩阵第一行乘上第二列，以此类推，直到把第一个矩阵每一行都乘完。

例如上面举例子的矩阵，结果为

由此可以写出代码：

s operator*(const s t)const{
        s ans(n,t.m);
        int r=n,c=t.m;
        for (int i=1;i<=r;i++)
            for (int j=1;j<=c;j++)
                for (int k=1;k<=m;k++)ans.a[i][j]+=a[i][k]*t.a[k][j];return ans;
}

矩阵快速幂

单位元

当在一个运算系统中，对于运算Δ，若存在x，对于任意y满足xΔy=yΔx=y，则称x为此运算系统的单位元。

例如：0为实数加法的单位元，1为实数乘法单位元。

在矩阵幂运算中，单位元为从左上到右下为1，其余都为0的矩阵。

例如：若是一个3*3的矩阵求幂，则他的单位元为

快速幂 

当知道单位元后，就可以按照快速幂的格式将矩阵快速幂写出了。

void unit() {
        memset(a,0,sizeof a);
        for(int i=1;i<=n;i++)a[i][i]=1;
}
s pow(int b){
        s ans(n,m),a=*this;ans.unit();
        while(b){if(b&1)ans=ans*a;a=a*a;b>>=1;}
        return ans;
}

其中unit()为给矩阵赋单位矩阵，pow(b)为求矩阵的b次幂。

带取模的矩阵快速幂

快速幂的优势就是便于取模，因此在写矩阵快速幂时，我们也需要将矩阵取模。

写一个函数将矩阵每个元素分别取模即可。

s operator%(const int t)const{
        s ans=*this;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)ans.a[i][j]%=t;
        return ans;
}

之后就可以把取模放到快速幂的模板里了。

s pow(int b,int p){
        s ans(n,m),a=*this;ans.unit();
        while(b){if(b&1)ans=ans*a%p;a=a*a%p;b>>=1;}
        return ans%p;
}

pow(b,p)表示矩阵的b次幂取模p。

模板

现在已经可以写出矩阵运算的模板了。

struct s{
    int a[105][105],r,c;
    s(int _r=0,int _c=0){
        r=_r,c=_c;memset(a,0,sizeof a);
        if(c==0)c=r;
    }
    void unit(){
        memset(a,0,sizeof a);
        for(int i=1;i<=r;i++)a[i][i]=1;
    }
    s operator*(const s t)const{
        s ans(r,t.c);int n=r,m=t.c;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
                for(int k=1;k<=c;k++)ans.a[i][j]+=a[i][k]*t.a[k][j];return ans;
    }
    s operator+(const s t)const{
        s ans(r,c);
        for(int i=1;i<=r;i++)
            for(int j=1;j<=c;j++)ans.a[i][j]=a[i][j]+t.a[i][j];
        return ans;
    }
    void print() {
        for(int i=1;i<=r;i++){
            for(int j=1;j<=c;j++)cout<<a[i][j]<<" ";cout << endl;
    }
    s pow(ll b) {
        s ans(r,c),a=*this;ans.unit();
        while(b){if(b&1)ans=ans*a;a=a*a;b>>=1;}
        return ans;
    }
    s operator%(const ll t)const{
        s ans=*this;
        for(int i=1;i<=r;i++)
            for (int j=1;j<=c;j++)ans.a[i][j]%=t;return ans;
    }
    s pow(ll b,ll p){
        s ans(r,c),a=*this;ans.unit();
        while(b){if(b&1)ans=ans*a%p;a=a*a%p;b>>=1;}
        return ans%p;
    }
};

题

斐波那契数列

题目描述

求斐波那契数列第n项Fn mod 1e9+7。

思路

若一个矩阵A乘上可以得出,那么只需要算出就可以得出

。

根据斐波那契数列定义可以推导出，A矩阵为。

AC代码