最简行阶梯矩阵生成器

最简行阶梯矩阵生成器

描述：给一个矩阵$A$，输出其最简行阶梯矩阵$R=rref（A）$.
输入：第一行两个正整数$N$和$M$，表示矩阵的行数和列数。
第$i（2<=i<=N+1）$行$M$个整数，表示$A$的第$i-1$行的$M$个数。整数在 $int$ 范围内。
输出：$R=rref（A）$，其中整数输出本身，若为分数则输出$u/v$，其中$u、v$分别为该项分子分母。
C++实现：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define MAXN 1005
#define re register int
struct shu{ll u,v;}a[MAXN][MAXN];
int n,m,ppos[MAXN];
ll gcd(ll x,ll y){x=abs(x),y=abs(y);return !y?x:gcd(y,x%y);}
void yuefen(ll &u,ll &v){ll t=gcd(u,v);u/=t;v/=t;if(v<0)u=-u,v=-v;}
void add(ll &x,ll &y,ll s,ll t){// x/y + = s/t
	ll tx=x,ty=y;x=t*tx+s*ty;y=t*ty;yuefen(x,y);}
void jian(ll &x,ll &y,ll s,ll t){add(x,y,-s,t);}
void times(ll &x,ll &y,ll s,ll t){// x/y * = s/t
    x*=s;y*=t;yuefen(x,y);}
void devide(ll &x,ll &y,ll s,ll t){times(x,y,t,s);}
int findp(int x){for(re i=1;i<=m;i++)if(a[x][i].u)return i;return m+1;}
void guiyi(int x){int t=findp(x);if(a[x][t].u==1&&a[x][t].v==1)return;for(re i=t+1;i<=m;i++)devide(a[x][i].u,a[x][i].v,a[x][t].u,a[x][t].v);a[x][t].u=a[x][t].v=1;}
void eij(int x,int y,ll k1,ll k2){//row x - = row y * k1/k2
	for(re i=1;i<=m;i++)jian(a[x][i].u,a[x][i].v,a[y][i].u*k1,a[y][i].v*k2);}
void swaprow(int x,int y){for(re i=1;i<=m;i++)swap(a[x][i].u,a[y][i].u),swap(a[x][i].v,a[y][i].v);}
void init(){ll x;for(re i=1;i<=n;i++)for(re j=1;j<=m;j++)scanf("%lld",&a[i][j].u),a[i][j].v=1;}
int main(){
	freopen("rref.in","r",stdin);freopen("rref.out","w",stdout);
	scanf("%d%d",&n,&m);
	init();
	for(re i=1;i<=n;i++){ppos[i]=findp(i);guiyi(i);}
	for(re i=1;i<n;i++)for(re j=i+1;j<=n;j++)if(ppos[i]>ppos[j])swaprow(i,j),swap(ppos[i],ppos[j]);
	for(re i=1;i<=n;i++){
		ppos[i]=findp(i);
		guiyi(i);
		for(re j=i+1;j<=n;j++){
		    guiyi(j);if(!a[j][ppos[i]].u)break;
		    eij(j,i,a[j][ppos[i]].u,a[j][ppos[i]].v);
		}
		for(re j=i+1;j<=n;j++)ppos[j]=findp(j);
		for(re j=i+1;j<n;j++)for(re k=j+1;k<=n;k++)if(ppos[j]>ppos[k])swaprow(j,k),swap(ppos[j],ppos[k]);
	}
	for(re i=n;i>1;i--){
		if(ppos[i]==m+1)continue;
		for(re j=i-1;j>=1;j--)eij(j,i,a[j][ppos[i]].u,a[j][ppos[i]].v);
	}
	for(re i=1;i<=n;i++){
		for(re j=1;j<=m;j++){
			if(a[i][j].v==1){printf("%lld ",a[i][j].u);continue;}
			printf("%lld/%lld ",a[i][j].u,a[i][j].v);
		}
		printf("\n");
	}
}