矩阵

luckydrawbox · 2021-08-15 12:45:41 · 算法·理论

矩阵结构体

变量

a[N][N]：表示矩阵。其中 N 为自定义常数。
h：行数。
l：列数。

函数

clean()：把结构体初始化。
build()：把结构体改成单位矩阵。若单位矩阵乘以某个 S，结果就是 S。
K_Recursion(k,*b)：构造 k 阶常系数线性递推项，b 是一个数组，构造的即是对于递推式 f_n=b_1f_{n-1}+b_2f_{n-2}+\cdots+b_kf_{n-k} 使 \begin{bmatrix}f_{n-1}&f_{n-2}&f_{n-3}&\cdots&f_{n-k}\end{bmatrix}\times A\Rrightarrow\begin{bmatrix}f_{n}&f_{n-1}&f_{n-2}&\cdots&f_{n-k+1}\end{bmatrix} 的矩阵 A，A 一般为 \begin{bmatrix}b_1&1&0&\cdots&0\\b_2&0&1&\cdots&0\\\vdots&\vdots&\vdots&\ddots&0\\b_{k-1}&0&0&\cdots&1\\b_k&0&0&\cdots&0\end{bmatrix}。

struct Juz{
    long long a[N][N],h,l;
    void clean(){
        memset(a,0,sizeof a);
        h=l=0;
    }
    long long & operator()(int x,int y){
        return a[x][y];
    }
    void build(){
        for(int i=1;i<N;++i)
            a[i][i]=1;
    }
    void K_Recursion(int k,int *b){
        clean();
        h=l=k;
        for(int i=1;i<=k;i++)
            a[i][1]=b[i];
        for(int i=1;i<k;i++)
            a[i][i+1]=1;
    }
};

矩阵加法

写法：a+b。

其中 a,b 都是 a.h\times a.l 的矩阵结构体，返回值也是 a.h\times a.l 的矩阵结构体。

Juz operator +(const Juz &x,const Juz &y){
    Juz z;
    z.clean();
    z.h=x.h;z.l=x.l;
    for(int i=1;i<=z.h;i++)
        for(int j=1;j<=z.l;j++)
            z.a[i][j]=x.a[i][j]+y.a[i][j];
    return z;
}

矩阵减法

写法：a-b。

其中 a,b 都是 a.h\times a.l 的矩阵结构体，返回值也是 a.h\times a.l 的矩阵结构体。

Juz operator -(const Juz &x,const Juz &y){
    Juz z;
    z.clean();
    z.h=x.h;z.l=x.l;
    for(int i=1;i<=z.h;i++)
        for(int j=1;j<=z.l;j++)
            z.a[i][j]=x.a[i][j]-y.a[i][j];
    return z;
}

矩阵乘法

写法：a*b。

其中 a,b 都是矩阵结构体，且 a.l=b.h，返回值是 a.h\times b.l 的矩阵结构体。

\begin{pmatrix}A&B&C\\D&E&F\end{pmatrix}\times\begin{pmatrix}a&b&c&d\\e&f&g&h\\i&j&k&l\end{pmatrix}=\begin{pmatrix}A\times a+B\times e+C\times i&A\times b+B\times f+C\times j&A\times c+B\times g+C\times k&A\times d+B\times h+C\times i\\D\times a+E\times e+F\times i&D\times b+E\times f+F\times j&D\times c+E\times g+F\times k&D\times d+E\times h+F\times i\end{pmatrix}

ans_{i,j}=\sum_{k=1}^{a.l}a_{i,k}\times b_{k,j}

答案矩阵的第 i 行第 j 列的数，就是由矩阵 a 第 i 行 a.l 个数与矩阵 b 第 j 列 a.l 个数分别相乘再相加得到的。

Juz operator *(const Juz &x,const Juz &y){
    Juz z;
    z.clean();
    z.h=x.h;z.l=y.l;
    for(int i=1;i<=x.h;++i)
        for(int j=1;j<=y.l;++j)
            for(int k=1;k<=x.l;++k)
                z.a[i][j]+=x.a[i][k]*y.a[k][j];
    //取模版
    //z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j]%mod)%mod; 
    return z;
}

矩阵快速幂

写法：Juzqmi(a,k)。

矩阵乘法符合结合律，按普通快速幂计算即可。

返回值是 a.h\times a.l 的矩阵结构体。

Juz Juzqmi(Juz a,long long k){ 
    Juz ans;
    ans.clean();
    ans.build();
    ans.h=a.h;ans.l=a.l;
    while(k){
        if(k&1)
            ans=ans*a;
        a=a*a;
        k>>=1; 
    }
    return ans;
}

矩阵求逆

将矩阵 x 的逆元存在 y 中，返回值为是否存在逆矩阵。

bool Juzinv(Juz &x,Juz &y){
    y.clean();
    y.h=x.h;y.l=x.l;
    y.build();
    for(int i=1;i<=x.h;i++){
        for(int j=i;j<=x.h;j++){
            if(x.a[j][i]){
                for(int k=1;k<=x.l;k++)
                    swap(x.a[i][k],x.a[j][k]);
                for(int k=1;k<=x.l;k++)
                    swap(y.a[i][k],y.a[j][k]);
                break;
            }
        }
        if(!x.a[i][i])return 0;
        ll rate=qmi(x.a[i][i],mod-2);
        for(int j=1;j<=x.l;j++){
            x.a[i][j]=x.a[i][j]*rate%mod;
            y.a[i][j]=y.a[i][j]*rate%mod;
        }
        for(int j=1;j<=x.h;j++){
            if(i==j||!x.a[j][i])continue;
            rate=x.a[j][i];
            for(int k=1;k<=x.l;k++){
                x.a[j][k]=(x.a[j][k]-rate*x.a[i][k]%mod+mod)%mod;
                y.a[j][k]=(y.a[j][k]-rate*y.a[i][k]%mod+mod)%mod;
            }
        }
    }
    return 1;
}

行列式求值

ll det(){
    ll f=1;
    for(int i=1;i<=h;i++){
        for(int j=i+1;j<=h;j++){
            while(a[i][i]){
                ll rate=a[j][i]/a[i][i];
                for(int k=1;k<=h;k++){
                    a[j][k]=(a[j][k]-rate*a[i][k]%mod+mod)%mod;
                    swap(a[j][k],a[i][k]);
                }
                f*=-1;
            }
            for(int k=1;k<=h;k++)swap(a[j][k],a[i][k]);f*=-1;
        }
    }
    f=(f+mod)%mod;
    for(int i=1;i<=h;i++)f=f*a[i][i]%mod;
    return f;
}

代码

//矩阵结构体
struct Juz{
    long long a[N][N],h,l;
    void clean(){
        memset(a,0,sizeof a);
        h=l=0;
    }
    void build(){
        for(int i=1;i<N;++i)
            a[i][i]=1;
    }
    void K_Recursion(int k,int *b){
        clean();
        h=l=k;
        for(int i=1;i<=k;i++)
            a[i][1]=b[i];
        for(int i=1;i<k;i++)
            a[i][i+1]=1;
    }
};
//矩阵加法
Juz operator +(const Juz &x,const Juz &y){
    Juz z;
    z.clean();
    z.h=x.h;z.l=x.l;
    for(int i=1;i<=z.h;i++)
        for(int j=1;j<=z.l;j++)
            z.a[i][j]=x.a[i][j]+y.a[i][j];
    return z;
}
//矩阵减法
Juz operator -(const Juz &x,const Juz &y){
    Juz z;
    z.clean();
    z.h=x.h;z.l=x.l;
    for(int i=1;i<=z.h;i++)
        for(int j=1;j<=z.l;j++)
            z.a[i][j]=x.a[i][j]-y.a[i][j];
    return z;
}
//矩阵乘法
Juz operator *(const Juz &x,const Juz &y){
    Juz z;
    z.clean();
    z.h=x.h;z.l=y.l;
    for(int i=1;i<=x.h;++i)
        for(int j=1;j<=y.l;++j)
            for(int k=1;k<=x.l;++k)
                z.a[i][j]+=x.a[i][k]*y.a[k][j];
    //取模版
    //z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j]%mod)%mod; 
    return z;
}
//矩阵快速幂，需要矩阵乘法
Juz Juzqmi(Juz a,long long k){ 
    Juz ans;
    ans.clean();
    ans.build();
    ans.h=a.h;ans.l=a.l;
    while(k){
        if(k&1)
            ans=ans*a;
        a=a*a;
        k>>=1; 
    }
    return ans;
}

back