// Copyright (C) 2012-2014 ChaosForge Ltd
// http://chaosforge.org/
//
// This file is part of Nova libraries. 
// For conditions of distribution and use, see copying.txt file in root folder.

#include "nv/core/random.hh"
#include "nv/core/time.hh"
#include "nv/stl/utility/common.hh"

using namespace nv;

static const uint32 mt_upper_mask = 0x80000000UL;
static const uint32 mt_lower_mask = 0x7FFFFFFFUL;
static const uint32 mt_full_mask  = 0xFFFFFFFFUL;
static const uint32 mt_matrix_a   = 0x9908B0DFUL;

#define NV_MT_MIXBITS(u, v) ( ( (u) & mt_upper_mask) | ( (v) & mt_lower_mask) )
#define NV_MT_TWIST(u, v)  ( (NV_MT_MIXBITS(u, v) >> 1) ^ ( (v) & 1UL ? mt_matrix_a : 0UL) )

void random::mt_init( uint32 seed )
{
	m_state[0] = static_cast<uint32_t>( seed & mt_full_mask );
	for ( int i = 1; i < mersenne_n; i++ )
	{
		m_state[i]  = ( 1812433253UL * ( m_state[i - 1] ^ ( m_state[i - 1] >> 30 ) ) + i );
		m_state[i] &= mt_full_mask;
	}

	m_remaining = 0;
	m_next      = nullptr;
	m_seeded    = 1;
}


void random::mt_update()
{
	uint32 *p = m_state;

	for ( int count = ( mersenne_n - mersenne_m + 1 ); --count; p++ )
		*p = p[mersenne_m] ^ NV_MT_TWIST( p[0], p[1] );

	for ( int count = mersenne_m; --count; p++ )
		*p = p[mersenne_m - mersenne_n] ^ NV_MT_TWIST( p[0], p[1] );

	*p = p[mersenne_m - mersenne_n] ^ NV_MT_TWIST( p[0], m_state[0] );

	m_remaining = mersenne_n;
	m_next      = m_state;
}


uint32 random::mt_uint32()
{
	uint32 r;

	if ( !m_remaining ) mt_update();

	r = *m_next++;
	m_remaining--;

	r ^= ( r >> 11 );
	r ^= ( r << 7 ) & 0x9D2C5680UL;
	r ^= ( r << 15 ) & 0xEFC60000UL;
	r ^= ( r >> 18 );

	r &= mt_full_mask;

	return r;
}

random::random( random::seed_type seed /*= 0 */ )
	: m_remaining( 0 ), m_next( nullptr ), m_seeded( 0 )
{
	mt_init( seed == 0 ? randomized_seed() : seed );
}

random::seed_type random::set_seed( random::seed_type seed /*= 0 */ )
{
	if ( seed == 0 ) seed = randomized_seed();
	mt_init( seed );
	return seed;
}

nv::random& random::get()
{
	static random default_rng;
	return default_rng;
}

random::result_type random::rand()
{
	return mt_uint32();
}

uint32 random::urand( uint32 val )
{
	uint32 x, max = mt_full_mask - ( mt_full_mask % val );
	while ( ( x = rand() ) >= max );
	return x / ( max / val );
}

random::seed_type random::randomized_seed()
{
	// TODO: this seems off, as it might often seed the same, use general time 
	// instead
	return narrow_cast< seed_type >( get_ticks() );
}

nv::vec2 nv::random::precise_unit_vec2()
{
	f32 angle = frand( math::pi<f32>() * 2.f );
	return vec2( math::cos( angle ), math::sin( angle ) );
}

nv::vec3 nv::random::precise_unit_vec3()
{
	f32 cos_theta = frange( -1.0f, 1.0f );
	f32 sin_theta = math::sqrt( 1.0f - cos_theta * cos_theta );
	f32 phi       = frand( 2 * math::pi<f32>() );
	return vec3( 
		sin_theta * math::sin(phi),
		sin_theta * math::cos(phi),
		cos_theta
		);
}

nv::vec2 nv::random::fast_disk_point()
{
	f32 r1 = frand();
	f32 r2 = frand();
	if ( r1 > r2 ) swap( r1, r2 );
	f32 rf = 2* math::pi<f32>()*(r1/r2);
	return vec2( r2*math::cos( rf ), r2*math::sin( rf ) );
}

nv::vec2 nv::random::precise_disk_point()
{
	f32 r = math::sqrt( frand() );
	f32 rangle = frand( math::pi<f32>() );
	return vec2( r*math::cos( rangle ), r*math::sin( rangle ) );
}

nv::vec3 nv::random::fast_sphere_point()
{
	f32 rad     = frand();
	f32 pi      = math::pi<f32>();
	f32 phi     = math::asin( frange( -1.0f, 1.0f ) ) + pi*.5f;
	f32 theta   = frange( 0.0f, 2 * math::pi<f32>() );
	f32 sin_phi = math::sin( phi );
	return vec3( 
		rad * math::cos(theta) * sin_phi,
		rad * math::sin(theta) * sin_phi,
		rad * math::cos(phi)
	);
}

nv::vec3 nv::random::precise_sphere_point()
{
	f32 radius = math::pow( frand(), 1.f/3.f );
	f32 cos_theta = frange( -1.0f, 1.0f );
	f32 sin_theta = math::sqrt( 1.0f - cos_theta * cos_theta );
	f32 phi       = frange( 0.0f, 2 * math::pi<f32>() );
	return vec3( 
		radius * sin_theta * math::sin(phi),
		radius * sin_theta * math::cos(phi),
		radius * cos_theta
		);
}

nv::vec2 nv::random::precise_ellipse_point( const vec2& radii )
{
	vec2 p = range( -radii, radii );
	vec2 inv_radii  = 1.f / radii;
	vec2 inv_radii2 = inv_radii * inv_radii;
	for ( uint32 i = 0; i < 12; ++i )
	{
		if ( p.x * p.x * inv_radii2.x + p.y * p.y * inv_radii2.y <= 1.f )
		{
			return p;
		}
	}
	return fast_disk_point() * radii;
}

nv::vec3 nv::random::precise_ellipsoid_point( const vec3& radii )
{
	vec3 p = range( -radii, radii );
	vec3 inv_radii  = 1.f / radii;
	vec3 inv_radii2 = inv_radii * inv_radii;
	for ( uint32 i = 0; i < 12; ++i )
	{
		if ( p.x * p.x * inv_radii2.x + p.y * p.y * inv_radii2.y + p.z * p.z * inv_radii2.z <= 1.f )
		{
			return p;
		}
	}
	return fast_sphere_point() * radii;
}

nv::vec2 nv::random::fast_hollow_disk_point( f32 iradius, f32 oradius )
{
	f32 idist2 = iradius * iradius;
	f32 odist2 = oradius * oradius;
	f32 rdist  = math::sqrt( frange( idist2, odist2 ) );
	return rdist * precise_unit_vec2();
}

nv::vec2 nv::random::precise_hollow_disk_point( f32 iradius, f32 oradius )
{
	return fast_hollow_disk_point( iradius, oradius );
}

nv::vec3 nv::random::fast_hollow_sphere_point( f32 iradius, f32 oradius )
{
	f32 idist3 = iradius * iradius * iradius;
	f32 odist3 = oradius * oradius * oradius;
	f32 rdist  = math::pow( frange( idist3, odist3 ), 1.f/3.f );
	return rdist * precise_unit_vec3();
}

nv::vec3 nv::random::precise_hollow_sphere_point( f32 iradius, f32 oradius )
{
	return fast_hollow_sphere_point( iradius, oradius );
}


nv::vec2 nv::random::fast_hollow_ellipse_point( const vec2& iradii, const vec2& oradii )
{
	vec2 iradii2 = iradii * iradii;
	vec2 opoint     = ellipse_edge( oradii );
	vec2 opoint2    = opoint * opoint;
	vec2 odir       = glm::normalize( opoint );
	f32 odist2    = opoint2.x + opoint2.y;

	f32 low    = iradii2.y * opoint2.x + iradii2.x * opoint2.y;
	f32 idist2 = ((iradii2.x * iradii2.y) / low ) * odist2;

	f32 rdist     = math::sqrt( frange( idist2, odist2 ) );
	return odir * rdist;	
}

nv::vec2 nv::random::precise_hollow_ellipse_point( const vec2& iradii, const vec2& oradii )
{
	return fast_hollow_ellipse_point( iradii, oradii );
}

nv::vec3 nv::random::fast_hollow_ellipsoid_point( const vec3& iradii, const vec3& oradii )
{
	vec3 iradii2 = iradii * iradii;
	vec3 opoint     = ellipsoid_edge( oradii );
	vec3 opoint2    = opoint * opoint;
	vec3 odir       = glm::normalize( opoint );
	f32 odist2    = opoint2.x + opoint2.y + opoint2.z;

	f32 low    = 
		iradii2.y * iradii2.z * opoint2.x + 
		iradii2.x * iradii2.z * opoint2.y +
		iradii2.x * iradii2.y * opoint2.z;
	f32 idist2 = ((iradii2.x * iradii2.y * iradii2.z) / low ) * odist2;

	f32 odist3 = odist2 * math::sqrt( odist2 );
	f32 idist3 = idist2 * math::sqrt( idist2 );

	f32 rdist     = math::pow( frange( idist3, odist3 ), 1.f/3.f );
	return odir * rdist;	
}

nv::vec3 nv::random::precise_hollow_ellipsoid_point( const vec3& iradii, const vec3& oradii )
{
	return fast_hollow_ellipsoid_point( iradii, oradii );
}