Index: trunk/src/random.cc
===================================================================
--- trunk/src/random.cc	(revision 307)
+++ trunk/src/random.cc	(revision 308)
@@ -90,2 +90,187 @@
 	return narrow_cast< seed_type >( get_ticks() );
 }
+
+nv::vec2 nv::random::precise_unit_vec2()
+{
+	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::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::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::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::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::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::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::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 )
+{
+	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 )
+{
+	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 )
+{
+	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 )
+{
+	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 );
+}
+