gistfile1.d

module toybox.math.vector;
private 
{
    import std.stdio, std.string, std.metastrings, std.typecons, std.traits;
    static import std.math, std.algorithm;
	import toybox.math.matrix;
}

template IsComponentType(T) { static if(is(T == float) || is(T == double) || is(T == int) || is(T == uint) || is(T == bool)) { enum IsComponentType = true; } else { enum IsComponentType = false; } }
template IsFloat(T) { static if(is(T == float) || is(T == double)) { enum IsFloat = true; } else { enum IsFloat = false; } }
template IsInteger(T) { static if(is(T == int) || is(T == uint)) { enum IsInteger = true; } else { enum IsInteger = false; } }
template IsScalar(T) { static if(is(T == float) || is(T == double) || is(T == int) || is(T == uint)) { enum IsScalar = true; } else { enum IsScalar = false; } }


struct VectorTN(T, int N) 
	if(IsComponentType!T && N > 1 && N <= 4)
{
	alias N dimensions;
    alias T memberType;


    /** Vector members as N dimensional static array of type T, by default initialize members to 0 */
    T[N] members = 0;



	/** Create a vector with every member set to value specified */
	this(U)(U value) pure nothrow
		if(IsComponentType!U)
	{
		auto valueT = cast(T)value;
		foreach(i; 0 .. N)
			members[i] = valueT;
	}

	/** Construct a vector with every member value provided in values */
	this(U)(const U[] values) pure nothrow
		if(IsComponentType!U)
    {
        foreach(i; 0 .. N < values.length ? N : values.length)
            members[i] = cast(T)values[i];
    }

    /** Construct a vector with members copied from other vector (cast to T) or 0 if other vector is smaller than this */
	this(U, int M)(const auto ref VectorTN!(U, M) other) pure nothrow
		if(IsComponentType!U)
    {
        foreach(i; 0 .. N < M ? N : M)
            members[i] = cast(T)other[i];
    }

	/** Construct a vector with two component values (rest are 0 if vector dimensions > 2) */
	this(U1, U2)(const U1 x, const U2 y) pure nothrow
		if(IsComponentType!U1 && IsComponentType!U2)
    {
        members[0] = cast(T)x;
        members[1] = cast(T)y;
    }

	static if(N > 2)
	{
		/** Construct a vector with three component values (rest are 0 if vector dimensions > 3) */
		this(U1, U2, U3)(const U1 x, const U2 y, const U3 z) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2 && IsComponentType!U3)
		{
			members[0] = cast(T)x;
			members[1] = cast(T)y;
			members[2] = cast(T)z;
		}

		/** Construct a vector with three component values (rest are 0 if vector dimensions > 3) */
		this(U1, U2)(const VectorTN!(U1, 2) xy, const U2 z) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2)
		{
			members[0] = cast(T)xy[0];
			members[1] = cast(T)xy[1];
			members[2] = cast(T)z;
		}

		/** Construct a vector with three component values (rest are 0 if vector dimensions > 3) */
		this(U1, U2)(const U1 x, const VectorTN!(U2, 2) yz) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2)
		{
			members[0] = cast(T)x;
			members[1] = cast(T)yz[0];
			members[2] = cast(T)yz[1];
		}
	}


	static if(N > 3)
	{
		/** Construct a vector with four component values */
		this(U1, U2, U3, U4)(const U1 x, const U2 y, const U3 z, const U4 w) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2 && IsComponentType!U3 && IsComponentType!U4)
		{
			members[0] = cast(T)x;
			members[1] = cast(T)y;
			members[2] = cast(T)z;
			members[3] = cast(T)w;
		}

		/** Construct a vector with four component values */
		this(U1, U2)(const auto ref VectorTN!(U1, 2) xy, const auto ref VectorTN!(U2, 2) zw) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2)
		{
			members[0] = cast(T)xy[0];
			members[1] = cast(T)xy[1];
			members[2] = cast(T)zw[0];
			members[3] = cast(T)zw[1];
		}

		/** Construct a vector with four component values */
		this(U1, U2, U3)(const auto ref VectorTN!(U1, 2) xy, const U2 z, const U3 w) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2 && IsComponentType!U3)
		{
			members[0] = cast(T)xy[0];
			members[1] = cast(T)xy[1];
			members[2] = cast(T)z;
			members[3] = cast(T)w;
		}

		/** Construct a vector with four component values */
		this(U1, U2, U3)(const U1 x, const auto ref VectorTN!(U2, 2) yz, const U3 w) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2 && IsComponentType!U3)
		{
			members[0] = cast(T)x;
			members[1] = cast(T)yz[0];
			members[2] = cast(T)yz[1];
			members[3] = cast(T)w;
		}

		/** Construct a vector with four component values */
		this(U1, U2, U3)(const U1 x, const U2 y, const auto ref VectorTN!(U3, 2) zw) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2 && IsComponentType!U3)
		{
			members[0] = cast(T)x;
			members[1] = cast(T)y;
			members[2] = cast(T)zw[0];
			members[3] = cast(T)zw[1];
		}

		/** Construct a vector with four component values */
		this(U1, U2)(const auto ref VectorTN!(U1, 3) xyz, const U2 w) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2)
		{
			members[0] = cast(T)xyz[0];
			members[1] = cast(T)xyz[1];
			members[2] = cast(T)xyz[2];
			members[3] = cast(T)w;
		}

		/** Construct a vector with four component values */
		this(U1, U2)(const U1 x, const auto ref VectorTN!(U2, 3) yzw) pure nothrow
			if(IsComponentType!U1 && IsComponentType!U2)
		{
			members[0] = cast(T)x;
			members[1] = cast(T)yzw[0];
			members[2] = cast(T)yzw[1];
			members[3] = cast(T)yzw[2];
		}
	}


	/** Return number of dimensions for this vector, implemented as a function this to be compatible with GLSL */
	@property int length() pure const nothrow { return N; }


    /** Index member access operator, simply returns a ref to members[index] */
	ref T opIndex(size_t index) pure nothrow { return members[index]; }

	/** Index member access operator, simply returns a const ref to members[index] */
	ref const(T) opIndex(size_t index) pure const nothrow { return members[index]; }


    /** Returns a pointer to the first element in vector */
    @property T* ptr() pure { return members.ptr; }

    /** Returns a const pointer to the first element in vector */
    @property const(T)* ptr() pure const { return members.ptr; }
    

	/** Binary functions only make sense when T is scalar */
	static if(IsScalar!T)
	{
		/** Binary vector with scalar operator, returns result of op applied to this with scalar for each member in a new vector */
		VectorTN!(T, N) opBinary(string op)(const T scalar) pure const nothrow
		{
			VectorTN!(T, N) result;
			foreach(i; 0 .. N)
				result[i] = mixin("members[i]" ~ op ~ "scalar");
			return result;
		}

		/** Binary vector with vector operator, returns the result of this op other for each member in a new vector */
		VectorTN!(T, N) opBinary(string op)(const auto ref VectorTN!(T, N) other) pure const nothrow
		{
			VectorTN!(T, N) result;
			foreach(i; 0 .. N)
				result[i] = mixin("members[i]" ~ op ~ "other[i]");
			return result;
		}

		/** Binary vector with matrix multiplication operator, implementing left side multiplication (row vector * matrix) */
		VectorTN!(T, N) opBinary(string op)(const auto ref MatrixTN!(T, N) other) pure const nothrow
			if(op == "*")
		{
			VectorTN!(T, N) result;
			foreach(i; 0 .. N)
				result[i] = dot(this, other.columns[i]);
			return result;
		}
	}


    // Internal template implementation details for opDispatch member accessors
    private
    {
        static int coordIndex(char c)() pure nothrow
        {
			switch(c)
			{
				case 'x': return 0;
				case 'y': return 1;
				case 'z': return N > 2 ? 2 : -1;
				case 'w': return N > 3 ? 3 : -1;
				default: 
					return -1;
			}
        }

		static bool coordValid(char c)() pure nothrow { return coordIndex!c != -1; }

        static bool coordsValidGetter(string coords)() pure nothrow
			if(coords.length > 0 && coords.length <= 4)
        {
			foreach(i; 0 .. coords.length)
				if(!coordValid!(coords[i])) 
					return false;
			return true;
        }

		static bool coordsValidSetter(string coords)() pure nothrow
		{
			if(!coordsValidGetter!(coords))
				return false;
			switch(coords.length)
			{
			    case 1: return true;
				case 2: return coords[0] != coords[1];
				case 3: return N > 2 && coords[0] != coords[1] && coords[0] != coords[2] && coords[1] != coords[2];
				case 4: return N > 3 && coords[0] != coords[1] && coords[0] != coords[2] && coords[0] != coords[3] && coords[1] != coords[2] && coords[1] != coords[3] && coords[2] != coords[3];
				default: return false;
			}
		}
    }


	/** Single member getter returns value ref not 1d vector */
    @property ref T opDispatch(string coords)() pure nothrow
        if(coords.length == 1 && coordValid!(coords[0]))
	{
		return members[coordIndex!(coords[0])];
	}

	/** Const version of single member getter returns value const ref not 1d vector */
    @property ref const(T) opDispatch(string coords)() pure const nothrow
        if(coords.length == 1 && coordValid!(coords[0]))
	{
		return members[coordIndex!(coords[0])];
	}


    /** Template swizzling property access */
    @property VectorTN!(T, coords.length) opDispatch(string coords)() pure const nothrow
        if(coords.length > 1 && coords.length <= 4 && coordsValidGetter!coords)
    {
        VectorTN!(T, coords.length) result;
        result[0] = members[coordIndex!(coords[0])];
        result[1] = members[coordIndex!(coords[1])];
        static if(coords.length > 2)
            result[2] = members[coordIndex!(coords[2])];
        static if(coords.length > 3)
            result[3] = members[coordIndex!(coords[3])];
        return result;
    }

    /** Subvector value setter, eg v4.xz = 1 will set x and z to 1, also implements named member setter (ie. vec.x = value, vec.y = value ...) */
    @property void opDispatch(string coords, U)(const U value) pure nothrow
        if(coords.length > 1 && coords.length <= 4 && coordsValidSetter!coords)
    {
        members[coordIndex!(coords[0])] = cast(T)value;
        static if(coords.length > 1)
            members[coordIndex!(coords[1])] = cast(T)value;
        static if(coords.length > 2)
            members[coordIndex!(coords[2])] = cast(T)value;
        static if(coords.length > 3)
            members[coordIndex!(coords[3])] = cast(T)value;
    }

    /** Swizzling copy setter eg. v3.xz = v3.zx will swap x and z,  */
    @property void opDispatch(string coords, U)(const VectorTN!(U, coords.length) other) pure nothrow
        if(coords.length > 1 && coords.length <= 4 && coordsValidSetter!coords)
    {
        members[coordIndex!(coords[0])] = cast(T)other[0];
        members[coordIndex!(coords[1])] = cast(T)other[1];
        static if(coords.length > 2)
            members[coordIndex!(coords[2])] = cast(T)other[2];
        static if(coords.length > 3)
            members[coordIndex!(coords[3])] = cast(T)other[3];
    }
}



//
// Alias VectorTN instances to corresponding GLSL types
//

alias VectorTN!(float, 2) vec2;
alias VectorTN!(float, 3) vec3;
alias VectorTN!(float, 4) vec4;

alias VectorTN!(double, 2) dvec2;
alias VectorTN!(double, 3) dvec3;
alias VectorTN!(double, 4) dvec4;

alias VectorTN!(int, 2) ivec2;
alias VectorTN!(int, 3) ivec3;
alias VectorTN!(int, 4) ivec4;

alias VectorTN!(uint, 2) uivec2;
alias VectorTN!(uint, 3) uivec3;
alias VectorTN!(uint, 4) uivec4;

alias VectorTN!(bool, 2) bvec2;
alias VectorTN!(bool, 3) bvec3;
alias VectorTN!(bool, 4) bvec4;



//
// GLSL functions that operate on vectors and scalars
//

// Internal implementation templates that generate functions that apply Fn to vectors per member and return a result vector
private
{
    mixin template FnX(alias Fn, alias TCons, bool Pure = true)
    {
        static if(Pure) 
        {
            T fnX(T)(const T x) nothrow pure if(TCons!T) { return Fn(x); }

            VectorTN!(T, N) fnX(T, int N)(const auto ref VectorTN!(T, N) x) nothrow pure
            {
                VectorTN!(T, N) result;
                foreach(i; 0 .. N)
                    result[i] = Fn(x[i]);
                return result;
            }
        }
        else 
        {
            T fnX(T)(const T x) nothrow if(TCons!T) { return Fn(x); }

            VectorTN!(T, N) fnX(T, int N)(const auto ref VectorTN!(T, N) x) nothrow
            {
                VectorTN!(T, N) result;
                foreach(i; 0 .. N)
                    result[i] = Fn(x[i]);
                return result;
            }
        }
    }


    mixin template FnXY(alias Fn, alias TCons, bool Pure = true)
    {
        static if(Pure) 
        {
            T fnXY(T) (const T x, const T y) pure nothrow if(TCons!T) { return Fn(x, y); }

            VectorTN!(T, N) fnXY(T, int N)(const auto ref VectorTN!(T, N) x, const T y) pure nothrow  
                if(TCons!T)
            {
                VectorTN!(T, N) result;
                foreach(i; 0 .. N)
                    result[i] = Fn(x[i], y);
                return result;
            }

            VectorTN!(T, N) fnXY(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y) pure nothrow  
                if(TCons!T)
            {
                VectorTN!(T, N) result;
                foreach(i; 0 .. N)
                    result[i] = Fn(x[i], y[i]);
                return result;
            }
        }
        else
        {
            T fnXY(T) (const T x, const T y) nothrow if(TCons!T) { return Fn(x, y); }

            VectorTN!(T, N) fnXY(T, int N)(const auto ref VectorTN!(T, N) x, const T y) nothrow  
                if(TCons!T)
            {
                VectorTN!(T, N) result;
                foreach(i; 0 .. N)
                    result[i] = Fn(x[i], y);
                return result;
            }

            VectorTN!(T, N) fnXY(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y) nothrow  
                if(TCons!T)
            {
                VectorTN!(T, N) result;
                foreach(i; 0 .. N)
                    result[i] = Fn(x[i], y[i]);
                return result;
            }
        }
    }
}



// Wrappers for most std.math functions are generated using macros 
mixin FnX!(std.math.sin, IsFloat) Sin; alias Sin.fnX sin;
mixin FnX!(std.math.cos, IsFloat) Cos; alias Cos.fnX cos;
mixin FnX!(std.math.tan, IsFloat) Tan; alias Tan.fnX tan;
mixin FnX!(std.math.asin, IsFloat) ASin; alias ASin.fnX asin;
mixin FnX!(std.math.acos, IsFloat) ACos; alias ACos.fnX acos;
mixin FnX!(std.math.atan, IsFloat) ATan; alias ATan.fnX atan;
mixin FnX!(std.math.sinh, IsFloat) SinH; alias SinH.fnX sinh;
mixin FnX!(std.math.cosh, IsFloat) CosH; alias CosH.fnX cosh;
mixin FnX!(std.math.tanh, IsFloat) TanH; alias TanH.fnX tanh;
mixin FnX!(std.math.asinh, IsFloat) ASinH; alias ASinH.fnX asinh;
mixin FnX!(std.math.acosh, IsFloat) ACosH; alias ACosH.fnX acosh;
mixin FnX!(std.math.atanh, IsFloat) ATanH; alias ATanH.fnX atanh;
mixin FnXY!(std.math.pow, IsFloat) Pow; alias Pow.fnXY pow;
mixin FnX!(std.math.exp, IsFloat) Exp; alias Exp.fnX exp;
mixin FnX!(std.math.log, IsFloat) Log; alias Log.fnX log;
mixin FnX!(std.math.exp2, IsFloat) Exp2; alias Exp2.fnX exp2;
mixin FnX!(std.math.log2, IsFloat) Log2; alias Log2.fnX log2;
mixin FnX!(std.math.sqrt, IsFloat) Sqrt; alias Sqrt.fnX sqrt;
mixin FnX!(std.math.abs, IsScalar) Abs; alias Abs.fnX abs;
mixin FnX!(std.math.sgn, IsScalar) Sign; alias Sign.fnX sign;
mixin FnX!(std.math.floor, IsFloat, false) Floor; alias Floor.fnX floor;
mixin FnX!(std.math.trunc, IsFloat, false) Trunc; alias Trunc.fnX trunc;
mixin FnX!(std.math.round, IsFloat, false) Round; alias Round.fnX round; alias Round.fnX roundEven;
mixin FnX!(std.math.ceil, IsFloat) Ceil; alias Ceil.fnX ceil;
mixin FnXY!(std.algorithm.min, IsScalar) Min; alias Min.fnXY min;
mixin FnXY!(std.algorithm.max, IsScalar) Max; alias Max.fnXY max;
mixin FnXY!(std.algorithm.fmod, IsScalar, false) Mod; alias Mod.fnXY mod;



/** Returns x - floor(x) */
T fract(T)(const T x) nothrow if(IsFloat!T) { return x - floor(x); }

/** Returns x - floor(x) */
VectorTN!(T, N) fract(T, int N)(const auto ref VectorTN!(T, N) x) nothrow if(IsFloat!T) { return x - floor(x); }


/** Template generalization for std.math.modf, returns integer part of x in ip and fraction in return value */
T modf(T)(const T x, out T ip) nothrow 
    if(IsFloat!T)
{
    real result, tmpip;
    result = std.math.modf(x, tmpip);
    ip = tmpip;
    return result;
}

/** Vector version of modf */
VectorTN!(T, N) modf(T, int N)(const auto ref VectorTN!(T, N) x, out VectorTN!(T, N) ip) nothrow
    if(IsFloat!T)
{
    VectorTN!(T, N) result;
    foreach(i; 0 .. N)
        result[i] = modf(x[i], ip[i]);
    return result;
}


/** Returns min(max(x, minVal), maxVal) */
VectorTN!(T, N) clamp(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) minVal, const auto ref VectorTN!(T, N) maxVal) pure nothrow
	if(IsScalar!T)
{
	return min(max(x, minVal), maxVal);
}

/** Returns min(max(x, minVal), maxVal) */
VectorTN!(T, N) clamp(T, int N)(const auto ref VectorTN!(T, N) x, const T minVal, const T maxVal) pure nothrow
	if(IsScalar!T)
{
	return min(max(x, minVal), maxVal);
}


/** Return a linear blend of x and y based on a, x * (1 - a) + y * a */
VectorTN!(T, N) mix(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y, const auto ref VectorTN!(T, N) a) pure nothrow
	if(IsFloat!T)
{
	return x * (VectorTN!(T, N)(1) - a) - y * a;
}

/** Return a linear blend of x and y based on a, x * (1 - a) + y * a */
VectorTN!(T, N) mix(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y, const T a) pure nothrow
	if(IsFloat!T)
{
	return x * (1 - a) - y * a;
}

/** Selects components from x and y based on values in a, if a[i] is false then result[i] = x[i], otherwise it's y[i] */
VectorTN!(T, N) mix(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y, const auto ref VectorTN!(bool, N) a) pure nothrow
	if(IsScalar!T)
{
	VectorTN!(T, N) result;
	foreach(i; 0 .. N)
		result[i] = a[i] ? y[i] : x[i];
	return result;
}


/** Returns 0 if x < edge, 1 otherwise */
VectorTN!(T, N) step(T, int N)(const auto ref VectorTN!(T, N) edge, const auto ref VectorTN!(T, N) x) pure nothrow
	if(IsScalar!T)
{
	VectorTN!(T, N) result;
	foreach(i; 0 .. N)
		result[i] = x[i] < edge[i] ? 0 : 1;
	return result;
}

/** Returns 0 if x < edge, 1 otherwise */
VectorTN!(T, N) step(T, int N)(const T edge, const auto ref VectorTN!(T, N) x) pure nothrow
	if(IsScalar!T)
{
	VectorTN!(T, N) result;
	foreach(i; 0 .. N)
		result[i] = x[i] < edge ? 0 : 1;
	return result;
}


/** Returns 0.0 if x <= edge0 and 1.0 if x >= edge1 and performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. */
VectorTN!(T, N) smoothstep(T, int N)(const auto ref VectorTN!(T, N) edge0, const auto ref VectorTN!(T, N) edge1, const auto ref VectorTN!(T, N) x) pure nothrow
	if(IsFloat!T)
{
	auto t = clamp((x - edge0) / (edge1 - edge0), 0, 1);
	return t * t * (VectorTN!(T, N)(3) - t * 2);
}

/** Returns 0.0 if x <= edge0 and 1.0 if x >= edge1 and performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. */
VectorTN!(T, N) smoothstep(T, int N)(const T edge0, const T edge1, const auto ref VectorTN!(T, N) x) pure nothrow
	if(IsFloat!T)
{
	auto t = clamp((x - edge0) / (edge1 - edge0), 0, 1);
	return t * t * (VectorTN!(T, N)(3) - t * 2);
}


/** Returns true if x is nan, false otherwise */
bool isnan(T)(const T x) pure nothrow if(IsFloat!T) { return std.math.isNaN(x); }

/** Return member is true if input member is nan, false otherwise. */
VectorTN!(bool, N) isnan(T, int N)(const auto ref VectorTN!(T, N) x) pure nothrow
	if(IsFloat!T) 
{
	VectorTN!(bool, N) result;
	foreach(i; 0 .. N)
		result[i] = isnan(x[0]);
	return result;
}


/** Return true if x is infinity (positive or negative), false otherwise */
bool isinfinity(T)(const T x) pure nothrow if(IsFloat!T) { return std.math.isInfinity(x); }

/** Return member is true if input member is infinity (positive or negative), false otherwise. */
VectorTN!(bool, N) isnan(T, int N)(const auto ref VectorTN!(T, N) x) pure nothrow
	if(IsFloat!T) 
{
	VectorTN!(bool, N) result;
	foreach(i; 0 .. N)
		result[i] = isinfinity(x[0]);
	return result;
}


/** Internal reinterpret cast template for vectors */
private
{
	mixin template reinterpretCast(T, U)
		if(IsScalar!T && IsScalar!U)
	{
		U reinterpretCast()(const T x) pure nothrow { return *cast(U*)&x; }
		VectorTN!(U, N) reinterpretCast(int N)(const auto ref VectorTN!(T, N) x) pure nothrow			
		{
			VectorTN!(U, N) result;
			foreach(i; 0 .. N)
				result[i] = reinterpretCast(x[i]);
			return result;
		}
	}
}

/** Reinterpret cast instances defined by GLSL  */
mixin reinterpretCast!(float, int) FloatBitsToInt; alias FloatBitsToInt.reinterpretCast floatBitsToInt;
mixin reinterpretCast!(float, uint) FloatBitsToUInt; alias FloatBitsToUInt.reinterpretCast floatBitsToUInt;
mixin reinterpretCast!(int, float) IntBitsToFloat; alias IntBitsToFloat.reinterpretCast intBitsToFloat;
mixin reinterpretCast!(uint, float) UIntBitsToFloat; alias UIntBitsToFloat.reinterpretCast uintBitsToFloat;


/** Multiply add, multiply x and y and add z returning result */
T fma(T)(const T a, const T b, const T c) pure nothrow if(IsFloat!T) { return a * b + c; }

/** Multiply add, multiply x and y and add z returning result */
T fma(T, int N)(const auto ref VectorTN!(T, N) a, const auto ref VectorTN!(T, N) b, const auto ref VectorTN!(T, N) c) pure nothrow if (IsFloat!T) { return a * b + c; }


/** Split x in to a significand in range [0.5, 1) and a int exponent of two such that : x = significand * 2 ^ exp */
T frexp(T)(const T x, out int exp) pure nothrow if(IsFloat!T) { return std.math.frexp(x, exp); }

/** Vector component version of frexp */
VectorTN!(T, N) frexp(T, int N)(const auto ref VectorTN!(T, N) x, out VectorTN!(int, N) exp) pure nothrow 
	if(IsFloat!T) 
{
	VectorTN!(T, N) result;
	foreach(i; 0 .. N)
		result[i] = frexp(x[i], exp[i]);
	return result;
}

/** Joins the results of frexp return original valie, ie. returns x * 2^exp */
T ldexp(T)(const T x, const int exp) pure nothrow if(IsFloat!T) { return std.math.ldexp(x, exp); }

/** Vector component version of ldexp */
VectorTN!(T, N) ldexp(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(int, N) exp) pure nothrow 
	if(IsFloat!T)  
{
	VectorTN!(T, N) result;
	foreach(i; 0 .. N)
		result[i] = ldexp(x[i], exp[i]);
	return result;
}


/** Convert degrees to radians */
T radians(T)(const T deg) pure nothrow if(IsFloat!T) { return deg * (std.math.PI / 180); }

/** Convert degrees to radians */
VectorTN!(T, N) radians(T, int N)(const auto ref VectorTN!(T, N) deg) pure nothrow if(IsFloat!T) { return deg * (std.math.PI / 180); }


/** Convert radians to degrees */
T degrees(T)(const T rad) pure nothrow if(IsFloat!T) { return rad * (180 / std.math.PI); }

/** Convert radians to degrees */
VectorTN!(T, N) degrees(T, int N)(const auto ref VectorTN!(T, N) rad) pure nothrow if(IsFloat!T) { return rad * (180 / std.math.PI); }


/** Dot product of two vectors */
T dot(T, int N) (const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y) pure nothrow
    if(IsFloat!T)
{
    T result = 0;
    foreach(i; 0 .. N)
        result += x[i] * y[i];
    return result;
}

/** Vector length */
T length(T, int N) (const auto ref VectorTN!(T, N) x) pure nothrow 
	if(IsFloat!T) 
{ 
	return sqrt(dot(vec, vec)); 
}


/** Distance between p0 and p1, equivalent to length(p0 - p1) */
T distance(T, int N) (const auto ref VectorTN!(T, N) p0, const auto ref VectorTN!(T, N) p1) pure nothrow 
	if(IsFloat!T) 
{ 
	return length(p0 - p1); 
}


/** Return a vector perpendicular to x and y */
VectorTN!(T, 3) cross(T) (const auto ref VectorTN!(T, 3) x, const auto ref VectorTN!(T, 3) y) pure nothrow 
    if(IsFloat!T)
{
    return VectorTN!(T, 3)(
        x.y * y.z - x.z * y.y,
        x.z * y.x - x.x * y.z,
        x.x * y.y - x.y * y.x);
}


/** Return a vector in same direction as x but with length of 1 */
VectorTN!(T, N) normalize(T, int N) (const auto ref VectorTN!(T, N) x) pure nothrow 
	if (IsFloat!T) 
{ 
	return x * (1 / length(x)); 
}


/** If dot(nref, ind) < 0 return norm, otherwise return –norm. */
VectorTN!(T, N) faceforward(T, int N)(const auto ref VectorTN!(T, N) norm, const auto ref VectorTN!(T, N) incd, const auto ref VectorTN!(T, N) nref) pure nothrow
    if(IsFloat!T)
{
    if(dot(nref, incd) < 0)
        return norm;
    else
        return -norm;
}


/** Reflection of incd vector arround normal */
VectorTN!(T, N) reflect(T, int N)(const auto ref VectorTN!(T, N) incd, const auto ref VectorTN!(T, N) norm) pure nothrow
    if(IsFloat!T)
{
    return incd - norm * dot(norm, incd) * 2;
}


/** Refract */
VectorTN!(T, N) refract(T, int N)(const auto ref VectorTN!(T, N) incd, const auto ref VectorTN!(T, N) norm, const T eta) pure nothrow
    if(IsFloat!T)
{
    T niprod = dot(norm, incd);
    T k = 1 - eta * eta * (1 - niprod * niprod);
    if(k < 0)
        return VectorTN!(T, N)(0);
    else
        return incd * eta - norm * (eta * niprod + sqrt(k));
}


/** Implement component based boolean functions using templates */
private 
{
	template boolFnXY(string op, alias TCons)
	{
		bool boolFnXY(T)(const T x, const T y) pure nothrow
			if(TCons!T)
		{
			return mixin("x " ~ op ~ " y");
		}

		VectorTN!(bool, N) boolFnXY(T, int N)(const auto ref VectorTN!(T, N) x, const auto ref VectorTN!(T, N) y) pure nothrow
			if(TCons!T)
		{
			VectorTN!(bool, N) result;
			foreach(i; 0 .. N)
				result[i] = boolFnXY(x[i], y[i]);
			return result;
		}

		VectorTN!(bool, N) boolFnXY(T, int N)(const auto ref VectorTN!(T, N) x, const T y) pure nothrow
			if(TCons!T)
		{
			VectorTN!(bool, N) result;
			foreach(i; 0 .. N)
				result[i] = boolFnXY(x[i], y);
			return result;
		}
	}
}

mixin boolFnXY!("<", IsScalar) LessThan; alias LessThan.boolFnXY lessThan;
mixin boolFnXY!("<=", IsScalar) LessThanEqual; alias LessThanEqual.boolFnXY lessThanEqual;
mixin boolFnXY!(">", IsScalar) GreaterThan; alias GreaterThan.boolFnXY greaterThan;
mixin boolFnXY!(">=", IsScalar) GreaterThanEqual; alias GreaterThanEqual.boolFnXY greaterThanEqual;
mixin boolFnXY!("==", IsScalar) Equal; alias Equal.boolFnXY equal;
mixin boolFnXY!("!=", IsScalar) NotEqual; alias NotEqual.boolFnXY notEqual;


/** Returns logical complement of x */
bool not()(bool x) pure nothrow { return !x; }

/** Returns the component-wise logical complement of x. */
VectorTN!(bool, N) not(int N)(const auto ref VectorTN!(bool, N) x) pure nothrow
{
	VectorTN!(bool, N) result;
	foreach(i; 0 .. N)
		result[i] = !x[i];
	return result;
}


/** Return true if any of vector components are true */
bool any(int N)(const auto ref VectorTN!(bool, N) x) pure nothrow
{
	bool result = false;
	foreach(i; 0 .. N)
		result |= x[i];
	return result;
}

/** Returns true only if all components of x are true. */
bool all(int N)(const auto ref VectorTN!(bool, N) x) pure nothrow
{
	bool result = true;
	foreach(i; 0 .. N)
		result &= x[i];
	return result;
}