矩阵问题（模板）

define MAXN 100

#define fabs(x) ((x)>0?(x):-(x))
#define zero(x) (fabs(x)<1e-10)

struct mat{
int n,m;
double data[MAXN][MAXN];
};

int mul(mat& c,const mat& a,const mat& b){
int i,j,k;
if (a.m!=b.n)
return 0;
c.n=a.n,c.m=b.m;
for (i=0;i<c.n;i++)
for (j=0;j<c.m;j++)
for (c.data[i][j]=k=0;k<a.m;k++)
c.data[i][j]+=a.data[i][k]*b.data[k][j];
return 1;
}

int inv(mat& a){
int i,j,k,is[MAXN],js[MAXN];
double t;
if (a.n!=a.m)
return 0;
for (k=0;k<a.n;k++){
for (t=0,i=k;i<a.n;i++)
for (j=k;j<a.n;j++)
if (fabs(a.data[i][j])>t)
t=fabs(a.data[is[k]=i][js[k]=j]);
if (zero(t))
return 0;
if (is[k]!=k)
for (j=0;j<a.n;j++)
t=a.data[k][j],a.data[k][j]=a.data[is[k]][j],a.data[is[k]][j]=t;
if (js[k]!=k)
for (i=0;i<a.n;i++)
t=a.data[i][k],a.data[i][k]=a.data[i][js[k]],a.data[i][js[k]]=t;
a.data[k][k]=1/a.data[k][k];
for (j=0;j<a.n;j++)
if (j!=k)
a.data[k][j]*=a.data[k][k];
for (i=0;i<a.n;i++)
if (i!=k)
for (j=0;j<a.n;j++)
if (j!=k)
a.data[i][j]-=a.data[i][k]*a.data[k][j];
for (i=0;i<a.n;i++)
if (i!=k)
a.data[i][k]*=-a.data[k][k];
}
for (k=a.n-1;k>=0;k--){
for (j=0;j<a.n;j++)
if (js[k]!=k)
t=a.data[k][j],a.data[k][j]=a.data[js[k]][j],a.data[js[k]][j]=t;
for (i=0;i<a.n;i++)
if (is[k]!=k)
t=a.data[i][k],a.data[i][k]=a.data[i][is[k]],a.data[i][is[k]]=t;
}
return 1;
}

double det(const mat& a){
int i,j,k,sign=0;
double b[MAXN][MAXN],ret=1,t;
if (a.n!=a.m)
return 0;
for (i=0;i<a.n;i++)
for (j=0;j<a.m;j++)
b[i][j]=a.data[i][j];
for (i=0;i<a.n;i++){
if (zero(b[i][i])){
for (j=i+1;j<a.n;j++)
if (!zero(b[j][i]))
break;
if (j==a.n)
return 0;
for (k=i;k<a.n;k++)
t=b[i][k],b[i][k]=b[j][k],b[j][k]=t;
sign++;
}
ret*=b[i][i];
for (k=i+1;k<a.n;k++)
b[i][k]/=b[i][i];
for (j=i+1;j<a.n;j++)
for (k=i+1;k<a.n;k++)
b[j][k]-=b[j][i]*b[i][k];
}
if (sign&1)
ret=-ret;
return ret;
}