// Copyright (C) 2012-2014 ChaosForge Ltd
// http://chaosforge.org/
//
// This file is part of NV Libraries.
// For conditions of distribution and use, see copyright notice in nv.hh

#include "nv/core/random.hh"
#include "nv/core/time.hh"
#include <random>

using namespace nv;

random::random( random::seed_type seed /*= 0 */ )
{
	static_assert( sizeof( std::mt19937 ) < sizeof( random::m_data ), "No room for mersenne twister!" );
	new (m_data)std::mt19937( seed == 0 ? randomized_seed() : seed );
}

random::seed_type random::randomize()
{
	std::mt19937& rng = *(( std::mt19937* )m_data);
	seed_type seed = randomized_seed();
	rng.seed( seed );
	return seed;
}

void random::set_seed( random::seed_type seed /*= 0 */ )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	rng.seed( seed == 0 ? randomized_seed() : seed );
}

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

random::result_type random::rand()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	return rng();
}

sint32 random::srand( sint32 val )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_int_distribution<sint32> dist( 0, val - 1 );
	return dist( rng );
}

uint32 random::urand( uint32 val )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_int_distribution<uint32> dist( 0, val - 1 );
	return dist( rng );
}

f32 random::frand( f32 val )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist( 0, val );
	return dist( rng );
}

sint32 random::srange( sint32 min, sint32 max )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_int_distribution<sint32> dist( min, max );
	return dist( rng );
}

uint32 random::urange( uint32 min, uint32 max )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_int_distribution<uint32> dist( min, max );
	return dist( rng );
}

f32 random::frange( f32 min, f32 max )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist( min, max );
	return dist( rng );
}

uint32 random::dice( uint32 count, uint32 sides )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_int_distribution<uint32> dist( 1, sides );
	uint32 result = 0;
	while (count-- > 0)
	{
		result += dist( rng );
	};
	return result;
}

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()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist( 0, glm::pi<float>() * 2.f );
	float angle = dist( rng );
	return vec2( glm::cos( angle ), glm::sin( angle ) );
}

nv::vec3 nv::random::precise_unit_vec3()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist11( -1.0f, 1.0f );
	std::uniform_real_distribution<f32> dist02pi( 0.0f, 2*glm::pi<float>() );
	float cos_theta = dist11( rng );
	float sin_theta = glm::sqrt( 1.0f - cos_theta * cos_theta );
	float phi       = dist02pi( rng );
	return vec3( 
		sin_theta * glm::sin(phi),
		sin_theta * glm::cos(phi),
		cos_theta
		);
}

nv::vec2 nv::random::fast_disk_point()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist( 0.0f, 1.0f );
	float r1 = dist( rng );
	float r2 = dist( rng );
	if ( r1 > r2 ) std::swap( r1, r2 );
	float rf = 2*glm::pi<float>()*(r1/r2);
	return vec2( r2*glm::cos( rf ), r2*glm::sin( rf ) );
}

nv::vec2 nv::random::precise_disk_point()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> unit( 0.0f, 1.0f );
	std::uniform_real_distribution<f32> angle( 0.0f, glm::pi<float>() );
	float r = glm::sqrt( unit( rng ) );
	float rangle = angle( rng );
	return vec2( r*glm::cos( rangle ), r*glm::sin( rangle ) );
}

nv::vec3 nv::random::fast_sphere_point()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist01( 0.0f, 1.0f );
	std::uniform_real_distribution<f32> dist11( -1.0f, 1.0f );
	std::uniform_real_distribution<f32> dist02pi( 0.0f, 2*glm::pi<float>() );
	float rad     = dist01( rng );
	float pi      = glm::pi<float>();
	float phi     = glm::asin( dist11( rng ) ) + pi*.5f;
	float theta   = dist02pi( rng );
	float sin_phi = glm::sin( phi );
	return vec3( 
		rad * glm::cos(theta) * sin_phi,
		rad * glm::sin(theta) * sin_phi,
		rad * glm::cos(phi)
	);
}

nv::vec3 nv::random::precise_sphere_point()
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> dist01( 0.0f, 1.0f );
	std::uniform_real_distribution<f32> dist11( -1.0f, 1.0f );
	std::uniform_real_distribution<f32> dist02pi( 0.0f, 2*glm::pi<float>() );
	float radius = std::pow( dist01( rng ), 1.f/3.f );
	float cos_theta = dist11( rng );
	float sin_theta = glm::sqrt( 1.0f - cos_theta * cos_theta );
	float phi       = dist02pi( rng );
	return vec3( 
		radius * sin_theta * glm::sin(phi),
		radius * sin_theta * glm::cos(phi),
		radius * cos_theta
		);
}

nv::vec2 nv::random::precise_ellipse_point( const vec2& radii )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> distx( -radii.x, radii.x );
	std::uniform_real_distribution<f32> disty( -radii.y, radii.y );
	vec2 inv_radii  = 1.f / radii;
	vec2 inv_radii2 = inv_radii * inv_radii;
	for ( uint32 i = 0; i < 12; ++i )
	{
		float x = distx( rng );
		float y = disty( rng );
		if ( x * x * inv_radii2.x + y * y * inv_radii2.y <= 1.f )
		{
			return vec2( x, y );
		}
	}
	return fast_disk_point() * radii;
}

nv::vec3 nv::random::precise_ellipsoid_point( const vec3& radii )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	std::uniform_real_distribution<f32> distx( -radii.x, radii.x );
	std::uniform_real_distribution<f32> disty( -radii.y, radii.y );
	std::uniform_real_distribution<f32> distz( -radii.z, radii.z );
	vec3 inv_radii  = 1.f / radii;
	vec3 inv_radii2 = inv_radii * inv_radii;
	for ( uint32 i = 0; i < 12; ++i )
	{
		float x = distx( rng );
		float y = disty( rng );
		float z = distz( rng );
		if ( x * x * inv_radii2.x + y * y * inv_radii2.y + z * z * inv_radii2.z <= 1.f )
		{
			return vec3( x, y, z );
		}
	}
	return fast_sphere_point() * radii;
}

nv::vec2 nv::random::fast_hollow_disk_point( float iradius, float oradius )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	float idist2 = iradius * iradius;
	float odist2 = oradius * oradius;
	float rdist  = glm::sqrt( std::uniform_real_distribution<f32>( idist2, odist2 )( rng ) );
	return rdist * precise_unit_vec2();
}

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

nv::vec3 nv::random::fast_hollow_sphere_point( float iradius, float oradius )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	float idist3 = iradius * iradius * iradius;
	float odist3 = oradius * oradius * oradius;
	float rdist  = std::pow( std::uniform_real_distribution<f32>( idist3, odist3 )( rng ), 1.f/3.f );
	return rdist * precise_unit_vec3();
}

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


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

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

	float rdist     = glm::sqrt( std::uniform_real_distribution<f32>( idist2, odist2 )( rng ) );
	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 )
{
	std::mt19937& rng = *( ( std::mt19937* )m_data );
	vec3 iradii2 = iradii * iradii;
	vec3 opoint     = ellipsoid_edge( oradii );
	vec3 opoint2    = opoint * opoint;
	vec3 odir       = glm::normalize( opoint );
	float odist2    = opoint2.x + opoint2.y + opoint2.z;

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

	float odist3 = odist2 * glm::sqrt( odist2 );
	float idist3 = idist2 * glm::sqrt( idist2 );

	float rdist     = std::pow( std::uniform_real_distribution<f32>( idist3, odist3 )( rng ), 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 );
}