C++實現線性代數矩陣行簡化

#include<iostream>
#include<cmath>
#include<algorithm>
#include<iomanip>
#include<cstdio>
#include<string>
#define f(i,l,r) for(i=(l);i<=(r);i++)
#define ff(i,r,l) for(i=(r);i>=(l);i--)
#define ll long long
#define EPS 1e-6
using namespace std;
const int MAXN=105;
int n,m;
char output[MAXN];
string s;
struct frac{
  int x,y=1;
  bool operator < (const frac &tmp)const{
    return 1.0*x/y<1.0*tmp.x/tmp.y;
  }
  frac operator - (const frac &tmp){
    frac ans;
    ans.x=x*tmp.y-y*tmp.x;
    ans.y=y*tmp.y;
    ans.sim();
    return ans;
  }
  frac operator * (const frac &tmp){
    frac ans;
    ans.x=x*tmp.x;
    ans.y=y*tmp.y;
    ans.sim();
    return ans;
  }
  frac operator / (const frac &tmp){
    frac ans;
    ans.x=x*tmp.y;
    ans.y=y*tmp.x;
    ans.sim();
    return ans;
  }
  int gcd(int a,int b){
    return b==0?a:gcd(b,a%b);
  }
  void sim(){
    int d=gcd(x,y);
    x/=d;
    y/=d;
    if(x<0&&y<0){
      x=-x;
      y=-y;
    }
  }
  void write()
  {
    int i;
    char tmp[MAXN];
    int num=0;
    int nx=x,ny=y;
    s="";
    if(!nx){
      s+="0";
    }
    else{
      if(nx<0||ny<0){
        s+='-';
        nx=abs(nx);
        ny=abs(ny);
      }
      if(nx%ny==0){
        nx/=ny;
        while(nx){
          tmp[++num]='0'+nx%10;
          nx/=10;
        }
        ff(i,num,1){
          s+=tmp[i];
        }
      }
      else{
        while(nx){
          tmp[++num]='0'+nx%10;
          nx/=10;
        }
        ff(i,num,1){
          s+=tmp[i];
        }
        s+='/';
        num=0;
        while(ny){
          tmp[++num]='0'+ny%10;
          ny/=10;
        }
        ff(i,num,1){
          s+=tmp[i];
        }
      }
    }
    f(i,s.length(),7){
      cout<<" ";
    }
    cout<<s;
  }
}a[MAXN][MAXN];
void out_f()
{
  int i,j;
  f(i,1,n){
    f(j,1,m){
      cout<<1.0*a[i][j].x/a[i][j].y<<" ";
    }
    cout<<endl;
  }
}
void out()
{
  int i,j;
  f(i,1,n){
    f(j,1,m){
      a[i][j].write();
      cout<<" ";
    }
    cout<<endl;
  }
}
int find(int r,int c)
{
  int i,t=-1;
  f(i,r,n){
    if(t==-1||a[t][c]<a[i][c]) t=i;
  }
  return t;
}
void interchange(int r1,int r2)
{
  int j;
  f(j,1,m){
    swap(a[r1][j],a[r2][j]);
  }
  return;
}
void scale(int r,int c)
{
  int j;
  ff(j,m,c){
    a[r][j]=a[r][j]/a[r][c];
  }
  return;
}
void muilt(int r,int c)
{
  int i,j;
  f(i,r+1,n){
    if(!a[i][c].x) continue;
    ff(j,m,c){
      a[i][j]=a[i][j]-(a[i][c]*a[r][j]);
    }
  }
  return;
}
void gause_1()
{
  int c,r=1;
  f(c,1,n){
    int pos=find(r,c);
    if(!a[pos][c].x) continue;
    interchange(r,pos);
    scale(r,c);
    muilt(r,c);
    r++;
  }
}
void gause_2()
{
  int i,j;
  ff(i,n,1){
    f(j,i+1,n){
      a[i][n+1]=a[i][n+1]-(a[i][j]*a[j][n+1]);
      a[i][j].x=0;
    }
  }
  return;
}
int main()
{
  int i,j;
  cout<<"輸入矩陣的行數和列數："<<endl;
  cin>>n>>m;
  cout<<"輸入矩陣元素："<<endl;
  f(i,1,n){
    f(j,1,m){
      cin>>a[i][j].x;
    }
  }
  gause_1();
  cout<<"階梯型為："<<endl;
  out();
  gause_2();
  cout<<"最簡型為："<<endl;
  out();
  return 0;
}