gauss.cpp

// g++ -o gauss gauss.cpp -DTYPE=Number

#include <cstdio>
#include <cstdlib>

// some helper functions and classes
static int gcd (int a, int b) {
	if (b == 0) return a;
	else return gcd (b, a%b);
}
static bool is_whitespace (char c) {
	return c == ' ' || c == '\t' || c == '\n';
}
static void getstr (char* str, int length) {
	int i = 0;
	char c = getchar ();
	while (c == ' ' || c == '\t' || c == '\n') {
		c = getchar ();
	}
	while (c != ' ' && c != '\t' && c != '\n' && i < length-1) {
		str[i++] = c;
		c = getchar ();
	}
	str[i] = '\0';
}
template <class T> class Array2D {
	int m, n;
	T* data;
public:
	Array2D (int m, int n): m(m), n(n) {
		data = (T*) malloc (m*n*sizeof(T));
	}
	Array2D (const Array2D& array): m(array.m), n(array.n) {
		data = (T*) malloc (m*n*sizeof(T));
		memcpy (data, array.data, m*n*sizeof(T));
	}
	~Array2D () {
		free (data);
	}
	Array2D& operator = (const Array2D& array) {
		if (m*n != array.m*array.n) {
			free (data);
			data = (T*) malloc (array.m*array.n*sizeof(T));
		}
		//memcpy
		m = array.m;
		n = array.n;
		return *this;
	}
	T* operator [] (int i) {
		return data + i * n;
	}
	const T* operator [] (int i) const {
		return data + i * n;
	}
};

// some fields
class Number {
	int numerator;
	int denominator;
public:
	static const Number ZERO;
	static const Number ONE;
	Number (int n): numerator(n), denominator(1) {
	
	}
	Number (int numerator, int denominator): numerator(numerator), denominator(denominator) {
		if (denominator < 0) {
			numerator = -numerator; // this
			denominator = -denominator;
		}
	}
	Number (const char* str): denominator(1) {
		numerator = atoi (str);
	}
	Number& cancel () {
		int gcd = ::gcd (abs(numerator), denominator);
		numerator /= gcd;
		denominator /= gcd;
		return *this;
	}
	Number operator - () const {
		return Number (-numerator, denominator);
	}
	Number operator + (Number n) const {
		if (denominator == n.denominator)
			return Number (numerator+n.numerator, denominator).cancel ();
		return Number (numerator*n.denominator+n.numerator*denominator, denominator*n.denominator).cancel ();
	}
	Number operator * (Number n) const {
		return Number (numerator*n.numerator, denominator*n.denominator).cancel ();
	}
	Number invert () const {
		return Number (denominator, numerator);
	}
	bool operator == (Number n) const {
		return numerator == n.numerator && denominator == n.denominator;
	}
	void print () const {
		if (denominator == 1)
			printf ("%d", numerator);
		else
			printf ("%d/%d", numerator, denominator);
	}
};
const Number Number::ZERO = 0;
const Number Number::ONE = 1;

template <int N> class Z {
	int value;
public:
	static const Z ZERO;
	static const Z ONE;
	Z (int value = 0): value(value%N) {
		
	}
	Z (const char* str) {
		value = atoi(str) % N;
	}
	Z operator - () const {
		return Z (N-value);
	}
	Z operator + (const Z& z) const {
		return Z (value + z.value);
	}
	Z operator * (const Z& z) const {
		return Z (value * z.value);
	}
	Z invert () const {
		for (int i=1; i<N; i++) {
			if ((value * i) % N == 1)
				return Z (i);
		}
	}
	bool operator == (Z z) const {
		return value == z.value;
	}
	void print () const {
		printf ("%d", value);
	}
};
template <int N> const Z<N> Z<N>::ZERO = 0;
template <int N> const Z<N> Z<N>::ONE = 1;

template <class T> class Matrix {
	int m, n;
	Array2D<T> data;
	int findNextRow (int j) {
		// find a row with value one
		for (int i = j; i < m; i++) {
			if (data[i][j] == T::ONE)
				return i;
		}
		// find a row with value != 0
		for (int i = j; i < m; i++) {
			if (!(data[i][j] == T::ZERO)) {
				multiply (i, data[i][j].invert());
				return i;
			}
		}
		return -1;
	}
public:
	Matrix (int m, int n): m(m), n(n), data(m,n) {
	
	}
	void scan () {
		for (int i=0; i<m; i++) {
			for (int j=0; j<n; j++) {
				char str[20];
				getstr (str, 20);
				data[i][j] = str;
			}
		}
	}
	void print () const {
		for (int i=0; i<m; i++) {
			data[i][0].print ();
			for (int j=1; j<n; j++) {
				putchar ('\t');
				data[i][j].print ();
			}
			putchar ('\n');
		}
	}
	void multiply (int i, const T& factor) {
		for (int j=0; j<n; j++) {
			data[i][j] = data[i][j] * factor;
		}
	}
	void add (int i1, int i2, const T& factor) {
		for (int j=0; j<n; j++) {
			data[i1][j] = data[i1][j] + data[i2][j] * factor;
		}
	}
	void swap (int i1, int i2) {
		for (int j=0; j<n; j++) {
			T tmp = data[i1][j];
			data[i1][j] = data[i2][j];
			data[i2][j] = tmp;
		}
	}
	void solve () {
		for (int j=0; j<n; j++) {
			int i = findNextRow (j);
			if (i == -1)
				break;
			if (i != j)
				swap (i, j);
			// normalize the other rows
			for (i = 0; i < j; i++) {
				add (i, j, -data[i][j]);
			}
			for (i = j+1; i < m; i++) {
				add (i, j, -data[i][j]);
			}
			print ();
			putchar ('\n');
		}
	}
};

int main (int argc, char** argv) {
	if (argc < 3) {
		printf ("Usage: %s [WIDTH] [HEIGHT]\n", argv[0]);
		return 0;
	}
	Complex c (0, -1);
	c.print();
	putchar ('\n');
	Matrix<TYPE> matrix (atoi(argv[1]), atoi(argv[2]));
	matrix.scan ();
	matrix.solve ();
	matrix.print ();
}