#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>

/**
 * Hilfsmakro für {@code STR(X)}
 */
#define _STR(X) #X
/**
 * Makro das den gegebenen Wert in einen String convertiert
 */
#define STR(X) _STR(X)

/**
 * Die maximal erlaubte Länge eines Dateipfads
 */
#define MAX_FILE_PATH_LENGTH FILENAME_MAX

/**
 * Die maximal erlaubte Größe einer Matrix
 */
#define MAX_MATRIX_SIZE 500
/**
 * Die maximale Anzahl von verketteten doubles in den Dateien
 */
#define MAX_MATRIX_IN_FILE_SIZE (MAX_MATRIX_SIZE + 2)

/**
 * Maximale Zahl der Iterationsschritte
 */
#define MAX_ITERATION_STEPS 100

/**
 * Die Struktur zur Speicherung von Matrizen. Speicherplatz für die Daten muss dynamisch als Pointer-Array allokiert werden.
 */
typedef struct {
	int n;
	double** data;
} Matrix;
/**
 * Erstellt eine neue Matrix-Struktur der Größe 0 mittels dynamischer Speicher-Allokation.
 * @return Ein Zeiger auf die neue Struktur.
 */
Matrix* createMatrix(void);
/**
 * Gibt den Speicherplatz einer dynamisch allokierten Matrix wieder frei.
 */
void freeMatrix(Matrix* matrix);
/**
 * Füllt eine Matrix mit {@code rows} Zeilen.
 */
void createMatrixRows(Matrix* matrix, int rows);
/**
 * Schreibt eine Matrix in die Konsole.
 */
void printMatrix(Matrix* matrix);

/**
 * Die Struktur zur Speicherung von Vektoren. Speicherplatz für die Daten muss dynamisch allokiert werden.
 */
typedef struct {
	int n;
	double* data;
} Vector;
/**
 * Erstellt eine neue Vektor-Struktur der Größe 0 mittels dynamischer Speicher-Allokation.
 * @see initVector
 * @return Ein Zeiger auf die neue Struktur
 */
Vector* createVector(void);
/**
 * Allokiert den Speicher für einen {@code size} großen Vektor.
 * @param vector ein Zeiger auf einen nicht initialisierten Vektor
 * @param size die Größe des Vektors
 */
void initVector(Vector* vector, int size);
/**
 * Gibt den Speicherplatz eines dynamisch allokierten Vektors wieder frei.
 * @param vector ein Zeiger auf den freizugebenden Vektor
 */
void freeVector(Vector* vector);
/**
 * Schreibt einen Vektor in die Konsole.
 * @param vector ein Zeiger auf den Vektor
 */
void printVector(Vector* vector);

/**
 * Die Auflistung aller erlaubten Berechnungs-Algorithmen.
 */
typedef enum {
	JACOBI = 0, GAUSS_SEIDEL = 1
} Method;

/**
 * Beseitigt nicht eingelesene Daten aus der Konsoleneingabe.
 */
void flushStdin(void);

/**
 * Hauptfunktion die eine Matrix und zwei Vektoren aus der gegebenen Datei lädt und in den durch Zeiger gegebenen Strukturen speichert.
 * @param filename Dateipfad zu der Datei. Unter Windows immer absolut, unter Linux auch relativ.
 * @param A Ein Zeiger auf eine leere Matrix
 * @param b Ein Zeiger auf einen leeren Vektor
 * @param x Ein Zeiger auf einen leeren Vektor
 * @return Gibt an ob die Daten erfolgreich geladen wurden
 */
bool load(const char* filename, Matrix* A, Vector* b, Vector* x);
/**
 * Hauptfunktion die das durch eine Matrix und einen Vektor gegebene LGS mit einem bestimmten Startwert löst.
 * @param mmethod Der Algorithmus, der zur Lösung verwendet werden soll.
 * @param A Die Koeffizienten der Unbekannten gegeben als Zeiger auf die Matrix
 * @param b Die konstanten Glieder der Gleichungen gegeben als Zeiger auf den Vektor
 * @param x Ein Zeiger auf einen Vektor mit Startwerten
 * @param e Die Genauigkeit, auf die die Lösung bestimmt werden soll
 * @return {@code NULL} bei Fehlern. Ansonsten ein Vektor-Array der Größe {@code MAX_ITERATION_STEPS + 1}. Das Array beginnt mit dem Startwert und endet vorzeitig mit einem Vektor der Größe {@code 0}.
 */
Vector* solve(Method method, Matrix* A, Vector* b, Vector* x, double e);

/**
 * Liest ein double-Array aus der gegebenen Datei.
 * @param file Ein Zeiger auf die Datei
 * @param matrixLine Ein double-Array der Länge {@code maxCols} oder größer
 * @param maxCols die maximal zu lesende Anzahl an doubles
 * @return Die Anzahl der gelesenen doubles. Negativ wenn ein Fehler aufgetreten ist
 */
int readMatrixLine(FILE* file, double* matrixLine, int maxCols);

int main(int argc, char* argv[]) {
	char filePath[MAX_FILE_PATH_LENGTH + 1]; // char-Array zum Speichern des Dateipfads (1 char mehr für '\0')

	if(argc >= 2) { // Wenn Argumente übergeben wurden
		strncpy(filePath, argv[1], MAX_FILE_PATH_LENGTH); // Das erste Argument wird (auf MAX_FILE_PATH_LENGTH gekürzt) in filePath kopiert
	} else {
		puts("Please enter the path of the file you'd like to open");
		int result = scanf("%" STR(MAX_FILE_PATH_LENGTH) "[^\n]", filePath); // Liest MAX_FILE_PATH_LENGTH Zeichen aus der Konsole bis zum ersten Zeilenumbruch
		if(result == EOF || result == 0 || filePath[0] == 0) { // Wenn Nichts gelesen werden konnte wird das Programm erfolgreich beendet
			return 0;
		}
		if(result != 1) { // Wenn ein Fehler beim Lesen aufgetreten ist wird das Programm mit Fehlermeldung beendet
			fputs("Couldn't read file path - stopping\n", stderr);
			return 1;
		}
		flushStdin(); // Die Konsoleneingabe wird bereinigt
	}

	// Speicherallokation für die Matrix und die Vektoren
	Matrix* matrix = createMatrix();
	Vector* b = createVector();
	Vector* x = createVector();

	int returnCode = 0; // Variable um den Rückgabewert an das Programmende zu liefern

	while(true) { // Endlosschleife, die durch Erfolg oder Fehler beendet wird
		if(load(filePath, matrix, b, x)) { // Versucht die Daten zu laden
			puts("Please enter the algorithm to use:\n\t0: Jacobi (default)\n\t1: Gauss-Seidel");
			int algorithm; // int-Variable um Eingabewerte zwischenzuspeichern

			while(true) { // Endlosschleife, die durch erfolgreiche Eingabe oder Programmfehler beendet wird
				algorithm = getchar(); // Nutzereingabe
				if(algorithm == EOF) { // Bei der Eingabe ist ein Fehler aufgetreten
					returnCode = 20;
					goto end;
				}
				if(algorithm == '\n' || algorithm == 0) { // Wenn der Nutzer nichts eingegeben hat
					puts("Defaulting to Jacobi");
					algorithm = JACOBI; // Standardalgoritmus: JACOBI
					break; // Verlassen der Schleife
				}
				flushStdin(); // Konsoleneingabe bereinigen
				if(algorithm == '0' || algorithm == '1') { // Wenn die Eingabe gültig ist
					algorithm -= '0'; // Algorithmusnummer wird bestimmt
					break; // Verlassen der Schleife
				}
				fputs("Please enter 0, 1 or leave empty to exit!\n", stderr); // Fehlermeldung
			}

			puts("Please enter the precision to use:");
			double e; // Variable für die Algorithmusgenauigkeit
			int scan; // Variable zum Speichern des scanf Rückgabewerts
			while(true) {
				scan = scanf("%lf", &e); // Eingabe eines doubles
				if(scan == EOF) { // Fehler bei der Eingabe
					returnCode = 20;
					goto end;
				}
				flushStdin(); // Konsoleneingabe bereinigen
				if(scan == 1) { // Bei erfolgreicher Eingabe die Schleife verlassen
					break;
				}
				fputs("Invalid input - please enter a valid floating point number!\n", stderr); // Fehlermeldung
			}

			Vector* result = solve(algorithm, matrix, b, x, e); // Ausführung des Algorithmus

			if(result == NULL) { // Beim Auftritt eines Fehlers
				puts("Given equation system doesn't converge!\nEnter 0 to repeat the program (default) or 1 to exit"); // Fehlermeldung und Auswahlmöglichkeiten
				algorithm = getchar(); // Benutzereingabe

				if(algorithm != '\n' && algorithm != 0) { // Wenn Eingabe nicht leer war
					flushStdin(); // Konsoleneingabe bereinigen
					if(algorithm == '1') { // Bei Auswahl der Beendigung Programm erfolgreich beenden
						goto end;
					}
				}
			} else {
				puts("Solution found!\nEnter 0 to just print the solution (default) or 1 to print the whole sequence of iteration steps!"); // Auswahlmöglichkeiten
				algorithm = getchar(); // Benutzereingabe
				if(algorithm != '\n' && algorithm != 0) // Wenn Eingabe nicht leer
					flushStdin(); // Konsoleneingabe bereinigen
				if(algorithm == '1') { // Alle Schritte ausgeben
					int i = 0; // Zählervariable
					do { // Solange wie der aktuelle Vektor nicht leer ist und das Array nicht zu Ende ist
						printVector(&result[i]); // Vektor ausgeben
						free(result[i].data); // Speicher freigeben
					} while(result[++i].n != 0 && i <= MAX_ITERATION_STEPS);
					free(result); // Schließlich den Speicher des Arrays freigeben
				} else {
					int i = 0; // Zählervariable (erstes Element wird durch kopfgesteurte Schleife nicht überprüft)
					while(result[++i].n != 0 && i < MAX_ITERATION_STEPS) { // Solange wie der aktuelle Vektor nicht leer ist und das Array nicht zu Ende ist
						free(result[i - 1].data); // Speicher des vorherigen Vektors 
					}
					printVector(&result[i-1]); // Ausgabe des letzten existierenden Vektors
					free(result[i-1].data); // Speicher der letzten existierenden Vektordaten wird freigegeben
					free(result); // Vektoren werden freigegeben
				}
				break; // Hauptschleife verlassen
			}
		} else { // Datei konnte nicht geladen werden
			fputs("Failed to load data from file.\n", stderr);
		}

		puts("Enter file path or leave empty to exit");
		int result = scanf("%" STR(MAX_FILE_PATH_LENGTH) "[^\n]", filePath); // Erneute Abfrage des Dateipfads
		if(result == EOF || result == 0 || filePath[0] == 0) { // Wenn Eingabe leer ist
			goto end;
		}
		if(result != 1) { // Wenn Fehler bei der Eingabe aufgetreten ist
			fputs("Couldn't read file path - stopping\n", stderr);
			returnCode = 1;
			goto end;
		}
		flushStdin(); // Konsoleneingabe bereinigen
	}

end:
	// Matrix- und Vektorspeicher freigeben
	freeMatrix(matrix);
	freeVector(b);
	freeVector(x);

	printf("\nReturning with code %i\n", returnCode); // Den Rückgabewert anzeigen
	return returnCode;
}

bool load(const char* filename, Matrix* matrix, Vector* b, Vector* x) {
	FILE* file = fopen(filename, "r"); // Die gegebene Datei zum Lesen öffnen
	if(file == NULL) { // Beim Auftritt eines Fehlers -> Laden ist fehlgeschlagen
		fprintf(stderr, "Failed to open file \"%s\"\n", filename);
		return false;
	} else {
		// Erste Zeile der Datei wird mit Maximalgröße gelesen
		double* firstLineBuffer = malloc(sizeof(double) * (MAX_MATRIX_IN_FILE_SIZE));
		int colsInFile = readMatrixLine(file, firstLineBuffer, MAX_MATRIX_IN_FILE_SIZE); // Zeile lesen und Anzahl der gelesenen Werte speichern

		if(colsInFile == -1) { // Es konnte kein double gelesen werden
			fputs("Unexpected input on line 0", stderr);
			free(firstLineBuffer); // Speicher freigeben
			goto failure; // Zum Fehlerlabel springen
		} else if(colsInFile == -2) { // Die Zeile war länger als die Maximalgröße erlaubt
			fprintf(stderr, "Exceeded maximum matrix size of %i\n", MAX_MATRIX_SIZE);
			free(firstLineBuffer);
			goto failure; // Zum Fehlerlabel springen
		} else { // Erfolgreich gelesen
			int cols = colsInFile - 1; // Tatsächliche Größe der Matrix/Vektoren ist um eins kleiner, da zumindest b mit eingelesen wurde
			initVector(b, cols); // Vektoren werden initialisiert
			initVector(x, cols);

			b->data[0] = firstLineBuffer[cols]; // b0 enthält den letzten Eintrag der ersten Zeile
			createMatrixRows(matrix, cols); // Matrix wird mit cols Zeilen angelegt (nur das Haupt-Array - die Einträge sind nicht initialisiert)
			matrix->data[0] = firstLineBuffer; // Erste Matrix-Zeile entspricht erster eingelesener Zeile. b-Wert ist irrelevant da durch Matrix-Größe die erste Zeile vorher endet

			int colsInLine; // Variable zum Speichern der pro Zeile eingelesenen Einträge
			for(int i = 1; i < cols; i++) { // Einlesen von cols Zeilen
				matrix->data[i] = malloc(sizeof(double) * (colsInFile)); // Matrix Zeile wird initialisiert
				colsInLine = readMatrixLine(file, matrix->data[i], colsInFile); // Nächste Zeile in Matrix einlesen
				if(colsInLine < 0) { // Wenn ein Fehler beim Einlesen aufetreten ist
					if(i == cols - 1) { // Wenn die vermeintlich letzte Zeile erreicht ist und ein Fehler auftritt, dann wird davon ausgegangen, dass x-Werte übergeben wurden
						puts("Optional parameters are being used");

						cols--; // Die tatsächliche Matrixgröße ist also eins geringer
						b->n--; // Ebenso für die Vektoren
						x->n--;
						matrix->n--;
						free(matrix->data[i]); // der bereits allokierte Speicher kann wieder freigegeben werden

						// Die letzte Spalte (derzeit b) ist eigentlich x, die Zeiger auf die Daten können also getauscht werden
						{
							double* temp = x->data;
							x->data = b->data;
							b->data = temp;
						}

						// b wird nun aus der zugehörigen Spalte in der Matrix ausgelesen
						for(int j = 0; j < cols; j++) {
							b->data[j] = matrix->data[j][cols];
						}
						goto success; // Ausführung erfolgreich beendet
					} else { // Unerwarteter Fehler ist aufgetreten
						fprintf(stderr, "Line %i contains illegal formatting - please fix!\n", i + 1); // Fehlermeldung
						goto failure; // Abbruch
					}
				} else if(colsInLine != colsInFile) { // Die Anzahl der gelesenen Spalten stimmt nicht mit der Erwartung überein
					fprintf(stderr, "Illegal line length found in line %i\n", i + 1); // Fehlermeldung
					goto failure; // Abbruch
				} else { // Erfolgreiches Einlesen
					b->data[i] = matrix->data[i][cols]; // b-Wert auslesen
				}
			} // Ende Schleife

			// Da an dieser Stelle keine x-Werte gelesen wurden, wird x mit Nullen initialisiert
			memset(x->data, 0, cols * sizeof(double));
		}
	}
success: // Erfolgslabel
	fclose(file); // Die Datei wird geschlossen
	return true;

failure: // Fehlerlabel
	fclose(file); // Die Datei wird geschlossen
	return false;
}


int readMatrixLine(FILE* file, double* matrixLine, int maxCols) {
	char nextChar = 0; // Variable für den nächsten zu verarbeitenden char
	int col = 0; // Die Anzahl gelesener Werte

	do { // Solange wie Kommas als Trennzeichen gelesen werden
		if(fscanf(file, "%lf", &matrixLine[col]) != 1) { // In das Array wird ein double-Wert eingelesen
			return -1; // Es konnte kein double gelesen werden
		}

		col++; // Gelesene Werte erhöhen
		if(col > maxCols) { // Wenn zu viele Werte gelesen wurden
			return -2; // Fehlerrückgabe
		}
		nextChar = fgetc(file); // nächstes Trennzeichen einlesen
	} while (nextChar == ',');
	return col; // Anzahl gelesener Werte zurückgeben
}

void flushStdin(void) {
	int c; // Buffer-Variable
	while((c = getchar()) != '\n' && c != EOF && c != 0); // Solange das Konsoleneingabe-Ende noch nicht erreicht ist werden einzelne chars eingelesen
}

Vector* solve(Method method, Matrix* A, Vector* b, Vector* x, double e) {
	Vector* vectors = malloc(sizeof(Vector) * (MAX_ITERATION_STEPS + 1)); // Historie der Iterationsschritte/Rückgabe-Array

	initVector(&vectors[0], b->n); // Erster Vektor wird initialisiert...
	memcpy(vectors[0].data, x->data, b->n * sizeof(double)); // ...und mit den Werten aus x gefüllt

	int vectorCount = 1; // Länge existierender Vektoren in vectors
    double temp; // Temporäre Variable zur Berechnung
	double delta; // Abweichung der Iterationsschritte - zur Genauigkeitsberechnung

	if (method == JACOBI){
		do { // Solange wie die Abweichung größer als die Genauigkeit ist
			delta = 0.0;
			initVector(&vectors[vectorCount], b->n); // Vektor an aktueller Position initialisieren

			for (int i = 0; i < b->n; i++) { // Iteration über alle Zeilen
				temp = 0.0; // temp wird zur Summenbildung verwendet
				for (int j = 0; j < b->n; j++) { // für alle Spalten, die nicht auf der Diagonale liegen
					if (j != i){
						temp = temp + A->data[i][j] * vectors[vectorCount - 1].data[j];
					}
				}
				vectors[vectorCount].data[i] = (b->data[i] - temp) / A->data[i][i]; // Berechneten Wert in Historie aufnehmen
				delta = fmax(fabs(vectors[vectorCount].data[i] - vectors[vectorCount - 1].data[i]), delta); // Berechnung der maximalen Abweichung durchführen
			}
			vectorCount++; // Größe der Historie erhöhen
			if (vectorCount >= MAX_ITERATION_STEPS) // Wenn zu viele Iterationsschritte durchgeführt wurden
				goto fail; // Fehlerbehandlung
		} while (delta > e);
	} else { // GAUSS-SEIDEL
		do { // Solange wie die Abweichung größer als die Genauigkeit ist
            delta = 0.0;
			initVector(&vectors[vectorCount], b->n); // Vektor an aktueller Position initialisieren
			memcpy(vectors[vectorCount].data, vectors[vectorCount - 1].data, b->n * sizeof(double)); // Werte des vorherigen Schrittes in den neuen Schritt kopieren

			for (int i = 0; i < b->n; i++) { // Iteration über alle Zeilen
				temp = 0.0; // temp wird zur Summenbildung verwendet
				for (int j = 0; j < b->n; j++) { // für alle Spalten, die nicht auf der Diagonale liegen
                    if (j != i){
                        temp = temp + A->data[i][j] * vectors[vectorCount].data[j];
                    }
                }
				vectors[vectorCount].data[i] = (b->data[i] - temp) / A->data[i][i]; // Berechneten Wert in Historie aufnehmen
				delta = fmax(fabs(vectors[vectorCount].data[i] - vectors[vectorCount - 1].data[i]), delta); // Berechnung der maximalen Abweichung durchführen
            }
			vectorCount++; // Größe der Historie erhöhen
			if (vectorCount >= MAX_ITERATION_STEPS) // Wenn zu viele Iterationsschritte durchgeführt wurden
				goto fail; // Fehlerbehandlung
        } while(delta > e);
	}

	vectors[vectorCount].n = 0; // Abschließenden Vektor auf Größe 0 setzen
	return vectors;

fail: // Fehlerlabel
	for (int i = 0; i < vectorCount; i++) { // Allokierten Daten-Speicher freigeben
		free(vectors[i].data);
	}
	
	free(vectors); // Vektoren-Array freigeben
	return NULL;
}

inline Matrix* createMatrix(void) {
	Matrix* matrix = malloc(sizeof(Matrix));
	matrix->n = 0;
	return matrix;
}

inline void createMatrixRows(Matrix* matrix, int rows) {
	matrix->n = rows;
	matrix->data = malloc(sizeof(double*) * rows);
}

inline void freeMatrix(Matrix* matrix) {
	for(int i = 0; i < matrix->n; i++) { // Gibt die Zeilen frei
		free(matrix->data[i]);
	}
	free(matrix->data); // Gibt das Zeilen-Array frei
	free(matrix); // Gibt das Matrix-Object frei
}

void printMatrix(Matrix* matrix) {
	for(int i = 0; i < matrix->n; i++) {
		for(int j = 0; j < matrix->n; j++) {
			printf("%lf, ", matrix->data[i][j]);
		}
		puts("");
	}
}

inline Vector* createVector(void) {
	Vector* vector = malloc(sizeof(Vector));
	vector->n = 0;
	return vector;
}

inline void initVector(Vector* vector, int size) {
	vector->n = size;
	vector->data = malloc(sizeof(double) * size);
}

inline void freeVector(Vector* vector) {
	free(vector->data);
	free(vector);
}

void printVector(Vector* vector) {
	for(int i = 0; i < vector->n; i++) {
		printf("%lf, ",  vector->data[i]);
	}
	puts("");
}