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

namespace glLibMathVec {
	void glLibInternal_math_init(void);
	void glLibInternal_math_deinit(void);


	//TODO: Rename fillones, fillzeroes/fillzeros, etc.
	//TODO: Add fill(num,vec) methods.
	//TODO: reciprocal
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num> inline void zero(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 zero2(type_num vec[2]) {
		vec[0] = vec[1] = (type_num)0;
	}
	template <typename type_num> inline void zero3(type_num vec[3]) {
		vec[0] = vec[1] = vec[2] = (type_num)0;
	}
	template <typename type_num> inline void zero4(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 one(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 one2(type_num vec[2]) {
		vec[0] = vec[1] = (type_num)1;
	}
	template <typename type_num> inline void one3(type_num vec[3]) {
		vec[0] = vec[1] = vec[2] = (type_num)1;
	}
	template <typename type_num> inline void one4(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 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, 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, 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, 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, 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];
	}
	#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(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(type_num vec[2], type_num result[2]) {
		result[0] = -vec[0];
		result[1] = -vec[1];
	}
	template <typename type_num> inline void negate3(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(type_num vec[4], type_num result[4]) {
		result[0] = -vec[0];
		result[1] = -vec[1];
		result[2] = -vec[2];
		result[3] = -vec[3];
	}
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num> inline void add(type_num vec1[], 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(type_num vec1[2], 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(type_num vec1[3], 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(type_num vec1[4], 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];
	}
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num> inline void scaleadd(type_num sc, type_num vec1[], 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, type_num vec1[2], 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, type_num vec1[3], 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, type_num vec1[4], 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];
	}
	//TODO: scalesubtract
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num> inline void subtract(type_num vec1[], 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(type_num vec1[2], 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(type_num vec1[3], 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(type_num vec1[4], 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];
	}
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num> inline void multiply(type_num vec1[], 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(type_num vec1[2], 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(type_num vec1[3], 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(type_num vec1[4], 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];
	}
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num> inline void divide(type_num vec1[], 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(type_num vec1[2], 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(type_num vec1[3], 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(type_num vec1[4], 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];
	}


	//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, type_float vec1[], 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, type_float vec1[2], 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, type_float vec1[3], 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, type_float vec1[4], 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];
	}


	#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> type_float length(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(type_float vec[2]) {
		return sqrt( vec[0]*vec[0]+vec[1]*vec[1] );
	}
	template <typename type_float> inline type_float length3(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(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 );
	}

	#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;
		}
		zero3(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(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(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(vec);
		return false;
	}
	template <typename type_float> bool normalize3(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;
		}
		zero3(vec);
		return false;
	}
	template <typename type_float> bool normalize4(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;
		}
		zero4(vec);
		return false;
	}


	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_float> type_float dot(type_float vec1[], 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(type_float vec1[2], type_float vec2[2]) {
		return vec1[0]*vec2[0] + vec1[1]*vec2[1];
	}
	template <typename type_float> inline type_float dot3(type_float vec1[3], 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(type_float vec1[4], type_float vec2[4]) {
		return vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2] + vec1[3]*vec2[3];
	}
	
	
	template <typename type_float> void cross(type_float vec1[3], type_float vec2[3], type_float result[3]) {
		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];
	}


	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_float> type_float min(type_float vec[], int vec_length) {
		type_float minimum = vec[0];
		int index = 1;
		while (index<vec_length) {
			minimum = std::min(minimum,vec[index++]);
		}
		return minimum;
	}
	#endif
	template <typename type_float> inline type_float min2(type_float vec[2]) {
		return std::min(vec[0],vec[1]);
	}
	template <typename type_float> inline type_float min3(type_float vec[3]) {
		return std::min(vec[0],std::min(vec[1],vec[2]));
	}
	template <typename type_float> inline type_float min4(type_float vec[4]) {
		return std::min(std::min(vec[0],vec[1]),std::min(vec[2],vec[3]));
	}
	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_float> type_float max(type_float vec[], int vec_length) {
		type_float maximum = vec[0];
		int index = 1;
		while (index<vec_length) {
			maximum = std::max(maximum,vec[index++]);
		}
		return maximum;
	}
	#endif
	template <typename type_float> inline type_float max2(type_float vec[2]) {
		return std::max(vec[0],vec[1]);
	}
	template <typename type_float> inline type_float max3(type_float vec[3]) {
		return std::max(vec[0],std::max(vec[1],vec[2]));
	}
	template <typename type_float> inline type_float max4(type_float vec[4]) {
		return std::max(std::max(vec[0],vec[1]),std::max(vec[2],vec[3]));
	}


	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num1, typename type_num2> void copy(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(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(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(type_num1 src[4], type_num2 dst[4]) {
		dst[0] = src[0];
		dst[1] = src[1];
		dst[2] = src[2];
		dst[3] = src[3];
	}


	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_num1, typename type_num2> bool equals(type_num1 vec1[], 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(type_num1 vec1[2], 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(type_num1 vec1[3], 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(type_num1 vec1[4], type_num2 vec2[4]) {
		return (vec1[0]==vec2[0]&&vec1[1]==vec2[1]&&vec1[2]==vec2[2]&&vec1[3]==vec2[3]) ? true : false;
	}


	#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);
	}


	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_float> void print(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(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(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(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]);
	}


	#ifndef GLLIBMATHVEC_DISALLOW_GENERICS
	template <typename type_float> void reflect(type_float incident[], 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(type_float incident[2], 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(type_float incident[3], 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(type_float incident[4], 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;
	}


	void rotate_rad(float  point[3], float  result[3], float  vec[3],         float  theta);
	void rotate_rad(double point[3], double result[3], double vec[3],         double theta);
	void rotate_deg(float  point[3], float  result[3], float  vec[3],         float  theta);
	void rotate_deg(double point[3], double result[3], double vec[3],         double theta);
	void rotate(float  point[3], float  result[3], glLibMathMatrix< float,char>* rot_matrix);
	void rotate(double point[3], double result[3], glLibMathMatrix<double,char>* rot_matrix);
};