矩阵

定义矩阵乘法：$C=AB$

如果$A$是$m*n$矩阵，$B$是$n*p$矩阵，那么$C$是一个$m*p$矩阵。

其中$C_{ik}=\sum_{k=1}^{n}A_{ij}B_{jk}$，如：

$$\begin{bmatrix} 1 & 2 & 3\\ 4 & 5 &6\\ \end{bmatrix} \begin{bmatrix} 1 & 2\\ 3 & 4\\ 5 & 6\\ \end{bmatrix}= \begin{bmatrix} 1+6+15 & 2+8+18\\ 4+15+30 & 8+20+36\\ \end{bmatrix}= \begin{bmatrix} 22 & 28\\ 49 & 64\\ \end{bmatrix}$$
矩阵乘法满足结合律，不满足交换律。

单位矩阵$In$：

$$In= \begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 1\\ \end{bmatrix}$$

矩阵乘法及快速幂：

struct Matrix{
    int m[M][M];
    Matrix(){
        memset(m,0,sizeof(m));
    }
    void Init(){//初始化单位矩阵 
        for(int i=1;i<=n;i++)
            m[i][i]=1;
    }
    Matrix operator *(const Matrix &a)const{
        Matrix res;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++){
                if(!m[i][j])continue;
                for(int k=1;k<=n;k++)
                    res.m[i][k]+=m[i][j]*a.m[j][k];
            }
        return res;
    }
};
Matrix Pow(Matrix x,int t){
    Matrix res;
    res.Init();
    while(t){
        if(t&1)res=res*x;
        x=x*x;
        t>>=1;
    }
    return res;
}

POJ3070 Fibonacci

Task：

求斐波那契数列第n位取模P的值。（P=10000，n<=1000000000）

Solution：
首先直接循环是O(n)的。

考虑矩阵快速幂：

$$\begin{bmatrix} f[n]\\ f[n-1]\\ \end{bmatrix}= \begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix} \begin{bmatrix} f[n-1]\\ f[n-2]\\ \end{bmatrix}$$

$$\begin{bmatrix} f[n]\\ f[n-1]\\ \end{bmatrix}= \begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix}^{n-1} \begin{bmatrix} f[1]\\ f[0]\\ \end{bmatrix}= \begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix}^{n-1} \begin{bmatrix} 1\\ 0\\ \end{bmatrix} 设为A^{n-1}B$$

于是将这个矩阵$A$快速幂，就可以在$O(logn)$时间内解决这道题了。

#include<stdio.h>
#include<string.h>
#define M 5
#define P 10000
const int n=2;
struct Matrix{
    int m[M][M];
    Matrix(){
        memset(m,0,sizeof(m));
    }
    void Init(){//初始化单位矩阵 
        for(int i=1;i<=n;i++)
            m[i][i]=1;
    }
    Matrix operator *(const Matrix &a)const{
        Matrix res;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++){
                if(!m[i][j])continue;
                for(int k=1;k<=n;k++)
                    res.m[i][k]=(res.m[i][k]+m[i][j]*a.m[j][k])%P;
            }
        return res;
    }
};
Matrix Pow(Matrix x,int t){
    Matrix res;
    res.Init();
    while(t){
        if(t&1)res=res*x;
        x=x*x;
        t>>=1;
    }
    return res;
}
int main(){
    Matrix B;
    B.m[1][1]=1;B.m[2][1]=0;
    while(1){
        int d;
        scanf("%d",&d);
        if(d==-1)break;
        if(!d){
            puts("0");
            continue;
        }
        Matrix A;
        A.m[1][1]=A.m[1][2]=A.m[2][1]=1;A.m[2][2]=0;
        A=Pow(A,d-1);
        A=A*B;
        printf("%d\n",A.m[1][1]);
    }
    return 0;
}