/*
 * eliminacja gaussa v3
 * specjalnie dla drg
 * (c) 2006 by tear
 * MDJCK
 */
 

#include <stdlib.h>
#include <stdio.h>
#define DATASET5

static double * matrix;
static int cols;
static int rows;


static double m(int r, int c)
{
	return matrix[r * cols + c];
}

static double * pm(int r, int c)
{
	return matrix + r * cols + c;
}

static void dump()
{
	int i, j;
	
	printf("---\n");
	for (i = 0; i < rows; i++) {
		for (j = 0; j < cols; j++) {
			printf("%.2f ", m(i, j));
		}
		printf("\n");
	}	
}

static int find_row(int n)
{
	int i;

	for (i = n; i < rows; i++) {
		if (m(i, n))
			break;
	}
	if (i >= rows)
		i = -1;
	if (i != n) {
		printf("ha!\n");
	}
	return i;
}


static void swap_rows(int a, int b)
{
	double x;
	int i;
	
	/* disclaimer: verylame */
	for (i = 0; i < cols; i++) {
		x = m(a, i);
		*pm(a, i) = m(b, i);
		*pm(b, i) = x;
	}
}

int main()
{
	int i, j, k;
	double * xp;

#ifdef DATASET1
	/* jedno rozwiazanie */

	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ =  2; *xp++ =  1; *xp++ =  3; *xp++ =  5;
	*xp++ =  4; *xp++ =  3; *xp++ =  1; *xp++ =  6;
	*xp++ =  3; *xp++ =  4; *xp++ =  5; *xp++ =  2;
#endif
#ifdef DATASET2
	/* jedno rozwiazanie */
	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ =  0; *xp++ =  6; *xp++ =  2; *xp++ = 4;
	*xp++ = -3; *xp++ = -6; *xp++ = -3; *xp++ = 5;
	*xp++ =  6; *xp++ =  3; *xp++ =  4; *xp++ = -9;
#endif
#ifdef DATASET3
	/* nieskonczenie wiele rozwiazan, zaleznych od jednej zmiennej */
	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ =  0; *xp++ =  6; *xp++ =  2; *xp++ = 4;
	*xp++ = -3; *xp++ = -6; *xp++ = -3; *xp++ = 5;
	*xp++ =  0; *xp++ = 12; *xp++ =  4; *xp++ = 8;
#endif
#ifdef DATASET4
	/* brak rozwiazan */
	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ =  0; *xp++ =  0; *xp++ =  2; *xp++ = 4;
	*xp++ = -3; *xp++ =  0; *xp++ = -3; *xp++ = 5;
	*xp++ =  6; *xp++ =  0; *xp++ =  4; *xp++ = -9;
#endif
#ifdef DATASET5
	/* nieskonczenie wiele rozwiazan, zmienna z kolumny 2 (liczac od 1) - dowolna */
	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ =  0; *xp++ =  0; *xp++ =  2; *xp++ = 4;
	*xp++ = -3; *xp++ =  0; *xp++ = -3; *xp++ = 5.0/2.0;
	*xp++ =  6; *xp++ =  0; *xp++ =  4; *xp++ = -9;
#endif
	if (!matrix)
		abort();
	dump();

	for (i = 0; i < rows; i++) {
		double t1;
		int t4;

		t4 = find_row(i);
		if (t4 == -1) {
			continue; /* nothing to do */
		}

		if (t4 != i) {
			swap_rows(t4, i);
		}

		t1 = m(i, i);

		for (j = i; j < cols; j++) {
			*pm(i, j) /= t1;
		}
		dump();
#if 1
		for (j = i + 1; j < rows; j++) {
			double t2;

			t2 = m(j, i);
			for (k = i; k < cols; k++) {
				*pm(j, k) -= t2 * m(i, k);
			}
		}
#endif
	}
	printf("---\n");
	printf("mid\n");
	for (i = rows - 1; i >= 0; i--) {

		if (m(i,i) == 0) {
			for (j = i + 1; j < cols; j++) {
				if (m(i, j)) {
					break;
				}
			}
			if (j != cols)
				break;
		}
		
		for (j = i - 1; j >= 0; j--) {
			double t3;
			
			t3 = m(j, i);

#if 0			
			if (!t3) {
				abort();
			}
#endif
			for (k = i; k < cols; k++) {
				*pm(j, k) -= t3 * m(i, k);
			}
		}
		dump();
	}
	
	if (i != -1) {
		printf("no solutions\n");
	}
	
	return 0;
}