Ignore:
Timestamp:
07/24/15 15:16:09 (10 years ago)
Author:
epyon
Message:
  • random - mersenne twister implementation
  • nv is officially STD free :P
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/core/random.cc

    r402 r443  
    77#include "nv/core/random.hh"
    88#include "nv/core/time.hh"
    9 #include <random>
     9#include "nv/stl/utility/common.hh"
    1010
    1111using namespace nv;
    1212
     13static const uint32 mt_upper_mask = 0x80000000UL;
     14static const uint32 mt_lower_mask = 0x7FFFFFFFUL;
     15static const uint32 mt_full_mask  = 0xFFFFFFFFUL;
     16static const uint32 mt_matrix_a   = 0x9908B0DFUL;
     17
     18#define NV_MT_MIXBITS(u, v) ( ( (u) & mt_upper_mask) | ( (v) & mt_lower_mask) )
     19#define NV_MT_TWIST(u, v)  ( (NV_MT_MIXBITS(u, v) >> 1) ^ ( (v) & 1UL ? mt_matrix_a : 0UL) )
     20
     21void random::mt_init( uint32 seed )
     22{
     23        m_state[0] = static_cast<uint32_t>( seed & mt_full_mask );
     24        for ( int i = 1; i < mersenne_n; i++ )
     25        {
     26                m_state[i]  = ( 1812433253UL * ( m_state[i - 1] ^ ( m_state[i - 1] >> 30 ) ) + i );
     27                m_state[i] &= mt_full_mask;
     28        }
     29
     30        m_remaining = 0;
     31        m_next      = nullptr;
     32        m_seeded    = 1;
     33}
     34
     35
     36void random::mt_update()
     37{
     38        uint32 *p = m_state;
     39
     40        for ( int count = ( mersenne_n - mersenne_m + 1 ); --count; p++ )
     41                *p = p[mersenne_m] ^ NV_MT_TWIST( p[0], p[1] );
     42
     43        for ( int count = mersenne_m; --count; p++ )
     44                *p = p[mersenne_m - mersenne_n] ^ NV_MT_TWIST( p[0], p[1] );
     45
     46        *p = p[mersenne_m - mersenne_n] ^ NV_MT_TWIST( p[0], m_state[0] );
     47
     48        m_remaining = mersenne_n;
     49        m_next      = m_state;
     50}
     51
     52
     53uint32 random::mt_uint32()
     54{
     55        uint32 r;
     56
     57        if ( !m_remaining ) mt_update();
     58
     59        r = *m_next++;
     60        m_remaining--;
     61
     62        r ^= ( r >> 11 );
     63        r ^= ( r << 7 ) & 0x9D2C5680UL;
     64        r ^= ( r << 15 ) & 0xEFC60000UL;
     65        r ^= ( r >> 18 );
     66
     67        r &= mt_full_mask;
     68
     69        return r;
     70}
     71
    1372random::random( random::seed_type seed /*= 0 */ )
    14 {
    15         static_assert( sizeof( std::mt19937 ) < sizeof( random::m_data ), "No room for mersenne twister!" );
    16         new (m_data)std::mt19937( seed == 0 ? randomized_seed() : seed );
    17 }
    18 
    19 random::seed_type random::randomize()
    20 {
    21         std::mt19937& rng = *( reinterpret_cast< std::mt19937* >( m_data ) );
    22         seed_type seed = randomized_seed();
    23         rng.seed( seed );
     73        : m_remaining( 0 ), m_next( nullptr ), m_seeded( 0 )
     74{
     75        mt_init( seed == 0 ? randomized_seed() : seed );
     76}
     77
     78random::seed_type random::set_seed( random::seed_type seed /*= 0 */ )
     79{
     80        if ( seed == 0 ) seed = randomized_seed();
     81        mt_init( seed );
    2482        return seed;
    25 }
    26 
    27 void random::set_seed( random::seed_type seed /*= 0 */ )
    28 {
    29         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    30         rng.seed( seed == 0 ? randomized_seed() : seed );
    3183}
    3284
    … …  
    3991random::result_type random::rand()
    4092{
    41         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    42         return rng();
    43 }
    44 
    45 sint32 random::srand( sint32 val )
    46 {
    47         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    48         std::uniform_int_distribution<sint32> dist( 0, val - 1 );
    49         return dist( rng );
     93        return mt_uint32();
    5094}
    5195
    5296uint32 random::urand( uint32 val )
    5397{
    54         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    55         std::uniform_int_distribution<uint32> dist( 0, val - 1 );
    56         return dist( rng );
    57 }
    58 
    59 f32 random::frand( f32 val )
    60 {
    61         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    62         std::uniform_real_distribution<f32> dist( 0, val );
    63         return dist( rng );
    64 }
    65 
    66 sint32 random::srange( sint32 min, sint32 max )
    67 {
    68         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    69         std::uniform_int_distribution<sint32> dist( min, max );
    70         return dist( rng );
    71 }
    72 
    73 uint32 random::urange( uint32 min, uint32 max )
    74 {
    75         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    76         std::uniform_int_distribution<uint32> dist( min, max );
    77         return dist( rng );
    78 }
    79 
    80 f32 random::frange( f32 min, f32 max )
    81 {
    82         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    83         std::uniform_real_distribution<f32> dist( min, max );
    84         return dist( rng );
    85 }
    86 
    87 uint32 random::dice( uint32 count, uint32 sides )
    88 {
    89         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    90         std::uniform_int_distribution<uint32> dist( 1, sides );
    91         uint32 result = 0;
    92         while (count-- > 0)
    93         {
    94                 result += dist( rng );
    95         };
    96         return result;
     98        uint32 x, max = mt_full_mask - ( mt_full_mask % val );
     99        while ( ( x = rand() ) >= max );
     100        return x / ( max / val );
    97101}
    98102
    … …  
    106110nv::vec2 nv::random::precise_unit_vec2()
    107111{
    108         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    109         std::uniform_real_distribution<f32> dist( 0, glm::pi<float>() * 2.f );
    110         float angle = dist( rng );
     112        float angle = frand( glm::pi<float>() * 2.f );
    111113        return vec2( glm::cos( angle ), glm::sin( angle ) );
    112114}
    … …  
    114116nv::vec3 nv::random::precise_unit_vec3()
    115117{
    116         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    117         std::uniform_real_distribution<f32> dist11( -1.0f, 1.0f );
    118         std::uniform_real_distribution<f32> dist02pi( 0.0f, 2*glm::pi<float>() );
    119         float cos_theta = dist11( rng );
     118        float cos_theta = frange( -1.0f, 1.0f );
    120119        float sin_theta = glm::sqrt( 1.0f - cos_theta * cos_theta );
    121         float phi       = dist02pi( rng );
     120        float phi       = frand( 2 * glm::pi<float>() );
    122121        return vec3(
    123122                sin_theta * glm::sin(phi),
    … …  
    129128nv::vec2 nv::random::fast_disk_point()
    130129{
    131         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    132         std::uniform_real_distribution<f32> dist( 0.0f, 1.0f );
    133         float r1 = dist( rng );
    134         float r2 = dist( rng );
    135         if ( r1 > r2 ) std::swap( r1, r2 );
     130        float r1 = frand();
     131        float r2 = frand();
     132        if ( r1 > r2 ) swap( r1, r2 );
    136133        float rf = 2*glm::pi<float>()*(r1/r2);
    137134        return vec2( r2*glm::cos( rf ), r2*glm::sin( rf ) );
    … …  
    140137nv::vec2 nv::random::precise_disk_point()
    141138{
    142         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    143         std::uniform_real_distribution<f32> unit( 0.0f, 1.0f );
    144         std::uniform_real_distribution<f32> angle( 0.0f, glm::pi<float>() );
    145         float r = glm::sqrt( unit( rng ) );
    146         float rangle = angle( rng );
     139        float r = glm::sqrt( frand() );
     140        float rangle = frand( glm::pi<float>() );
    147141        return vec2( r*glm::cos( rangle ), r*glm::sin( rangle ) );
    148142}
    … …  
    150144nv::vec3 nv::random::fast_sphere_point()
    151145{
    152         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    153         std::uniform_real_distribution<f32> dist01( 0.0f, 1.0f );
    154         std::uniform_real_distribution<f32> dist11( -1.0f, 1.0f );
    155         std::uniform_real_distribution<f32> dist02pi( 0.0f, 2*glm::pi<float>() );
    156         float rad     = dist01( rng );
     146        float rad     = frand();
    157147        float pi      = glm::pi<float>();
    158         float phi     = glm::asin( dist11( rng ) ) + pi*.5f;
    159         float theta   = dist02pi( rng );
     148        float phi     = glm::asin( frange( -1.0f, 1.0f ) ) + pi*.5f;
     149        float theta   = frange( 0.0f, 2 * glm::pi<float>() );
    160150        float sin_phi = glm::sin( phi );
    161151        return vec3(
    … …  
    168158nv::vec3 nv::random::precise_sphere_point()
    169159{
    170         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    171         std::uniform_real_distribution<f32> dist01( 0.0f, 1.0f );
    172         std::uniform_real_distribution<f32> dist11( -1.0f, 1.0f );
    173         std::uniform_real_distribution<f32> dist02pi( 0.0f, 2*glm::pi<float>() );
    174         float radius = std::pow( dist01( rng ), 1.f/3.f );
    175         float cos_theta = dist11( rng );
     160        float radius = glm::pow( frand(), 1.f/3.f );
     161        float cos_theta = frange( -1.0f, 1.0f );
    176162        float sin_theta = glm::sqrt( 1.0f - cos_theta * cos_theta );
    177         float phi       = dist02pi( rng );
     163        float phi       = frange( 0.0f, 2 * glm::pi<float>() );
    178164        return vec3(
    179165                radius * sin_theta * glm::sin(phi),
    … …  
    185171nv::vec2 nv::random::precise_ellipse_point( const vec2& radii )
    186172{
    187         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    188         std::uniform_real_distribution<f32> distx( -radii.x, radii.x );
    189         std::uniform_real_distribution<f32> disty( -radii.y, radii.y );
     173        vec2 p = range( -radii, radii );
    190174        vec2 inv_radii  = 1.f / radii;
    191175        vec2 inv_radii2 = inv_radii * inv_radii;
    192176        for ( uint32 i = 0; i < 12; ++i )
    193177        {
    194                 float x = distx( rng );
    195                 float y = disty( rng );
    196                 if ( x * x * inv_radii2.x + y * y * inv_radii2.y <= 1.f )
     178                if ( p.x * p.x * inv_radii2.x + p.y * p.y * inv_radii2.y <= 1.f )
    197179                {
    198                         return vec2( x, y );
     180                        return p;
    199181                }
    200182        }
    … …  
    204186nv::vec3 nv::random::precise_ellipsoid_point( const vec3& radii )
    205187{
    206         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    207         std::uniform_real_distribution<f32> distx( -radii.x, radii.x );
    208         std::uniform_real_distribution<f32> disty( -radii.y, radii.y );
    209         std::uniform_real_distribution<f32> distz( -radii.z, radii.z );
     188        vec3 p = range( -radii, radii );
    210189        vec3 inv_radii  = 1.f / radii;
    211190        vec3 inv_radii2 = inv_radii * inv_radii;
    212191        for ( uint32 i = 0; i < 12; ++i )
    213192        {
    214                 float x = distx( rng );
    215                 float y = disty( rng );
    216                 float z = distz( rng );
    217                 if ( x * x * inv_radii2.x + y * y * inv_radii2.y + z * z * inv_radii2.z <= 1.f )
     193                if ( p.x * p.x * inv_radii2.x + p.y * p.y * inv_radii2.y + p.z * p.z * inv_radii2.z <= 1.f )
    218194                {
    219                         return vec3( x, y, z );
     195                        return p;
    220196                }
    221197        }
    … …  
    225201nv::vec2 nv::random::fast_hollow_disk_point( float iradius, float oradius )
    226202{
    227         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    228203        float idist2 = iradius * iradius;
    229204        float odist2 = oradius * oradius;
    230         float rdist  = glm::sqrt( std::uniform_real_distribution<f32>( idist2, odist2 )( rng ) );
     205        float rdist  = glm::sqrt( frange( idist2, odist2 ) );
    231206        return rdist * precise_unit_vec2();
    232207}
    … …  
    239214nv::vec3 nv::random::fast_hollow_sphere_point( float iradius, float oradius )
    240215{
    241         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    242216        float idist3 = iradius * iradius * iradius;
    243217        float odist3 = oradius * oradius * oradius;
    244         float rdist  = std::pow( std::uniform_real_distribution<f32>( idist3, odist3 )( rng ), 1.f/3.f );
     218        float rdist  = glm::pow( frange( idist3, odist3 ), 1.f/3.f );
    245219        return rdist * precise_unit_vec3();
    246220}
    … …  
    254228nv::vec2 nv::random::fast_hollow_ellipse_point( const vec2& iradii, const vec2& oradii )
    255229{
    256         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    257230        vec2 iradii2 = iradii * iradii;
    258231        vec2 opoint     = ellipse_edge( oradii );
    … …  
    264237        float idist2 = ((iradii2.x * iradii2.y) / low ) * odist2;
    265238
    266         float rdist     = glm::sqrt( std::uniform_real_distribution<f32>( idist2, odist2 )( rng ) );
     239        float rdist     = glm::sqrt( frange( idist2, odist2 ) );
    267240        return odir * rdist;   
    268241}
    … …  
    275248nv::vec3 nv::random::fast_hollow_ellipsoid_point( const vec3& iradii, const vec3& oradii )
    276249{
    277         std::mt19937& rng = *( reinterpret_cast<std::mt19937*>( m_data ) );
    278250        vec3 iradii2 = iradii * iradii;
    279251        vec3 opoint     = ellipsoid_edge( oradii );
    … …  
    291263        float idist3 = idist2 * glm::sqrt( idist2 );
    292264
    293         float rdist     = std::pow( std::uniform_real_distribution<f32>( idist3, odist3 )( rng ), 1.f/3.f );
     265        float rdist     = glm::pow( frange( idist3, odist3 ), 1.f/3.f );
    294266        return odir * rdist;   
    295267}
Note: See TracChangeset for help on using the changeset viewer.