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

namespace glLib { namespace Math { namespace Vec {

void glLibInternal_math_init(void);
void glLibInternal_math_deinit(void);

#pragma region FILL
//TODO: Test fill(num,vec) methods.
//TODO: reciprocal test
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void fillzeros(type_num vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) { vec[i]=(type_num)0; }
}
#endif
template <typename type_num> inline void fillzeros2(type_num vec[2]) {
	vec[0] = vec[1] = (type_num)0;
}
template <typename type_num> inline void fillzeros3(type_num vec[3]) {
	vec[0] = vec[1] = vec[2] = (type_num)0;
}
template <typename type_num> inline void fillzeros4(type_num vec[4]) {
	vec[0] = vec[1] = vec[2] = vec[3] = (type_num)0;
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void fillones(type_num vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) { vec[i]=(type_num)1; }
}
#endif
template <typename type_num> inline void fillones2(type_num vec[2]) {
	vec[0] = vec[1] = (type_num)1;
}
template <typename type_num> inline void fillones3(type_num vec[3]) {
	vec[0] = vec[1] = vec[2] = (type_num)1;
}
template <typename type_num> inline void fillones4(type_num vec[4]) {
	vec[0] = vec[1] = vec[2] = vec[3] = (type_num)1;
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void fill(type_num num, type_num vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) { vec[i]=num; }
}
#endif
template <typename type_num> inline void fill2(type_num num, type_num vec[2]) {
	vec[0] = vec[1] = num;
}
template <typename type_num> inline void fill3(type_num num, type_num vec[3]) {
	vec[0] = vec[1] = vec[2] = num;
}
template <typename type_num> inline void fill4(type_num num, type_num vec[4]) {
	vec[0] = vec[1] = vec[2] = vec[3] = num;
}
#pragma endregion

#pragma region SCALE
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void scale(type_num sc, type_num vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) { vec[i]=sc*vec[i]; }
}
#endif
template <typename type_num> inline void scale2(type_num sc, type_num vec[2]) {
	vec[0] = sc*vec[0];
	vec[1] = sc*vec[1];
}
template <typename type_num> inline void scale3(type_num sc, type_num vec[3]) {
	vec[0] = sc*vec[0];
	vec[1] = sc*vec[1];
	vec[2] = sc*vec[2];
}
template <typename type_num> inline void scale4(type_num sc, type_num vec[4]) {
	vec[0] = sc*vec[0];
	vec[1] = sc*vec[1];
	vec[2] = sc*vec[2];
	vec[3] = sc*vec[3];
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void scale(type_num sc, const type_num vec[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=sc*vec[i]; }
}
#endif
template <typename type_num> inline void scale2(type_num sc, const type_num vec[2], type_num result[2]) {
	result[0] = sc*vec[0];
	result[1] = sc*vec[1];
}
template <typename type_num> inline void scale3(type_num sc, const type_num vec[3], type_num result[3]) {
	result[0] = sc*vec[0];
	result[1] = sc*vec[1];
	result[2] = sc*vec[2];
}
template <typename type_num> inline void scale4(type_num sc, const type_num vec[4], type_num result[4]) {
	result[0] = sc*vec[0];
	result[1] = sc*vec[1];
	result[2] = sc*vec[2];
	result[3] = sc*vec[3];
}
#pragma endregion
#pragma region NEGATE
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void negate(type_num vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) { vec[i]=-vec[i]; }
}
#endif
template <typename type_num> inline void negate2(type_num vec[2]) {
	vec[0] = -vec[0];
	vec[1] = -vec[1];
}
template <typename type_num> inline void negate3(type_num vec[3]) {
	vec[0] = -vec[0];
	vec[1] = -vec[1];
	vec[2] = -vec[2];
}
template <typename type_num> inline void negate4(type_num vec[4]) {
	vec[0] = -vec[0];
	vec[1] = -vec[1];
	vec[2] = -vec[2];
	vec[3] = -vec[3];
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void negate(const type_num vec[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=-vec[i]; }
}
#endif
template <typename type_num> inline void negate2(const type_num vec[2], type_num result[2]) {
	result[0] = -vec[0];
	result[1] = -vec[1];
}
template <typename type_num> inline void negate3(const type_num vec[3], type_num result[3]) {
	result[0] = -vec[0];
	result[1] = -vec[1];
	result[2] = -vec[2];
}
template <typename type_num> inline void negate4(const type_num vec[4], type_num result[4]) {
	result[0] = -vec[0];
	result[1] = -vec[1];
	result[2] = -vec[2];
	result[3] = -vec[3];
}
#pragma endregion
#pragma region ADD
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void add(const type_num vec1[], const type_num vec2[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=vec1[i]+vec2[i]; }
}
#endif
template <typename type_num> inline void add2(const type_num vec1[2], const type_num vec2[2], type_num result[2]) {
	result[0] = vec1[0] + vec2[0];
	result[1] = vec1[1] + vec2[1];
}
template <typename type_num> inline void add3(const type_num vec1[3], const type_num vec2[3], type_num result[3]) {
	result[0] = vec1[0] + vec2[0];
	result[1] = vec1[1] + vec2[1];
	result[2] = vec1[2] + vec2[2];
}
template <typename type_num> inline void add4(const type_num vec1[4], const type_num vec2[4], type_num result[4]) {
	result[0] = vec1[0] + vec2[0];
	result[1] = vec1[1] + vec2[1];
	result[2] = vec1[2] + vec2[2];
	result[3] = vec1[3] + vec2[3];
}
#pragma endregion
#pragma region SCALEADDSUBTRACT
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void scaleadd(type_num sc, const type_num vec1[], const type_num vec2[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=sc*vec1[i]+vec2[i]; }
}
#endif
template <typename type_num> inline void scaleadd2(type_num sc, const type_num vec1[2], const type_num vec2[2], type_num result[2]) {
	result[0] = sc*vec1[0] + vec2[0];
	result[1] = sc*vec1[1] + vec2[1];
}
template <typename type_num> inline void scaleadd3(type_num sc, const type_num vec1[3], const type_num vec2[3], type_num result[3]) {
	result[0] = sc*vec1[0] + vec2[0];
	result[1] = sc*vec1[1] + vec2[1];
	result[2] = sc*vec1[2] + vec2[2];
}
template <typename type_num> inline void scaleadd4(type_num sc, const type_num vec1[4], const type_num vec2[4], type_num result[4]) {
	result[0] = sc*vec1[0] + vec2[0];
	result[1] = sc*vec1[1] + vec2[1];
	result[2] = sc*vec1[2] + vec2[2];
	result[3] = sc*vec1[3] + vec2[3];
}
#pragma endregion
#pragma region SUBTRACT
//TODO: scalesubtract
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void subtract(const type_num vec1[], const type_num vec2[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=vec1[i]-vec2[i]; }
}
#endif
template <typename type_num> inline void subtract2(const type_num vec1[2], const type_num vec2[2], type_num result[2]) {
	result[0] = vec1[0] - vec2[0];
	result[1] = vec1[1] - vec2[1];
}
template <typename type_num> inline void subtract3(const type_num vec1[3], const type_num vec2[3], type_num result[3]) {
	result[0] = vec1[0] - vec2[0];
	result[1] = vec1[1] - vec2[1];
	result[2] = vec1[2] - vec2[2];
}
template <typename type_num> inline void subtract4(const type_num vec1[4], const type_num vec2[4], type_num result[4]) {
	result[0] = vec1[0] - vec2[0];
	result[1] = vec1[1] - vec2[1];
	result[2] = vec1[2] - vec2[2];
	result[3] = vec1[3] - vec2[3];
}
#pragma endregion
#pragma region MULTIPLY
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void multiply(const type_num vec1[], const type_num vec2[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=vec1[i]*vec2[i]; }
}
#endif
template <typename type_num> inline void multiply2(const type_num vec1[2], const type_num vec2[2], type_num result[2]) {
	result[0] = vec1[0] * vec2[0];
	result[1] = vec1[1] * vec2[1];
}
template <typename type_num> inline void multiply3(const type_num vec1[3], const type_num vec2[3], type_num result[3]) {
	result[0] = vec1[0] * vec2[0];
	result[1] = vec1[1] * vec2[1];
	result[2] = vec1[2] * vec2[2];
}
template <typename type_num> inline void multiply4(const type_num vec1[4], const type_num vec2[4], type_num result[4]) {
	result[0] = vec1[0] * vec2[0];
	result[1] = vec1[1] * vec2[1];
	result[2] = vec1[2] * vec2[2];
	result[3] = vec1[3] * vec2[3];
}
#pragma endregion
#pragma region DIVIDE
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void divide(const type_num vec1[], const type_num vec2[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=vec1[i]/vec2[i]; }
}
#endif
template <typename type_num> inline void divide2(const type_num vec1[2], const type_num vec2[2], type_num result[2]) {
	result[0] = vec1[0] / vec2[0];
	result[1] = vec1[1] / vec2[1];
}
template <typename type_num> inline void divide3(const type_num vec1[3], const type_num vec2[3], type_num result[3]) {
	result[0] = vec1[0] / vec2[0];
	result[1] = vec1[1] / vec2[1];
	result[2] = vec1[2] / vec2[2];
}
template <typename type_num> inline void divide4(const type_num vec1[4], const type_num vec2[4], type_num result[4]) {
	result[0] = vec1[0] / vec2[0];
	result[1] = vec1[1] / vec2[1];
	result[2] = vec1[2] / vec2[2];
	result[3] = vec1[3] / vec2[3];
}
#pragma endregion

#pragma region LERP
//Not using glLibMathMisc::lerp because we can cache 1.0-x here.
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> inline type_float lerp(type_float x, const type_float vec1[], const type_float vec2[], type_float result[], int vec_length) {
	type_float x1 = (type_float)(1) - x;
	for (int i=0;i<vec_length;++i) result[i] = x1*vec1[i] + x*vec2[i];
}
#endif
template <typename type_float> inline type_float lerp2(type_float x, const type_float vec1[2], const type_float vec2[2], type_float result[2]) {
	type_float x1 = (type_float)(1) - x;
	result[0] = x1*vec1[0] + x*vec2[0];
	result[1] = x1*vec1[1] + x*vec2[1];
}
template <typename type_float> inline type_float lerp3(type_float x, const type_float vec1[3], const type_float vec2[3], type_float result[3]) {
	type_float x1 = (type_float)(1) - x;
	result[0] = x1*vec1[0] + x*vec2[0];
	result[1] = x1*vec1[1] + x*vec2[1];
	result[2] = x1*vec1[2] + x*vec2[2];
}
template <typename type_float> inline type_float lerp4(type_float x, const type_float vec1[4], const type_float vec2[4], type_float result[4]) {
	type_float x1 = (type_float)(1) - x;
	result[0] = x1*vec1[0] + x*vec2[0];
	result[1] = x1*vec1[1] + x*vec2[1];
	result[2] = x1*vec1[2] + x*vec2[2];
	result[3] = x1*vec1[3] + x*vec2[3];
}
#pragma endregion

#pragma region ABS
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> inline type_float abs(type_float vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) vec[i] = abs(vec[i]);
}
#endif
template <typename type_float> inline type_float abs2(type_float vec[2]) {
	vec[0] = abs(vec[0]);
	vec[1] = abs(vec[1]);
}
template <typename type_float> inline type_float abs3(type_float vec[3]) {
	vec[0] = abs(vec[0]);
	vec[1] = abs(vec[1]);
	vec[2] = abs(vec[2]);
}
template <typename type_float> inline type_float abs4(type_float vec[4]) {
	vec[0] = abs(vec[0]);
	vec[1] = abs(vec[1]);
	vec[2] = abs(vec[2]);
	vec[3] = abs(vec[3]);
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> inline type_float abs(const type_float vec[], type_float result[], int vec_length) {
	for (int i=0;i<vec_length;++i) vec[i] = abs(vec[i]);
}
#endif
template <typename type_float> inline type_float abs2(const type_float vec[2], type_float result[2]) {
	vec[0] = abs(vec[0]);
	vec[1] = abs(vec[1]);
}
template <typename type_float> inline type_float abs3(const type_float vec[3], type_float result[3]) {
	vec[0] = abs(vec[0]);
	vec[1] = abs(vec[1]);
	vec[2] = abs(vec[2]);
}
template <typename type_float> inline type_float abs4(const type_float vec[4], type_float result[4]) {
	vec[0] = abs(vec[0]);
	vec[1] = abs(vec[1]);
	vec[2] = abs(vec[2]);
	vec[3] = abs(vec[3]);
}
#pragma endregion

#pragma region LENGTH
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> type_float length(const type_float vec[], int vec_length) {
	type_float squared_distance = 0.0;
	for (int i=0;i<vec_length;++i) { squared_distance += vec[i]*vec[i]; }
	return sqrt(squared_distance);
}
#endif
template <typename type_float> inline type_float length2(const type_float vec[2]) {
	return sqrt( vec[0]*vec[0]+vec[1]*vec[1] );
}
template <typename type_float> inline type_float length3(const type_float vec[3]) {
	return sqrt( vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2] );
}
template <typename type_float> inline type_float length4(const type_float vec[4]) {
	return sqrt( vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3] );
}

template <typename type_float> inline type_float length2(type_float a, type_float b) {
	return sqrt( a*a + b*b );
}
template <typename type_float> inline type_float length3(type_float a, type_float b, type_float c) {
	return sqrt( a*a + b*b + c*c );
}
template <typename type_float> inline type_float length4(type_float a, type_float b, type_float c, type_float d) {
	return sqrt( a*a + b*b + c*c + d*d );
}
#pragma endregion
#pragma region NORMALIZE
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> bool normalize(type_float vec[], int vec_length) {
	type_float l = length(vec,vec_length);
	if (l!=0.0) {
		scale<type_float>( ((type_float)(1))/l, vec, vec_length );
		return true;
	}
	for (int i=0;i<vec_length;++i) { vec[i]=0.0; }
	return false;
}
#endif
template <typename type_float> bool normalize2(type_float vec[2]) {
	type_float l = length2(vec);
	if (l!=0.0) {
		scale2<type_float>( ((type_float)(1))/l, vec );
		return true;
	}
	zero2(vec);
	return false;
}
template <typename type_float> bool normalize3(type_float vec[3]) {
	type_float l = length3(vec);
	if (l!=0.0) {
		scale3<type_float>( ((type_float)(1))/l, vec );
		return true;
	}
	fillzeros3(vec);
	return false;
}
template <typename type_float> bool normalize4(type_float vec[4]) {
	type_float l = length4(vec);
	if (l!=0.0) {
		scale4<type_float>( ((type_float)(1))/l, vec );
		return true;
	}
	zero4(vec);
	return false;
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> bool normalize(const type_float vec[], type_float result[], int vec_length) {
	type_float l = length(vec,vec_length);
	if (l!=0.0) {
		scale<type_float>( ((type_float)(1))/l, vec, result, vec_length );
		return true;
	}
	for (int i=0;i<vec_length;++i) { result[i]=0.0; }
	return false;
}
#endif
template <typename type_float> bool normalize2(const type_float vec[2], type_float result[2]) {
	type_float l = length2(vec);
	if (l!=0.0) {
		scale2<type_float>( ((type_float)(1))/l, vec, result );
		return true;
	}
	zero2(result);
	return false;
}
template <typename type_float> bool normalize3(const type_float vec[3], type_float result[3]) {
	type_float l = length3(vec);
	if (l!=0.0) {
		scale3<type_float>( ((type_float)(1))/l, vec, result );
		return true;
	}
	fillzeros3(result);
	return false;
}
template <typename type_float> bool normalize4(const type_float vec[4], type_float result[4]) {
	type_float l = length4(vec);
	if (l!=0.0) {
		scale4<type_float>( ((type_float)(1))/l, vec, result );
		return true;
	}
	fillzeros4(result);
	return false;
}
#pragma endregion

#pragma region DOT
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> type_float dot(const type_float vec1[], const type_float vec2[], int vec_length) {
	type_float result = 0.0;
	for (int i=0;i<vec_length;++i) {
		result += vec1[i]*vec2[i];
	}
	return result;
}
#endif
template <typename type_float> inline type_float dot2(const type_float vec1[2], const type_float vec2[2]) {
	return vec1[0]*vec2[0] + vec1[1]*vec2[1];
}
template <typename type_float> inline type_float dot3(const type_float vec1[3], const type_float vec2[3]) {
	return vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2];
}
template <typename type_float> inline type_float dot4(const type_float vec1[4], const type_float vec2[4]) {
	return vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2] + vec1[3]*vec2[3];
}
#pragma endregion
#pragma region CROSS
template <typename type_float> void cross(const type_float vec1[3], const type_float vec2[3], type_float result[3]) {
	#ifdef GLLIB_ERROR_MODE_REPORT_WARNINGS
	if (vec1==result||vec2==result) warning("Cross product implementation does not operate in-place!");
	#endif
	result[0] = vec1[1]*vec2[2]-vec1[2]*vec2[1];
	result[1] = vec1[2]*vec2[0]-vec1[0]*vec2[2];
	result[2] = vec1[0]*vec2[1]-vec1[1]*vec2[0];
}
#pragma endregion

#pragma region MINMAX
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> type_float min(const type_float vec[], int vec_length) {
	type_float minimum = vec[0];
	int index = 1;
	while (index<vec_length) {
		minimum = ::min(minimum,vec[index++]);
	}
	return minimum;
}
#endif
template <typename type_float> inline type_float min2(const type_float vec[2]) {
	return ::min(vec[0],vec[1]);
}
template <typename type_float> inline type_float min3(const type_float vec[3]) {
	return ::min(vec[0],::min(vec[1],vec[2]));
}
template <typename type_float> inline type_float min4(const type_float vec[4]) {
	return ::min(::min(vec[0],vec[1]),::min(vec[2],vec[3]));
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> type_float max(const type_float vec[], int vec_length) {
	type_float maximum = vec[0];
	int index = 1;
	while (index<vec_length) {
		maximum = ::max(maximum,vec[index++]);
	}
	return maximum;
}
#endif
template <typename type_float> inline type_float max2(const type_float vec[2]) {
	return ::max(vec[0],vec[1]);
}
template <typename type_float> inline type_float max3(const type_float vec[3]) {
	return ::max(vec[0],::max(vec[1],vec[2]));
}
template <typename type_float> inline type_float max4(const type_float vec[4]) {
	return ::max(::max(vec[0],vec[1]),::max(vec[2],vec[3]));
}
#pragma endregion
//TODO: stuff like max(vec1,vec2,resultvec)

#pragma region COPY
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num1, typename type_num2> void copy(const type_num1 src[], type_num2 dst[], int vec_length) {
	for (int i=0;i<vec_length;++i) { dst[i] = src[i]; }
}
#endif
template <typename type_num1, typename type_num2> inline void copy2(const type_num1 src[2], type_num2 dst[2]) {
	dst[0] = src[0];
	dst[1] = src[1];
}
template <typename type_num1, typename type_num2> inline void copy3(const type_num1 src[3], type_num2 dst[3]) {
	dst[0] = src[0];
	dst[1] = src[1];
	dst[2] = src[2];
}
template <typename type_num1, typename type_num2> inline void copy4(const type_num1 src[4], type_num2 dst[4]) {
	dst[0] = src[0];
	dst[1] = src[1];
	dst[2] = src[2];
	dst[3] = src[3];
}
#pragma endregion

#pragma region EQUALS
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num1, typename type_num2> bool equals(const type_num1 vec1[], const type_num2 vec2[], int vec_length) {
	for (int i=0;i<vec_length;++i) { if (vec1[i]!=vec2[i]) return false; }
	return true;
}
#endif
template <typename type_num1, typename type_num2> inline bool equals2(const type_num1 vec1[2], const type_num2 vec2[2]) {
	return (vec1[0]==vec2[0]&&vec1[1]==vec2[1]) ? true : false;
}
template <typename type_num1, typename type_num2> inline bool equals3(const type_num1 vec1[3], const type_num2 vec2[3]) {
	return (vec1[0]==vec2[0]&&vec1[1]==vec2[1]&&vec1[2]==vec2[2]) ? true : false;
}
template <typename type_num1, typename type_num2> inline bool equals4(const type_num1 vec1[4], const type_num2 vec2[4]) {
	return (vec1[0]==vec2[0]&&vec1[1]==vec2[1]&&vec1[2]==vec2[2]&&vec1[3]==vec2[3]) ? true : false;
}
#pragma endregion

#pragma region CLAMP
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> void clamp(type_num vec[], type_num low, type_num high, int vec_length) {
	for (int i=0;i<vec_length;++i) vec[i] = std::max(std::min(vec[i],high),low);
}
#endif
template <typename type_num> inline void clamp2(type_num vec[2], type_num low, type_num high) {
	vec[0] = std::max(std::min(vec[0],high),low);
	vec[1] = std::max(std::min(vec[1],high),low);
}
template <typename type_num> inline void clamp3(type_num vec[3], type_num low, type_num high) {
	vec[0] = std::max(std::min(vec[0],high),low);
	vec[1] = std::max(std::min(vec[1],high),low);
	vec[2] = std::max(std::min(vec[2],high),low);
}
template <typename type_num> inline void clamp4(type_num vec[4], type_num low, type_num high) {
	vec[0] = std::max(std::min(vec[0],high),low);
	vec[1] = std::max(std::min(vec[1],high),low);
	vec[2] = std::max(std::min(vec[2],high),low);
	vec[3] = std::max(std::min(vec[3],high),low);
}
#pragma endregion
//TODO: these with result vectors.  Check this whole module for similar problems (can we add result vector overloads?)

#pragma region RECIPROCAL
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void reciprocal(type_num vec[], int vec_length) {
	for (int i=0;i<vec_length;++i) { vec[i]=(type_num)(1)/vec[i]; }
}
#endif
template <typename type_num> inline void reciprocal2(type_num vec[2]) {
	vec[0] = (type_num)(1)/vec[0];
	vec[1] = (type_num)(1)/vec[1];
}
template <typename type_num> inline void reciprocal3(type_num vec[3]) {
	vec[0] = (type_num)(1)/vec[0];
	vec[1] = (type_num)(1)/vec[1];
	vec[2] = (type_num)(1)/vec[2];
}
template <typename type_num> inline void reciprocal4(type_num vec[4]) {
	vec[0] = (type_num)(1)/vec[0];
	vec[1] = (type_num)(1)/vec[1];
	vec[2] = (type_num)(1)/vec[2];
	vec[3] = (type_num)(1)/vec[3];
}
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_num> inline void reciprocal(const type_num vec[], type_num result[], int vec_length) {
	for (int i=0;i<vec_length;++i) { result[i]=(type_num)(1)/vec[i]; }
}
#endif
template <typename type_num> inline void reciprocal2(const type_num vec[2], type_num result[2]) {
	result[0] = (type_num)(1)/vec[0];
	result[1] = (type_num)(1)/vec[1];
}
template <typename type_num> inline void reciprocal3(const type_num vec[3], type_num result[3]) {
	result[0] = (type_num)(1)/vec[0];
	result[1] = (type_num)(1)/vec[1];
	result[2] = (type_num)(1)/vec[2];
}
template <typename type_num> inline void reciprocal4(const type_num vec[4], type_num result[4]) {
	result[0] = (type_num)(1)/vec[0];
	result[1] = (type_num)(1)/vec[1];
	result[2] = (type_num)(1)/vec[2];
	result[3] = (type_num)(1)/vec[3];
}
#pragma endregion

#pragma region PRINT
#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> void print(const type_float vec[], int vec_length, bool add_newline=true) {
	printf("[");
	for (int i=0;i<vec_length-1;++i) {
		printf("%f, ",vec[i]);
	}
	printf("%f",vec[vec_length-1]);
	(add_newline) ? printf("]\n") : printf("]");
}
#endif
template <typename type_float> inline void print2(const type_float vec[2], bool add_newline=true) {
	(add_newline) ? printf("[%f, %f]\n",vec[0],vec[1]) : printf("[%f, %f]",vec[0],vec[1]);
}
template <typename type_float> inline void print3(const type_float vec[3], bool add_newline=true) {
	(add_newline) ? printf("[%f, %f, %f]\n",vec[0],vec[1],vec[2]) : printf("[%f, %f, %f]",vec[0],vec[1],vec[2]);
}
template <typename type_float> inline void print4(const type_float vec[4], bool add_newline=true) {
	(add_newline) ? printf("[%f, %f, %f, %f]\n",vec[0],vec[1],vec[2],vec[3]) : printf("[%f, %f, %f, %f]",vec[0],vec[1],vec[2],vec[3]);
}
#pragma endregion

#pragma region REFLECTREFRACT
template <typename type_float> inline void reflect_point_plane(type_float point[3], const type_float normal[3]) {
	reflect_point_plane(point,normal[0],normal[1],normal[2]);
}
//TODO: a version also with d
template <typename type_float> void reflect_point_plane(type_float point[3], type_float a, type_float b, type_float c) {
	//http://en.wikipedia.org/wiki/Transformation_matrix#Reflection
	type_float out_x =  (1.0-2.0*a*a)*point[0] -2.0*a*b      *point[1] -2.0*a*c      *point[2];
	type_float out_y = -2.0*a*b      *point[0] +(1.0-2.0*b*b)*point[1] -2.0*b*c      *point[2];
	type_float out_z = -2.0*a*c      *point[0] -2.0*b*c      *point[1] +(1.0-2.0*c*c)*point[2];
	point[0]=out_x; point[1]=out_y; point[2]=out_z;
}

#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> void reflect(const type_float incident[], const type_float normal[], type_float result[], int vec_length) {
	type_float* temp = new type_float[vec_length];
	scale(2*dot(incident,normal,vec_length),normal,temp,vec_length);
	subtract(temp,incident,result,vec_length);
	delete temp;
}
#endif
template <typename type_float> void reflect2(const type_float incident[2], const type_float normal[2], type_float result[2]) {
	type_float temp[2];
	scale2(2*dot2(incident,normal),normal,temp);
	subtract2(temp,incident,result);
}
template <typename type_float> void reflect3(const type_float incident[3], const type_float normal[3], type_float result[3]) {
	type_float temp[3];
	scale3(2*dot3(incident,normal),normal,temp);
	subtract3(temp,incident,result);
}
template <typename type_float> void reflect4(const type_float incident[4], const type_float normal[4], type_float result[4]) {
	type_float temp[4];
	scale4(2*dot4(incident,normal),normal,temp);
	subtract4(temp,incident,result);
}

#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
template <typename type_float> bool refract(type_float incident[], type_float normal[], type_float result[], type_float eta, int vec_length) {
	//http://www.opengl.org/sdk/docs/manglsl/xhtml/refract.xml

	//printf("incident = <%f,%f,%f>\n",incident[0],incident[1],incident[2]);
	//printf("normal = <%f,%f,%f>\n",normal[0],normal[1],normal[2]);
	//printf("eta = %f\n",eta);
	
	//k = 1.0 - eta * eta * (1.0 - dot(N, I) * dot(N, I));
	type_float n_dot_i = -dot(normal,incident,vec_length);
	//printf("n * i = %f\n",n_dot_i);
	type_float k = 1.0 - eta * eta * (1.0 - n_dot_i * n_dot_i);
	//printf("k = %f\n",k);
	if (k<0.0) return false;

	//R = eta * I - (eta * dot(N, I) + sqrt(k)) * N;
	type_float* eta_i = new type_float[vec_length];
	scale(-eta,incident,eta_i,vec_length);
	
	//printf("eta_i = <%f,%f,%f>\n",eta_i[0],eta_i[1],eta_i[2]);
	scale( eta*n_dot_i+sqrt(k), normal, result, vec_length); //result is used as a temp here
	//printf("result = <%f,%f,%f>\n",result[0],result[1],result[2]);
	subtract(eta_i,result,result,vec_length);
	//printf("result = <%f,%f,%f>\n\n",result[0],result[1],result[2]);
	delete eta_i;

	return true;
}
#endif
template <typename type_float> bool refract2(type_float incident[2], type_float normal[2], type_float result[2], type_float eta) {
	type_float n_dot_i = -dot2(normal,incident);
	type_float k = 1.0 - eta * eta * (1.0 - n_dot_i * n_dot_i);
	if (k<0.0) return false;

	type_float eta_i[2];
	scale2(-eta,incident,eta_i);

	scale2( eta*n_dot_i+sqrt(k), normal, result); //result is used as a temp here
	subtract2(eta_i,result,result);

	return true;
}
template <typename type_float> bool refract3(type_float incident[3], type_float normal[3], type_float result[3], type_float eta) {
	type_float n_dot_i = -dot3(normal,incident);
	type_float k = 1.0 - eta * eta * (1.0 - n_dot_i * n_dot_i);
	if (k<0.0) return false;

	type_float eta_i[3];
	scale3(-eta,incident,eta_i);

	scale3( eta*n_dot_i+sqrt(k), normal, result); //result is used as a temp here
	subtract3(eta_i,result,result);

	return true;
}
template <typename type_float> bool refract4(type_float incident[4], type_float normal[4], type_float result[4], type_float eta) {
	type_float n_dot_i = -dot4(normal,incident);
	type_float k = 1.0 - eta * eta * (1.0 - n_dot_i * n_dot_i);
	if (k<0.0) return false;

	type_float eta_i[4];
	scale4(-eta,incident,eta_i);

	scale4( eta*n_dot_i+sqrt(k), normal, result); //result is used as a temp here
	subtract4(eta_i,result,result);

	return true;
}
#pragma endregion

#pragma region ROTATE
void rotate_rad(const  float point[3],  float result[3], const  float vec[3], float  theta);
void rotate_rad(const double point[3], double result[3], const double vec[3], double theta);
void rotate_deg(const  float point[3],  float result[3], const  float vec[3], float  theta);
void rotate_deg(const double point[3], double result[3], const double vec[3], double theta);
void     rotate(const  float point[3],  float result[3], const Math::Matrix::matrix< float,char>* rot_matrix);
void     rotate(const double point[3], double result[3], const Math::Matrix::matrix<double,char>* rot_matrix);
#pragma endregion

};};};