source: trunk/src/core/random.cc @ 437

// 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 <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 = *( reinterpret_cast< 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 = *( reinterpret_cast<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 = *( reinterpret_cast<std::mt19937*>( m_data ) );
        return rng();
}

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

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

f32 random::frand( f32 val )
{
        std::mt19937& rng = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 = *( reinterpret_cast<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 );
}