线性代数#0 板子

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define lf double
const ll MAXN=105;
struct matrix{
	lf a[MAXN][MAXN];ll n,m;
	matrix(){
		memset(a,0,sizeof a);
	}
	inline void build(){   
		for(ll i=1;i<=n;++i)a[i][i]=1;
	}
}a;
matrix operator *(const matrix &x,const matrix &y){ 
	matrix z;
    z.n=x.n;z.m=y.m;
	for(ll k=1;k<=x.m;++k)
		for(ll i=1;i<=x.n;++i)
			for(ll j=1;j<=y.m;++j)
				z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j]);
	return z;
}
matrix operator +(const matrix &x,const matrix &y){ 
	matrix z;
    z.n=x.n;z.m=x.m;
	for(ll i=1;i<=x.n;++i)
		for(ll j=1;j<=x.m;++j)
			z.a[i][j]=x.a[i][j]+y.a[i][j];
	return z;
}
matrix operator -(const matrix &x,const matrix &y){ 
	matrix z;
    z.n=x.n;z.m=x.m;
	for(ll i=1;i<=x.n;++i)
		for(ll j=1;j<=x.m;++j)
			z.a[i][j]=x.a[i][j]-y.a[i][j];
	return z;
}

matrix matrixquickpower(ll k,matrix a)
{
	matrix ans;
	ans.n=ans.m=a.n;
	ans.build();
	do{ 
		if(k&1)ans=ans*a;
		a=a*a;k>>=1;
	}while(k);
	return ans;
}
lf caldet(matrix a)
{
	lf res=1,w=1,b[MAXN][MAXN];
	for(ll i=1;i<=a.n;i++)
		for(ll j=1;j<=a.m;j++)
			b[i][j]=a.a[i][j];
	for(ll i=1;i<=a.n;i++)
	{ 
		for(ll j=i+1;j<=a.n;++j)
		{
    		while(b[i][i])
			{
     	    	lf div=b[j][i]/b[i][i];
        		for(ll k=i;k<=a.n;++k)
        		    b[j][k]=(b[j][k]-div*b[i][k]);
        		swap(b[i],b[j]);w=-w;
    		}
    		swap(b[i],b[j]);w=-w;
		}
	}
	for(ll i=1;i<=a.n;i++)res=b[i][i]*res;
	res=w*res;
	return res;
}


ll jem(matrix a,matrix &ans)//Jordanian elimination method
{
    ans=a;
    ll n=a.n,m=a.m;
    for(ll i=1;i<=n;++i)
	{
		ll mx=i;
		for( ll j=i+1;j<=n;++j)
			if(fabs(ans.a[j][i])>fabs(ans.a[mx][i]))
				mx=j;
		for( ll j=1;j<=m;++j)
			swap(ans.a[i][j],ans.a[mx][j]);
		if(!ans.a[i][i])
			return -1;
		for( ll j=1;j<=n;++j)
		{
			if(j!=i)
			{
				lf temp=ans.a[j][i]/ans.a[i][i];
				for( ll k=i+1;k<=m;++k)
					ans.a[j][k]-=ans.a[i][k]*temp;
			}
		}
	}
	for(ll i=1;i<=n;i++)
		for(ll j=n+1;j<=m;j++)
			ans.a[i][j]/=ans.a[i][i];
    return 1;
}

ll Matrix_inversion(matrix a,matrix &ans)
{
    matrix b;ll n=a.n;
    b.n=a.n;b.m=a.n*2;
    ans.n=ans.m=n;
    for(ll i=1;i<=n;i++)
        for(ll j=1;j<=n;j++)
            b.a[i][j]=a.a[i][j];
    for(ll i=1;i<=n;i++)b.a[i][n+i]=1;
    matrix tmp;
    ll flag=jem(b,tmp);
    if(flag==-1)return -1;
    for(ll i=1;i<=n;i++)
        for(ll j=1;j<=n;j++)
            ans.a[i][j]=tmp.a[i][j+n];
    return 1;
}