// Gleichungslöser nach Gauss,
// mit dyn. angelegtem Array von Pointern auf Arrays
//
// Aufruf: gauss inputfile
//
// Klaus Kusche, 2013

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>

using namespace std;

void input(istream &infile, double **a, int varAnz);
bool solve(double **a, double *x, int varAnz);
void writeMatrix(double **a, int varAnz);
void writeResult(double *x, int varAnz);

// lies die Koeffizienten eines Gleichungssystems mit varAnz Variablen
// vom File infile in die Matrix a
void input(istream &infile, double **a, int varAnz)
{
  for (int i = 0; i < varAnz; ++i) {         // für alle Zeilen
    for (int j = 0; j <= varAnz; ++j) {       // ... und Spalten
      infile >> a[i][j];
    }
  }
}

// Berechne die Lösung x[j] aus der Koeffizientenmatrix a[i][j]
// Ergebnis true: Gelungen
// Ergebnis false: Keine oder keine eindeutige Lösung
bool solve(double **a, double *x, int varAnz)
{
  double t;

  // 1.) Gleichungssystem in Dreiecksform bringen
  // für alle Zeilen ...
  for (int i = 0; i < varAnz; ++i) {
    // i-ter Koeffizient in Zeile i darf nicht 0 sein,
    // sonst passiert im nächsten Schritt eine Division durch 0
    if (a[i][i] == 0) {
      for (int j = i + 1; ; ++j) {   // schau die Zeilen darunter an
        if (j == varAnz) {       // alle Zeilen durch, nichts gefunden
          return false;          // gib auf, nicht eindeutig lösbar
        }
        if (a[j][i] != 0) {      // Zeile j hat keine 0 in Spalte i
          // vertausche Zeile j und Zeile i elementweise, und zwar ab Spalte i
          // (Elemente vor Spalte i sind sowieso in beiden Zeilen 0)
          for (int k = i; k <= varAnz; ++k) { 
            t = a[i][k];
            a[i][k] = a[j][k];
            a[j][k] = t;
          }
          break;                 // fertig, restliche Zeilen nicht mehr prüfen
        }
      }
    }
    // eleminiere Spalte i in allen Zeilen unterhalb von Zeile i
    for (int j = i + 1; j < varAnz; ++j) {
      // Multiplikationsfaktor für Zeile i,
      // damit in Zeile j Spalte i (= a[j][i]) wegfällt
      t = a[j][i] / a[i][i];
      // Ziehe t mal Zeile i von Zeile j ab (elementweise)
      // für alle Spalten ab der i. Spalte, incl. der rechten Spalte
      for (int k = i; k <= varAnz; ++k) {    
        a[j][k] -= a[i][k] * t;
      }
    }
  }

  // 2.) Variablen von rechts unten nach links oben einzeln ausrechnen
  for (int i = varAnz - 1; i >= 0; --i) {
    // berechne x[i] (mit Koeffizient a[i][i]) 
    // bilde in t die Summe aller a[i][k] * x[k] in Zeile i,
    //   die *rechts* von x[i] stehen
    // alle a[i][k] links von x[i] sind 0, tragen also nichts bei!
    t = 0;
    for (int k = i + 1; k < varAnz; ++k)
      t += a[i][k] * x[k];
    // bring die Summe t nach rechts und rechne daraus x[i] aus
    x[i] = (a[i][varAnz] - t) / a[i][i];
  }

  return true;
}

// Gib das Gleichungssystem mit varAnz Variablen und den Koeffizienten a aus
void writeMatrix(double **a, int varAnz)
{
  for (int i = 0; i < varAnz; ++i) {
    for (int j = 0; j < varAnz; ++j) {
      if (j > 0) cout << " + ";
      cout << fixed << setw(10) << setprecision(4) << a[i][j] << "*x" << j;
    }
    cout << " = " << fixed << setw(10) << setprecision(4) << a[i][varAnz] << endl;
  }
  cout << endl;
}

// Gib die Lösungswerte für varAnz Variablen in x aus
void writeResult(double *x, int varAnz)
{
  for (int j = 0; j < varAnz; ++j) {
    if (j > 0) cout << ",  ";
    cout << fixed << setw(10) << setprecision(4) << x[j] << "=x" << j;
  }
  cout << endl;
}

int main(int argc, const char *argv[])
{
  // Mit wie vielen Variablen haben wir es wirklich zu tun?
  int varAnz;
  // Matrix mit dem Gleichungssystem = Eingabe:
  // Ist ein eindimensionales Array von varAnz vielen Pointern,
  // die auf varAnz viele eindimensionale Arrays mit je varAnz+1 doubles zeigen!
  // matrix[i][j] enthält den Koeffizienten der j. Variable in der i. Gleichung
  // matrix[i][varAnz] enthält das konstante Glied (rechte Seite) der i. Gleichung
  double **matrix;
  // Vektor der Variablenwerte = Ausgabe:
  // result[j] enthält die Lösung für die j. Variable
  double *result;

  if (argc != 2) {
    cerr << "Aufruf: " << argv[0] << " infilename" << endl;
    exit(EXIT_FAILURE);
  }
  ifstream inf(argv[1]);
  if (!inf) {
    cerr << argv[0] << ": Cannot open " << argv[1] << " for reading" << endl;
    exit(EXIT_FAILURE);
  }
  inf >> varAnz;
  if (varAnz <= 0) {
    cerr << "Die Anzahl der Variablen (" << varAnz <<
            ") muss größer 1 sein" << endl;
    exit(EXIT_FAILURE);
  }

  matrix = new double*[varAnz];
  for (int i = 0; i < varAnz; ++i) {
    matrix[i] = new double[varAnz + 1];
  }  
  result = new double[varAnz];

  input(inf, matrix, varAnz);
  writeMatrix(matrix, varAnz);
  if (solve(matrix, result, varAnz)) {
    writeMatrix(matrix, varAnz);
    writeResult(result, varAnz);
  } else {
    writeMatrix(matrix, varAnz);
    cout << "Keine / keine eindeutige Lösung" << endl;
  }

  for (int i = 0; i < varAnz; ++i) {
    delete [] matrix[i];
  }
  delete [] matrix;
  delete [] result;
  
  exit(EXIT_SUCCESS);
}