Index: trunk/src/core/random.cc
===================================================================
--- trunk/src/core/random.cc	(revision 442)
+++ trunk/src/core/random.cc	(revision 443)
@@ -7,26 +7,78 @@
 #include "nv/core/random.hh"
 #include "nv/core/time.hh"
-#include <random>
+#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 */ )
-{
-	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 );
+	: 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;
-}
-
-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 );
 }
 
@@ -39,60 +91,12 @@
 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 );
+	return mt_uint32();
 }
 
 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;
+	uint32 x, max = mt_full_mask - ( mt_full_mask % val );
+	while ( ( x = rand() ) >= max );
+	return x / ( max / val );
 }
 
@@ -106,7 +110,5 @@
 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 );
+	float angle = frand( glm::pi<float>() * 2.f );
 	return vec2( glm::cos( angle ), glm::sin( angle ) );
 }
@@ -114,10 +116,7 @@
 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 cos_theta = frange( -1.0f, 1.0f );
 	float sin_theta = glm::sqrt( 1.0f - cos_theta * cos_theta );
-	float phi       = dist02pi( rng );
+	float phi       = frand( 2 * glm::pi<float>() );
 	return vec3( 
 		sin_theta * glm::sin(phi),
@@ -129,9 +128,7 @@
 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 r1 = frand();
+	float r2 = frand();
+	if ( r1 > r2 ) swap( r1, r2 );
 	float rf = 2*glm::pi<float>()*(r1/r2);
 	return vec2( r2*glm::cos( rf ), r2*glm::sin( rf ) );
@@ -140,9 +137,6 @@
 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 );
+	float r = glm::sqrt( frand() );
+	float rangle = frand( glm::pi<float>() );
 	return vec2( r*glm::cos( rangle ), r*glm::sin( rangle ) );
 }
@@ -150,12 +144,8 @@
 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 rad     = frand();
 	float pi      = glm::pi<float>();
-	float phi     = glm::asin( dist11( rng ) ) + pi*.5f;
-	float theta   = dist02pi( rng );
+	float phi     = glm::asin( frange( -1.0f, 1.0f ) ) + pi*.5f;
+	float theta   = frange( 0.0f, 2 * glm::pi<float>() );
 	float sin_phi = glm::sin( phi );
 	return vec3( 
@@ -168,12 +158,8 @@
 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 radius = glm::pow( frand(), 1.f/3.f );
+	float cos_theta = frange( -1.0f, 1.0f );
 	float sin_theta = glm::sqrt( 1.0f - cos_theta * cos_theta );
-	float phi       = dist02pi( rng );
+	float phi       = frange( 0.0f, 2 * glm::pi<float>() );
 	return vec3( 
 		radius * sin_theta * glm::sin(phi),
@@ -185,16 +171,12 @@
 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 p = range( -radii, radii );
 	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 )
+		if ( p.x * p.x * inv_radii2.x + p.y * p.y * inv_radii2.y <= 1.f )
 		{
-			return vec2( x, y );
+			return p;
 		}
 	}
@@ -204,18 +186,12 @@
 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 p = range( -radii, radii );
 	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 )
+		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 vec3( x, y, z );
+			return p;
 		}
 	}
@@ -225,8 +201,7 @@
 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 ) );
+	float rdist  = glm::sqrt( frange( idist2, odist2 ) );
 	return rdist * precise_unit_vec2();
 }
@@ -239,8 +214,7 @@
 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 );
+	float rdist  = glm::pow( frange( idist3, odist3 ), 1.f/3.f );
 	return rdist * precise_unit_vec3();
 }
@@ -254,5 +228,4 @@
 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 );
@@ -264,5 +237,5 @@
 	float idist2 = ((iradii2.x * iradii2.y) / low ) * odist2;
 
-	float rdist     = glm::sqrt( std::uniform_real_distribution<f32>( idist2, odist2 )( rng ) );
+	float rdist     = glm::sqrt( frange( idist2, odist2 ) );
 	return odir * rdist;	
 }
@@ -275,5 +248,4 @@
 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 );
@@ -291,5 +263,5 @@
 	float idist3 = idist2 * glm::sqrt( idist2 );
 
-	float rdist     = std::pow( std::uniform_real_distribution<f32>( idist3, odist3 )( rng ), 1.f/3.f );
+	float rdist     = glm::pow( frange( idist3, odist3 ), 1.f/3.f );
 	return odir * rdist;	
 }