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;
}
Matrix problem (Template)