#ifndef __LIBMATH_H__
#define __LIBMATH_H__

#include <math.h>
#include <algorithm>
using std::swap;

#include "core.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#endif

#define EPSILON 1e-6
#define INFINITY 1e30

// single precision mathematical functions

#if _MSC_VER < 1310
inline float ceil( float x ) { return ceilf( x ); }
inline float floor( float x ) { return floorf( x ); }
inline float fabs( float x ) { return fabsf( x ); }
inline float sqrt( float x ) { return sqrtf( x ); }
inline float sin( float x ) { return sinf( x ); }
inline float cos( float x ) { return cosf( x ); }
inline float tan( float x ) { return tanf( x ); }
inline float acos( float x ) { return acosf( x ); }
inline float asin( float x ) { return asinf( x ); }
inline float exp( float x ) { return expf( x ); }
inline float pow( float x, float y ) { return powf( x, y ); }
inline float atan( float yoverx ) { return atanf( yoverx ); }
#endif

inline float sqr( float x ) { return x * x; }
inline float log( float x, float base ) { return log( x ) / log( base ); }

inline float atan( float y, float x ) { return atan2f( y, x ); }
template <typename T> T lerp( float alpha, const T &a, const T &b ) { return a + alpha * ( b - a ); }

#undef min
inline float min ( float x, float y ) { return x < y ? x : y; }
#undef max
inline float max ( float x, float y ) { return x > y ? x : y; }

inline unsigned long urand();
float random();

inline bool outside( float value, float min, float max )
{
	return ( min < max ) ? ( value < min || value > max ) : ( value > min && value < max );
}

//inline void swap( float &a, float &b ) { float tmp = a; a = b; b = tmp; }

inline bool 
quadratic( float A, float B, float C, float *t0, float *t1 ) 
{
	// Find quadratic discriminant
	float discrim = B * B - 4.0 * A * C;
	if (discrim < 0.) return false;
	float rootDiscrim = sqrt( discrim );
	// Compute quadratic [[t]] values
	float q;
	if ( B < 0 ) 
		q = -.5f * ( B - rootDiscrim );
	else       
		q = -.5f * ( B + rootDiscrim );
	*t0 = q / A;
	*t1 = C / q;
	if (*t0 > *t1) std::swap( *t0, *t1 );
	return true;
}

template < class T, int n > T 
deCasteljau( T Q[], float t )
{
	for ( int i = 1; i < n; i++ )
		for ( int j = 0; j < n - i; j++ )
			Q[j] = lerp( t, Q[j], Q[j + 1] );

	return Q[0];
}

template < class T > T 
deCasteljau( T Q[], int order, float t )
{
	for ( int i = 1; i < order; i++ )
		for ( int j = 0; j < order - i; j++ )
			Q[j] = lerp( t, Q[j], Q[j + 1] );

	return Q[0];
}

template < class T, int n > T 
deCasteljau( const T values[], float t )
{
	T Q[n-1];

	for ( int j = 0; j < n - 1; j++ )
		Q[j] = lerp( t, values[j], values[j + 1] );
	for ( int i = 2; i < n; i++ )
		for ( int j = 0; j < n - i; j++ )
			Q[j] = lerp( t, Q[j], Q[j + 1] );

	return Q[0];
}

template < class T > inline bool
solveLinearSystem2x2( const float A[2][2],
		const T B[2], T x[2] ) 
{
	float det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
	if ( abs(det) < 1e-5 )
		return false;
	float invDet = 1.0f/det;
	x[0] = (A[1][1]*B[0] - A[0][1]*B[1]) * invDet;
	x[1] = (A[0][0]*B[1] - A[1][0]*B[0]) * invDet;
	return true;
}

#endif //__LIBMATH_H__