矩阵

矩阵结构体

变量

• 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