线性方程组求解器 [Linear Equation Solver]
求解线性方程Ax=B，其中A是N*N的矩阵，B是有N个元素的列向量。若可解则输出唯一解，或输出”multiple”，无解时输出”inconsistent”。
题目链接：
https://open.kattis.com/problems/equationsolver
如果必然有解，用行列式法很容易。行列式的定义就符合递归，非常简洁地就可以实现对行列式的求解.例如用Python用二维列表a表示行列式，则求行列式的值det(a)就可以设计成：
def det(a):
if len(a)==1:
return a[0][0]
res = 0
sign = 1
for i in range(len(a)):
subMat = []
for j in range(1,len(a)):
subLine = []
for k in range(len(a)):
if k!=i:
subLine.append(a[j][k])
subMat.append(subLine)
res += sign*a[0][i]*det(subMat)
sign = -sign
return res但是此题可能出现行列式为0的情况，当行列式为0时，可能无解，可能有多解，这就需要处理矩阵A，在矩阵变换过程中，若出现某一行全为0，若此时这一行对应的增广矩阵的对应元素不是0，则必然无解，对应元素是0则是多解（如果有多行都是全0，只要有一行对应增广矩阵不是0就是无解的）
几乎是无可避免的要处理矩阵A，按照如下方法：将A通过交换与行变换形成上三角矩阵（在此过程中可以判断是否是多解或无解，只有唯一解才能化成三角矩阵），然后将矩阵单位化。
main()函数可以写为：
#include <cstdio>
#include <cmath>
#define MAXN 107
double Mat[MAXN*MAXN];
double BMat[MAXN];
int main(void)
{
int N;
char zeroColumnFlag;
while(1)
{
scanf("%d",&N);
if(N==0)
break;
readMat(N);
zeroColumnFlag = solveMat(N);
if(zeroColumnFlag)
{
InvalidStatus(N);
}
else
{
identityMatrix(N);
outputResult(N);
}
}
return 0;
}然后逐步把各函数按照数学方法设计完成即可。
void readMat(int N)
{
int i;
for(i=0;i<N*N;i++)
{
scanf("%lf",Mat+i);
}
for(i=0;i<N;i++)
{
scanf("%lf",BMat+i);
}
}
void outputMat(int N)
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
printf("%9.3f",Mat[i*N+j]);
}
printf("\t%9.3f\n",BMat[i]);
}
printf("\n");
}
int isZero(double x)
{
const double epsilon = 0.00000001;
if(fabs(x)<epsilon)
return 1;
else
return 0;
}
void swapElem(double *a,int i,int j)
{
double tmp = a[i];
a[i]=a[j];
a[j]=tmp;
}
void swapLine(int N,int k1,int k2)
{
int i;
for(i=0;i<N;i++)
{
swapElem(Mat,k1*N+i,k2*N+i);
}
}
void upperTriangularMat(int i,int N)
{
int k1,k2;
for (k1=i+1;k1<N;k1++)
{
double rate = Mat[k1*N+i]/Mat[i*N+i];
BMat[k1] -= rate *BMat[i];
for(k2=i;k2<N;k2++)
{
Mat[k1*N+k2] -= rate*Mat[i*N+k2];
}
}
}
void identityMatrix(int N)
{
int i,j;
for ( i=N-1;i>=0;i--)
{
double rate = 1/Mat[i*N+i];
BMat[i] = rate*BMat[i];
Mat[i*N+i] = 1;
for (j=0;j<i;j++)
{
BMat[j] -= BMat[i]*Mat[j*N+i];
Mat[j*N+i]=0;
}
}
}
void InvalidStatus(int N)
{
int i,j;
double s,*p;
for (i=0;i<N;i++)
{
s=0;
p = Mat+i*N;
for(j=0;j<N;j++)
{
s+=p[j]*p[j];
}
if(isZero(s) && !isZero(BMat[i]))
{
printf("inconsistent\n");
return ;
}
}
printf("multiple\n");
return ;
}
int solveMat(int N)
{
int i,k1;
int zeroColumnFlag=0;
for(i=0;i<N;i++)
{
zeroColumnFlag = 1;
for (k1 = i;k1<N;k1++)
{
if(!isZero(Mat[k1*N+i]))
{
swapLine(N,k1,i);
swapElem(BMat,k1,i);
zeroColumnFlag = 0;
break;
}
}
if(zeroColumnFlag)
{
break;
}
else
{
upperTriangularMat(i,N);
}
}
return zeroColumnFlag;
}
void outputResult(int N)
{
for (int i=0;i<N;i++)
{
if(i>0)
printf(" ");
printf("%.5lf",BMat[i]);
}
printf("\n");
}