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

// 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>( seed & mt_full_mask );
        for ( uint32 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_next( nullptr ), m_remaining( 0 ), 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( cos( angle ), sin( angle ) );
}

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

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

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

nv::vec3 nv::random::precise_sphere_point()
{
        f32 radius = pow( frand(), 1.f/3.f );
        f32 cos_theta = frange( -1.0f, 1.0f );
        f32 sin_theta = sqrt( 1.0f - cos_theta * cos_theta );
        f32 phi       = frange( 0.0f, 2 * math::pi<f32>() );
        return vec3( 
                radius * sin_theta * sin(phi),
                radius * sin_theta * 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  = 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  = 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       = math::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     = 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       = math::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 * sqrt( odist2 );
        f32 idist3 = idist2 * sqrt( idist2 );

        f32 rdist     = 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 );
}