#pragma once
#include "../../stdafx.h"
#include "../glLibLocals.h"

#include "../glLibError.h"
#include "glLibMathFraction.h"

namespace glLib { namespace Math { namespace Matrix {

#define  MATRIX_TEMPLATE template <typename  type_element, typename  type_int>
#define MATRIX_TEMPLATE2 template <typename type_element2, typename type_int2>
#define  TYPE_MATRIX matrix< type_element, type_int>
#define TYPE_MATRIX2 matrix<type_element2,type_int2>

MATRIX_TEMPLATE class matrix;

#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
MATRIX_TEMPLATE void glLibInternal_check_3x3(const TYPE_MATRIX* matrix) {
	if ((matrix->_rows==3)&&(matrix->_columns==3)) return;
	throw GLLIB_ERROR_TYPE_ERROR("Cannot compute operation; matrix must be 3x3.");
}
MATRIX_TEMPLATE void glLibInternal_check_4x4(const TYPE_MATRIX* matrix) {
	if ((matrix->_rows==4)&&(matrix->_columns==4)) return;
	throw GLLIB_ERROR_TYPE_ERROR("Cannot compute operation; matrix must be 4x4.");
}
MATRIX_TEMPLATE void glLibInternal_check_3x3or4x4(const TYPE_MATRIX* matrix) {
	if ((matrix->_rows==3)&&(matrix->_columns==3)) return;
	if ((matrix->_rows==4)&&(matrix->_columns==4)) return;
	throw GLLIB_ERROR_TYPE_ERROR("Cannot compute operation; matrix must be 3x3 or 4x4.");
}
#endif

MATRIX_TEMPLATE void       add(const TYPE_MATRIX* m1, const TYPE_MATRIX* m2, TYPE_MATRIX* result) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (m1->_rows!=m2->_rows||m1->_columns!=m2->_columns) throw GLLIB_ERROR_TYPE_ERROR("Cannot add; dimensions do not match!");
	#endif
	for(type_int i=0; i<m1->_rows; ++i)
		for(type_int j=0; j<m1->_columns; ++j)
			result->values[i][j] = m1->values[i][j] + m2->values[i][j];
}
MATRIX_TEMPLATE void  subtract(const TYPE_MATRIX* m1, const TYPE_MATRIX* m2, TYPE_MATRIX* result) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (m1->_rows!=m2->_rows||m1->_columns!=m2->_columns) throw GLLIB_ERROR_TYPE_ERROR("Cannot add; dimensions do not match!");
	#endif
	for(type_int i=0; i<m1->_rows; ++i)
		for(type_int j=0; j<m1->_columns; ++j)
			result->values[i][j] = m1->values[i][j] - m2->values[i][j];
}
MATRIX_TEMPLATE void  multiply(const TYPE_MATRIX* m1, const TYPE_MATRIX* m2, TYPE_MATRIX* result) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	//http://mathworld.wolfram.com/MatrixMultiplication.html
	if (m1->_columns!=m2->_rows) {
		throw GLLIB_ERROR_TYPE_ERROR("Cannot multiply; dimensions are incompatible!");
	}
	if ((m1->_rows!=result->_rows)||(m2->_columns!=result->_columns)) {
		throw GLLIB_ERROR_TYPE_ERROR("Result matrix's dimensions are incorrect for the multiplication!");
	}
	#endif
	for (type_int i=0; i<m1->_rows; ++i) {
		for (type_int j=0; j<m2->_columns; ++j) {
			type_element sum = 0.0;
			for(type_int k=0; k<m1->_columns; ++k) { //m1->_columns should == m2->_rows
				sum += m1->values[i][k] * m2->values[k][j];
			}
			result->values[i][j] = sum;
		}
	}
}
MATRIX_TEMPLATE void transpose(TYPE_MATRIX* m) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (m->_rows!=m->_columns) throw GLLIB_ERROR_TYPE_ERROR("Cannot transpose a matrix that is not square!");
	#endif

	for (type_int i=0; i<m->_rows; ++i) {
		for (type_int j=i; j<m->_columns; ++j) {
			type_element temp = m->values[i][j];
			m->values[i][j] = m->values[j][i];
			m->values[j][i] = temp;
		}
	}
}
MATRIX_TEMPLATE void transpose(const TYPE_MATRIX* m, TYPE_MATRIX* result) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (m==result) throw GLLIB_ERROR_TYPE_ERROR("This overload does not operate in-place.");
	if (m->_rows!=result->_columns||m->_columns!=result->_rows) throw GLLIB_ERROR_TYPE_ERROR("Matrices' dimensions incompatible!");
	if (m->_rows!=m->_columns) throw GLLIB_ERROR_TYPE_ERROR("Cannot transpose a matrix that is not square!");
	#endif
	for (type_int i=0; i<m->_rows; ++i) {
		for (type_int j=0; j<m->_columns; ++j) {
			result->values[i][j] = m->values[j][i];
		}
	}
}
//TODO: determinant function, minor function
MATRIX_TEMPLATE inline type_element         det(const TYPE_MATRIX* matrix) { return determinant(matrix); }
MATRIX_TEMPLATE        type_element determinant(const TYPE_MATRIX* matrix) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (matrix->_rows!=matrix->_columns) throw GLLIB_ERROR_TYPE_ERROR("Cannot find the determinant of a non-square matrix!");
	#endif
	throw GLLIB_ERROR_TYPE_ERROR("Not implemented!");
}

/*static void minor(const TYPE_MATRIX* matrix, type_int column, type_int row, TYPE_MATRIX* result) {
	throw GLLIB_ERROR_TYPE_ERROR("Not implemented!");
}*/

// in matrix A, add x times row r2 to row r1
MATRIX_TEMPLATE void row_add_multiple(TYPE_MATRIX* m, type_element x, type_int r1, type_int r2) {
	for (type_int c=0; c<m->_columns; ++c) {
		m->values[r1][c] += x * m->values[r2][c];
	}
}
//TODO: swap column, column_add_multiple
MATRIX_TEMPLATE void row_swap(TYPE_MATRIX* m, type_int r1, type_int r2) {
	type_element* temp = m->values[r1];
	m->values[r1] = m->values[r2];
	m->values[r2] = temp;
}
MATRIX_TEMPLATE void row_scale(TYPE_MATRIX* m, type_element sc, type_int row) {
	for (type_int c=0; c<m->_columns; ++c) {
		m->values[row][c] *= sc;
	}
}
//TODO: ref
MATRIX_TEMPLATE void rref(TYPE_MATRIX* m) {
	type_int lead = 0;

	for (type_int r=0; r<m->_rows; ++r) {
		if (lead>=m->_columns) return;
		type_int i = r;
		while (m->values[i][lead]==0) {
			i++;
			if (i==m->_rows) {
				i = r;
				lead++;
				if (lead==m->_columns) {
					return;
				}
			}
		}
		row_swap(m, i,r);
		row_scale(m, (type_element)(1.0)/m->values[r][lead],r);
		for (i=0; i<m->_rows; ++i) {
			if (i!=r) {
				row_add_multiple(m, -m->values[i][lead], i, r);
			}
		}
		lead++;
	}
}

MATRIX_TEMPLATE void invert(const TYPE_MATRIX* m, TYPE_MATRIX* result) {
	throw GLLIB_ERROR_TYPE_ERROR("Not implemented!");
}

MATRIX_TEMPLATE void copy(const TYPE_MATRIX* src, TYPE_MATRIX* dst) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if ((src->_rows!=dst->_rows)||(src->_columns!=dst->_columns)) throw GLLIB_ERROR_TYPE_ERROR("Cannot copy; dimensions do not match!");
	#endif

	for(type_int i=0; i<src->_rows; ++i) {
		for(type_int j=0; j<src->_columns; ++j) {
			dst->values[i][j] = src->values[i][j];
		}
	}
}
MATRIX_TEMPLATE void copy(
	const TYPE_MATRIX* src, TYPE_MATRIX* dst,
	type_int src_left, type_int src_right, type_int src_top, type_int src_bottom,
	type_int dst_left, type_int dst_right, type_int dst_top, type_int dst_bottom)
{
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (  src_left<0) throw GLLIB_ERROR_TYPE_ERROR(  "Source's left index cannot be negative!");
	if ( src_right<0) throw GLLIB_ERROR_TYPE_ERROR( "Source's right index cannot be negative!");
	if (   src_top<0) throw GLLIB_ERROR_TYPE_ERROR(   "Source's top index cannot be negative!");
	if (src_bottom<0) throw GLLIB_ERROR_TYPE_ERROR("Source's bottom index cannot be negative!");
	if (  dst_left<0) throw GLLIB_ERROR_TYPE_ERROR(  "Destination's left index cannot be negative!");
	if ( dst_right<0) throw GLLIB_ERROR_TYPE_ERROR( "Destination's right index cannot be negative!");
	if (   dst_top<0) throw GLLIB_ERROR_TYPE_ERROR(   "Destination's top index cannot be negative!");
	if (dst_bottom<0) throw GLLIB_ERROR_TYPE_ERROR("Destination's bottom index cannot be negative!");

	//TODO
	/*if (left> right) throw GLLIB_ERROR_TYPE_ERROR("Left index cannot be greater than the right index!");
	if ( top>bottom) throw GLLIB_ERROR_TYPE_ERROR("Top index cannot be greater than the bottom index!");

	if ( right+1>src->_columns) throw GLLIB_ERROR_TYPE_ERROR("Requested index is outside the bounds of the source matrix!");
	if (bottom+1>src->   _rows) throw GLLIB_ERROR_TYPE_ERROR("Requested index is outside the bounds of the source matrix!");

	if ( right-left+1>dst->_columns) throw GLLIB_ERROR_TYPE_ERROR("Requested submatrix will not fit in destination matrix!");
	if (bottom- top+1>dst->   _rows) throw GLLIB_ERROR_TYPE_ERROR("Requested submatrix will not fit in destination matrix!");

	if ((bottom-top+1!=dst->_rows)||(right-left+1!=dst->_columns)) throw GLLIB_ERROR_TYPE_ERROR("Cannot copy; dimensions do not match!");
	*/
	#endif

	for (type_int i1=src_top,i2=dst_top; i1<=src_bottom; ++i1,++i2) {
		for (type_int j1=src_left,j2=dst_left; j1<=src_right; ++j1,++j2) {
			//printf("Setting (%d,%d) to (%d,%d)\n",i2,j2,i1,j1);
			dst->values[i2][j2] = src->values[i1][j1];
		}
	}
}

MATRIX_TEMPLATE void scale(type_element scalar,       TYPE_MATRIX* m) {
	for (type_int i=0; i<m->_rows; ++i)
		for (type_int j=0; j<m->_columns; ++j)
			m->values[i][j] = scalar * m->values[i][j];
}
MATRIX_TEMPLATE void scale(type_element scalar, const TYPE_MATRIX* m, TYPE_MATRIX* result) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	if (m->_rows!=result->_rows||m->_columns!=result->_columns) throw GLLIB_ERROR_TYPE_ERROR("Cannot scale; source and destination dimensions do not match!");
	#endif
	for (type_int i=0; i<m->_rows; ++i)
		for (type_int j=0; j<m->_columns; ++j)
			result->values[i][j] = scalar * m->values[i][j];
}

MATRIX_TEMPLATE        void set_rotation_matrix_rad(TYPE_MATRIX* rot_matrix, const type_element vec[3], type_element theta) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	glLibInternal_check_3x3or4x4(rot_matrix);
	#endif
	type_element c = cos(theta);
	type_element s = sin(theta);
	type_element C = 1.0 - c;
	type_element x = vec[0]; type_element y = vec[1]; type_element z = vec[2];
	type_element xs = x*s;   type_element ys = y*s;   type_element zs = z*s;
	type_element xC = x*C;   type_element yC = y*C;   type_element zC = z*C;
	type_element xyC = x*yC; type_element yzC = y*zC; type_element zxC = z*xC;
	type_element** rot_matrix_values = rot_matrix->values;
	rot_matrix_values[0][0] = x*xC+c; rot_matrix_values[0][1] = xyC-zs; rot_matrix_values[0][2] = zxC+ys;
	rot_matrix_values[1][0] = xyC+zs; rot_matrix_values[1][1] = y*yC+c; rot_matrix_values[1][2] = yzC-xs;
	rot_matrix_values[2][0] = zxC-ys; rot_matrix_values[2][1] = yzC+xs; rot_matrix_values[2][2] = z*zC+c;
}
MATRIX_TEMPLATE inline void set_rotation_matrix_deg(TYPE_MATRIX* rot_matrix, const type_element vec[3], type_element theta) {
	set_rotation_matrix_rad(rot_matrix,vec,radians(theta));
}

//Creates a matrix that transforms the unit i, j, and k vectors onto new vectors
MATRIX_TEMPLATE void set_transform(TYPE_MATRIX* matrix, const type_element new_i[3], const type_element new_j[3], const type_element new_k[3]) {
	//See http://www.gamedev.net/topic/604893-general-transformation-matrix/page__pid__4827470#entry4827470
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	glLibInternal_check_4x4(matrix);
	#endif
	type_element** values = matrix->values;
	values[0][0] = new_i[0]; values[0][1] = new_j[0]; values[0][2] = new_k[0]; values[0][3] = 0;
	values[1][0] = new_i[1]; values[1][1] = new_j[1]; values[1][2] = new_k[1]; values[1][3] = 0;
	values[2][0] = new_i[2]; values[2][1] = new_j[2]; values[2][2] = new_k[2]; values[2][3] = 0;
	values[3][0] =        0; values[3][1] =        0; values[3][2] =        0; values[3][3] = 1;
}

MATRIX_TEMPLATE void conjugate_gradient(const TYPE_MATRIX* A, const TYPE_MATRIX* x, const TYPE_MATRIX* b, int iterations) {
	#ifdef GLLIB_ERROR_MODE_REPORT_WARNINGS
	if (iterations>A->_rows) GLLIB_ERROR_TYPE_WARNING("Conjugate gradient should converge in no more than n iterations for A(n x n).");
	#endif

	TYPE_MATRIX A_times_x(b->_rows,1);
	TYPE_MATRIX A_times_p(b->_rows,1);

	TYPE_MATRIX r(b->_rows,1);
	TYPE_MATRIX r_next(b->_rows,1);

	TYPE_MATRIX p(b->_rows,1);
	TYPE_MATRIX alpha_p(b->_rows,1);
	TYPE_MATRIX beta_p(b->_rows,1);

	multiply(A,x,&A_times_x);

	subtract(b,&A_times_x,&r);
	copy(&r,&p);

	for (int k=0;k<iterations;++k) {
		multiply(A,&p, &A_times_p);
		type_element alpha = dot(&r,&r) / dot(&p,&A_times_p);

		scale(alpha,&p, &alpha_p);
		add(x,&alpha_p, x);

		multiply(A,&p, &A_times_p);
		scale(alpha,&A_times_p);
		subtract(&r,&A_times_p, &r_next);

		//if r_next is sufficiently small then break

		type_element beta = dot(&r_next,&r_next) / dot(&r,&r);

		scale(beta,&p, &beta_p);
		add(&r_next,&beta_p, &p);

		copy(&r_next,&r);
	}
}

MATRIX_TEMPLATE type_element dot(const TYPE_MATRIX* a, const TYPE_MATRIX* b) {
	/*printf("Attempting to dot a:\n");
	a->print();
	printf("With b:\n");
	b->print();*/
	type_element result = 0.0;
	if (a->_columns==1) {
		if (b->_columns==1) {
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			if (a->_rows!=b->_rows) goto ERROR_MATCH;
			#endif
			/*for (type_int i=0;i<a->_rows;++i) {
				type_element value1 = a->values[i][0];
				type_element value2 = b->values[i][0];
				type_element value = value1*value2;
				printf("Adding (a[%d,%d]=%f)*(b[%d,%d]=%f) = %f\n",0,i,value1,0,i,value2,value);
				result += value;
			}*/
			for (type_int i=0;i<a->_rows;++i) result += a->values[i][0]*b->values[i][0];
		}
		#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
		else if (b->_rows==1) {
		#else
		else {
		#endif
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			if (a->_rows!=b->_columns) goto ERROR_MATCH;
			#endif
			for (type_int i=0;i<a->_rows;++i) result += a->values[i][0]*b->values[0][i];
		}
	}
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	else if (a->_rows==1) {
	#else
	else {
	#endif
		if (b->_columns==1) {
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			if (a->_columns!=b->_rows) goto ERROR_MATCH;
			#endif
			for (type_int i=0;i<a->_rows;++i) result += a->values[0][i]*b->values[i][0];
		}
		#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
		else if (b->_rows==1) {
		#else
		else {
		#endif
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			if (a->_columns!=b->_columns) goto ERROR_MATCH;
			#endif
			for (type_int i=0;i<a->_rows;++i) result += a->values[0][i]*b->values[0][i];
		}
	}
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	else goto ERROR_DIM;
	#endif

	return result;

	ERROR_MATCH: throw GLLIB_ERROR_TYPE_ERROR("Cannot compute dot product; dimensions incompatible!");
	ERROR_DIM: throw GLLIB_ERROR_TYPE_ERROR("Cannot computer dot product; matrix does not emulate a 2D array!");
}

//TODO: in-place transpose method, jacobi, inversion, copy_to(GLfloat[16]), operators :-(
//TODO: make work nicely with glLibFraction
//TODO: glLibMath -> glLibExtMath
MATRIX_TEMPLATE class matrix {
	//Represents matrices internally as an array of rows:
	//.values[row][column]

	#pragma region FRIENDS
	friend class glLibInternal_glLibView;
	MATRIX_TEMPLATE2 friend void glLibInternal_check_3x3(const TYPE_MATRIX2* matrix);
	MATRIX_TEMPLATE2 friend void glLibInternal_check_4x4(const TYPE_MATRIX2* matrix);
	MATRIX_TEMPLATE2 friend void glLibInternal_check_3x3or4x4(const TYPE_MATRIX2* matrix);

	MATRIX_TEMPLATE2 friend void       add(const TYPE_MATRIX2* m1, const TYPE_MATRIX2* m2, TYPE_MATRIX2* result);
	MATRIX_TEMPLATE2 friend void  subtract(const TYPE_MATRIX2* m1, const TYPE_MATRIX2* m2, TYPE_MATRIX2* result);
	MATRIX_TEMPLATE2 friend void  multiply(const TYPE_MATRIX2* m1, const TYPE_MATRIX2* m2, TYPE_MATRIX2* result);
	MATRIX_TEMPLATE2 friend void transpose(TYPE_MATRIX2* m);
	MATRIX_TEMPLATE2 friend void transpose(const TYPE_MATRIX2* m, TYPE_MATRIX2* result);

	MATRIX_TEMPLATE2 friend type_element2 determinant(const TYPE_MATRIX2* matrix);
	MATRIX_TEMPLATE2 friend void row_add_multiple(TYPE_MATRIX2* m, type_element2 x, type_int2 r1, type_int2 r2);
	MATRIX_TEMPLATE2 friend void row_swap(TYPE_MATRIX2* m, type_int2 r1, type_int2 r2);
	MATRIX_TEMPLATE2 friend void row_scale(TYPE_MATRIX2* m, type_element2 sc, type_int2 row);
	MATRIX_TEMPLATE2 friend void rref(TYPE_MATRIX2* m);
	MATRIX_TEMPLATE2 friend void invert(const TYPE_MATRIX2* m, TYPE_MATRIX2* result);
	MATRIX_TEMPLATE2 friend void copy(const TYPE_MATRIX2* src, TYPE_MATRIX2* dst);
	MATRIX_TEMPLATE2 friend void copy(
		const TYPE_MATRIX2* src, TYPE_MATRIX2* dst,
		type_int2 src_left, type_int2 src_right, type_int2 src_top, type_int2 src_bottom,
		type_int2 dst_left, type_int2 dst_right, type_int2 dst_top, type_int2 dst_bottom);
	MATRIX_TEMPLATE2 friend void scale(type_element2 scalar,       TYPE_MATRIX2* m);
	MATRIX_TEMPLATE2 friend void scale(type_element2 scalar, const TYPE_MATRIX2* m, TYPE_MATRIX2* result);
	MATRIX_TEMPLATE2 friend void set_rotation_matrix_rad(TYPE_MATRIX2* rot_matrix, const type_element2 vec[3], type_element2 theta);
	MATRIX_TEMPLATE2 friend void set_transform(TYPE_MATRIX2* matrix, const type_element2 new_i[3], const type_element2 new_j[3], const type_element2 new_k[3]);
	MATRIX_TEMPLATE2 friend void conjugate_gradient(const TYPE_MATRIX2* A, const TYPE_MATRIX2* x, const TYPE_MATRIX2* b, int iterations);
	MATRIX_TEMPLATE2 friend type_element2 dot(const TYPE_MATRIX2* a, const TYPE_MATRIX2* b);
	#pragma endregion

	//Heavily adapted from http://www.acims.arizona.edu/PUBLICATIONS/PDF/RajnikanthTemplates04.pdf
	public:
		type_element** values;
	private:
		type_int _rows; //TODO: numof_rows, numof_columns
		type_int _columns;
	private:
		void glLibInternal_init(type_int rows, type_int columns) {
			_rows = rows;
			_columns = columns;

			values = new type_element*[_rows];
			for (type_int i=0;i<_rows;++i) {
				values[i] = new type_element[_columns];
			}
		}
	public:
		//Make an empty, identity, zeroes, or ones, matrix
		matrix(type_int rows, type_int columns, const glLibFlagMatrix flag) {
			glLibInternal_init(rows,columns);
			if      (flag==GLLIB_MATRIX_EMPTY);
			else if (flag==GLLIB_MATRIX_IDENTITY) {
				for (type_int i=0;i<_rows;++i) {
					for (type_int j=0;j<_columns;++j) {
						if (i==j) values[i][j] = 1;
						else values[i][j] = 0;
					}
				}
			}
			else if (flag==GLLIB_MATRIX_ZEROS) { //implicitly includes GLLIB_MATRIX_ZEROES, with "E"
				for (type_int i=0;i<_rows;++i) {
					for (type_int j=0;j<_columns;++j) {
						values[i][j] = 0;
					}
				}
			}
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			else if (flag==GLLIB_MATRIX_ONES) {
			#else
			else {
			#endif
				for (type_int i=0;i<_rows;++i) {
					for (type_int j=0;j<_columns;++j) {
						values[j][j] = 1;
					}
				}
			}
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			else {
				throw error("Matrix type not recognized!");
			}
			#endif
		}
		//Make matrix from data
		matrix(type_int rows, type_int columns, type_element* data, const glLibFlagMatrix row_column_flag) {
			glLibInternal_init(rows,columns);

			if (row_column_flag==GLLIB_MATRIX_ROW) {
				type_int index = 0;
				for (type_int i=0;i<_columns;++i) {
					for (type_int j=0;j<_rows;++j) {
						values[i][j] = data[index++];
					}
				}
			}
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			else if (row_column_flag==GLLIB_MATRIX_COLUMN) {
			#else
			else {
			#endif
				//Matrix creation from array goes down columns instead of across rows.
				//Good for direct and intuitive copying to and from OpenGL.
				type_int index = 0;
				for (type_int i=0;i<_columns;++i) {
					for (type_int j=0;j<_rows;++j) {
						values[j][i] = data[index++];
					}
				}
			}
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			else {
				throw GLLIB_ERROR_TYPE_ERROR("Array interpretation not recognized!  Must be GLLIB_MATRIX_ROW or GLLIB_MATRIX_COLUMN.");
			}
			#endif
		}
		//Make matrix from data (array is interpreted as array of rows, just as internal state)
		matrix(type_int rows, type_int columns, type_element** data) {
			glLibInternal_init(rows,columns);

			for (type_int i=0;i<_rows;++i) {
				for (type_int j=0;j<_columns;++j) {
					values[i][j] = data[i][j];
				}
			}
		}
		//Make an augmented matrix
		//TODO: add overloads these two functions?
		matrix(const TYPE_MATRIX* mat1, const TYPE_MATRIX* mat2) {
			#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
			if (mat1->_rows!=mat2->_rows) throw GLLIB_ERROR_TYPE_ERROR("Row dimensions do not match!  Cannot construct augmented matrix!");
			#endif

			glLibInternal_init(mat1->_rows,mat1->_columns+mat2->_columns);

			for (type_int i=0;i<_rows;++i) {
				for (type_int j=0;j<mat1->_columns;++j) {
					values[i][j] = mat1->values[i][j];
				}
				for (type_int j=mat1->_columns;j<_columns;++j) {
					values[i][j] = mat2->values[i][j-mat1->_columns];
				}
			}
		}
		//Make a submatrix
		matrix(const TYPE_MATRIX* mat, type_int left, type_int right, type_int top, type_int bottom) {
			throw GLLIB_ERROR_TYPE_ERROR("Not implemented!");
		}
		//Copy from a matrix
		matrix(const TYPE_MATRIX& a) {
			glLibInternal_init(a._rows,a._columns);

			for (type_int i=0;i<_rows;++i) {
				for (type_int j=0;j<_columns;++j) {
					values[i][j] = a.values[i][j];
				}
			}
		}
		~matrix(void) {
			for (type_int i=0;i<_rows;++i) {
				delete [] values[i];
			}
			delete [] values;
		}

		void print(void) const {
			//TODO: make #ifdef to use unicode and non-unicode characters
			//TODO: format arguments
			//TODO: newline

			/*printf("\xC9"); //╔ 201, 2554
			for (type_int i=0; i<_columns; ++i) printf("         ");
			printf("\xBB\n"); //╗, 187, 2557*/
			printf("\xDA"); //┌
			for (type_int i=0; i<_columns; ++i) printf("         ");
			printf("\xBF\n"); //┐

			for (type_int i=0; i<_rows; ++i) {
				//printf("\xBA"); //║ 186, 2551
				printf("\xB3"); //│
				for (type_int j=0; j<_columns; ++j) {
					glLibInternal_print_element(values[i][j]);
				}
				//printf("\xBA\n"); //║ 186, 2551
				printf("\xB3\n"); //│
			}

			/*printf("\xC8"); //╚ 200, 255A
			for (type_int i=0; i<_columns; ++i) printf("         ");
			printf("\xBC\n"); //╝ 188, 255D*/
			printf("\xC0"); //└
			for (type_int i=0; i<_columns; ++i) printf("         ");
			printf("\xD9\n"); //┘
		}
	private:
		template <typename type_int2> void glLibInternal_print_element(const ::glLib::Math::Fraction::fraction<type_int2> value) const {
			printf(" (");
			value.glLibInternal_print_numerator();
			printf("/");
			value.glLibInternal_print_denominator();
			printf(") ");
		}
		template <typename type_element2> void glLibInternal_print_element(type_element2 value) const {
			printf("% 8.4f ",value);
			//printf("% f\t",values[i][j]);
		}
	public:

		//TODO: copy_from_OpenGL({GL_MODELVIEW|GL_PROJECTION|GL_TEXTURE|etc.})
		//TODO: copy_to_OpenGL({GL_MODELVIEW|GL_PROJECTION|GL_TEXTURE|etc.})

		//TODO: test these and set_at
		//TODO: row column, not column row (ray tracer uses this incorrectly, change)

		//TODO: overload [] operator
		inline type_element get_at(type_int column, type_int row) const { //column, row, corresponds to x, y
			return values[row][column];
		}
		inline type_element get_at(type_int n) const { //Chooses the nth element (by mathematical standard; reads along rows)
			return values[n/_rows][n%_rows];
		}
		//TODO: set_at
};

template <typename type_int> void glMultMatrixf(const matrix<float,type_int>* matrix) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	glLibInternal_check_4x4(matrix);
	#endif

	GLfloat values[16];
	char index = 0; //OpenGL is backwards from this, hence j, i.
	for (char i=0; i<4; ++i) { for (char j=0; j<4; ++j) { values[index++] = matrix->values[j][i]; } }

	::glMultMatrixf(values);
}
template <typename type_int> void glMultMatrixd(const matrix<double,type_int>* matrix) {
	#ifdef GLLIB_ERROR_MODE_REPORT_ERRORS
	glLibInternal_check_4x4(matrix);
	#endif

	GLdouble values[16];
	char index = 0; //OpenGL is backwards from this, hence j, i.
	for (char i=0; i<4; ++i) { for (char j=0; j<4; ++j) { values[index++] = matrix->values[j][i]; } }

	::glMultMatrixd(values);
}

};};};