矩阵和线性方程组学习笔记

矩阵方法还是要会的

矩阵

矩阵的基本运算

矩阵加法
最简单了，就是对应数的相加.
如果用i,j来表示行和列
那么就有

$C_{(i,j)} = A_{(i,j)} + B_{(i,j)}$

减法同理，不作太多的阐述。唯一需要注意事项，当两个矩阵的行和列不相等的时候，求他们的矩阵和或者矩阵差是毫无意义的。只有r1=r2,s1=s2才有加减法的定义。

矩阵乘法（划重点）

稍显复杂，矩阵乘法必须在n行m列和m行p列的两个矩阵中进行运算，否则一样没有意义。

得到一个n行p列的矩阵。

有以下的递推式:

$C_{i,k} = \sum_{j=1}^nA_{i,j} *B_{j,k}$

具体细节

Markdown

根据这个递推式，我们不难得到一个O(n^3)的算法

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int maxn = 510;
const int mod = 1e9+7;

struct Matrix
{
	int m[maxn][maxn];
	int r, s;
};

int N,M,P;
Matrix A,B,ans;
signed main(signed argc, char *argv[])
{
	scanf("%lld %lld %lld",&N,&P,&M);
	for (int i = 1; i <= N; i++)
		for (int j = 1; j <= P; j++)
		scanf("%lld",&A.m[i][j]);
	system("pause");
	for (int i = 1; i <= P; i++)
		for (int j = 1; j <= M; j++)
		scanf("%lld",&B.m[i][j]);
		
	
	
	for (int i = 1; i <= N; i++)
		for (int k = 1; k <= P; k++)
			for (int j = 1; j <= M; j++)
			{
				ans.m[i][j] = (ans.m[i][j] + mod) % mod;
				ans.m[i][j] += (A.m[i][k]%mod)* (B.m[k][j]%mod);
				ans.m[i][j] = (ans.m[i][j] + mod) % mod;
			}
			
	for (int i = 1; i <= N; i++)
	{
		for (int j = 1; j <= M; j++)
		printf("%lld ",(ans.m[i][j]%mod+mod)%mod);
		printf("\n");
	}
	return 0;
}

尽管联赛的要求就是O(n^3)

但是还可以更快

可以试想，当A(i,j)等于0时,那里一整块的结果都是0，不需要继续运算

for (int i = 1; i <= N; i++)
		for (int k = 1; k <= P; k++)
			if (A.m[i][k])
			for (int j = 1; j <= M; j++)
			{
				ans.m[i][j] = (ans.m[i][j] + mod) % mod;
				ans.m[i][j] += (A.m[i][k]%mod)* (B.m[k][j]%mod);
				ans.m[i][j] = (ans.m[i][j] + mod) % mod;
			}

还有神仙操作的可以优化到O(nm)的操作

贴下代码慢慢理解啦。

高斯费马树上开花

我们

俯身欣赏.

#include<bits/stdc++.h>
typedef long long i64;
const int P=1e9+7;
char rb[8000000],*rp=rb,ob[4000000],*op=ob;
int _(){
	int x=0,f=1;
	while(*rp<48)*rp++=='-'?f=-1:0;
	while(*rp>47)x=x*10+*rp++-48;
	return x*f;
}
void pr(int x){
	x+=(x>>31&P);
	int ss[15],sp=0;
	do ss[++sp]=x%10,x/=10;while(x);
	while(sp)*op++=ss[sp--]+48;
}
int n,p,m;
int a[507][507],b[507][507];
i64 c[507][507];
int main(){
	fread(rb,1,sizeof(rb),stdin);
	n=_(),p=_(),m=_();
	for(int i=1;i<=n;++i)
		for(int j=1;j<=p;++j)a[i][j]=_();
	for(int i=1;i<=p;++i)
		for(int j=1;j<=m;++j)b[i][j]=_();
	int m1=m;
	while(m1&3)++m1;
	for(int i=1;i<=n;++i){
		i64*z=c[i];
		for(int k=1;k<=p;++k){
			int x=a[i][k],*y=b[k];
			for(int j=4;j<=m1;j+=4){
				z[j-3]+=i64(x)*y[j-3];
				z[j-2]+=i64(x)*y[j-2];
				z[j-1]+=i64(x)*y[j-1];
				z[j  ]+=i64(x)*y[j  ];
			}
			if(!(k&7))
			for(int j=1;j<=m;++j)z[j]-=(z[j]>>30)*P;
		}
	}
	for(int i=1;i<=n;++i){
		for(int j=1;j<=m;++j){
			pr(c[i][j]%P);
			*op++=32;
		}
		*op++=10;
	}
	fwrite(ob,1,op-ob,stdout);
	return 0;
}

高斯消元

我们假设有这样一个n阶的可逆矩阵，例如

$2x+y-z =8(L_1)\\\\-3x-y+2z=-11(L_2)\\\\-2x+y+2z=-3(L_3)$

这样的一个三元一次方程组写成矩阵的形式为

$\begin{bmatrix} 2 & 1 & -1 &| &8 \\\\ -3 & -1 & 2 & | & -11\\\\ -2 & 1 & 2 & | & -3 \end{bmatrix}$

这个矩阵被称作增广矩阵（AUGMENTED MATRIX）,最后一列事实上不是系数矩阵而是常数。

这个矩阵可以经过这样的处理。

对于每一行要使得Aii≠0且Aji(j>i)均为0,具体操作

1.来到当前行以上方矩阵为例子，现在要处理第一行。我们先往下面找到绝对值最大的A(i,r),这里找到是-3,于是交换第1行和第2行.（矩阵变成这个样子，就是交换了一二两行）

$\begin{bmatrix} -3 & -1 & 2 & | & -11\\\\ 2 & 1 & -1 &| &8 \\\\ -2 & 1 & 2 & | & -3 \end{bmatrix}$

对于L1，L2 将L1 3 +L2 2来加减消元，对于第三行，和第二行一样，直接与现在的L1操作。

将L2、L3的首项消成0后继续操作第二行。

$\begin{bmatrix} -3 & -1 & 2 & | & -11\\\\ \space & \frac{1}{3} & \frac{1}{3} & | &\frac{2}{3} \\\\ \space & \frac{5}{3} & \frac{2}{3} & | &\frac{13}{3} \\\\ \end{bmatrix}$

以此类推最终得到一个三角矩阵。

差不多就长这个样子。（滑稽

$\begin{bmatrix} -3 & -1 & 2 & | & -11\\\\ \space & \frac{5}{3} & \frac{2}{3} & | &\frac{13}{3} \\\\ \space & \space & \frac{1}{5} & | &-\frac{1}{5} \\\\ \end{bmatrix}$

这样其实最后一列已经告诉了我们z=-1。

自下向上迭代即可。

高斯消元代码如下。

void Gauss_Eli(){
    int num;
    for (int i=1;i<n;i++){
        num=i;
        for (int j=i+1;j<=n;j++)
          if (fabs(A[j][i])>fabs(A[num][i]))
            num=j;
        for (int j=1;j<=n;j++) swap(A[i][j],A[num][j]);
        swap(b[num],b[i]);
        for (int j=i+1;j<=n;j++){
            if (fabs(A[j][i])<=eps) continue;
            double t=A[j][i]/A[i][i];
            for (int k=1;k<=n;k++)
              A[j][k]-=A[i][k]*t;
            b[j]-=b[i]*t;
        }
    }
    for (int i=n;i>=1;i--){
        ans[i]=b[i]/A[i][i];
        for (int j=i;j>=1;j--)
          b[j]-=A[j][i]*ans[i];
    }
}

对于异或^方程组可以用bitset优化.

bitset<250>A[250]; \\常数项b[i]贴在A[i][n+1]中.
void Gauss_Eli(){
    int num;
    for (int i=1;i<=n;i++){
        if (A[i][i]==0){
            num=i;
            for (int j=i+1;j<=n;j++)
              if (A[j][i]!=0){num=j;break;}
            swap(A[num],A[i]);
        }
        for (int j=i+1;j<=n;j++)
          if (A[j][i]!=0) A[j]^=A[i];
    }
    for (int i=n;i>=1;i--){
        ans[i]=A[i][n+1];
        for (int j=i;j>=1;j--)
          if (A[j][i]!=0)
            A[j][n+1]=A[j][n+1]^ans[i];
    }
}

题目传送门:

UVa-10870
LA-3704
UVa-11542
POJ-1222