+static void eval_bezier3( v3f p0, v3f p1, v3f p2, f32 t, v3f p )
+{
+ f32 u = 1.0f-t;
+
+ v3_muls( p0, u*u, p );
+ v3_muladds( p, p1, 2.0f*u*t, p );
+ v3_muladds( p, p2, t*t, p );
+}
+
+/*
+ * -----------------------------------------------------------------------------
+ * Section 5.f Volumes
+ * -----------------------------------------------------------------------------
+ */
+
+static float vg_sphere_volume( float radius ){
+ float r3 = radius*radius*radius;
+ return (4.0f/3.0f) * VG_PIf * r3;
+}
+
+/*
+ * -----------------------------------------------------------------------------
+ * Section 6.a PSRNG and some distributions
+ * -----------------------------------------------------------------------------
+ */
+
+/* An implementation of the MT19937 Algorithm for the Mersenne Twister
+ * by Evan Sultanik. Based upon the pseudocode in: M. Matsumoto and
+ * T. Nishimura, "Mersenne Twister: A 623-dimensionally
+ * equidistributed uniform pseudorandom number generator," ACM
+ * Transactions on Modeling and Computer Simulation Vol. 8, No. 1,
+ * January pp.3-30 1998.
+ *
+ * http://www.sultanik.com/Mersenne_twister
+ * https://github.com/ESultanik/mtwister/blob/master/mtwister.c
+ */
+
+#define MT_UPPER_MASK 0x80000000
+#define MT_LOWER_MASK 0x7fffffff
+#define MT_TEMPERING_MASK_B 0x9d2c5680
+#define MT_TEMPERING_MASK_C 0xefc60000
+
+#define MT_STATE_VECTOR_LENGTH 624
+
+/* changes to STATE_VECTOR_LENGTH also require changes to this */
+#define MT_STATE_VECTOR_M 397
+
+struct {
+ u32 mt[MT_STATE_VECTOR_LENGTH];
+ i32 index;
+}
+static vg_rand;
+
+static void vg_rand_seed( unsigned long seed )
+{
+ /* set initial seeds to mt[STATE_VECTOR_LENGTH] using the generator
+ * from Line 25 of Table 1 in: Donald Knuth, "The Art of Computer
+ * Programming," Vol. 2 (2nd Ed.) pp.102.
+ */
+ vg_rand.mt[0] = seed & 0xffffffff;
+ for( vg_rand.index=1; vg_rand.index<MT_STATE_VECTOR_LENGTH; vg_rand.index++){
+ vg_rand.mt[vg_rand.index] =
+ (6069 * vg_rand.mt[vg_rand.index-1]) & 0xffffffff;
+ }
+}
+
+/*
+ * Generates a pseudo-randomly generated long.
+ */
+static u32 vg_randu32(void)
+{
+ u32 y;
+ /* mag[x] = x * 0x9908b0df for x = 0,1 */
+ static u32 mag[2] = {0x0, 0x9908b0df};
+ if( vg_rand.index >= MT_STATE_VECTOR_LENGTH || vg_rand.index < 0 ){
+ /* generate STATE_VECTOR_LENGTH words at a time */
+ int kk;
+ if( vg_rand.index >= MT_STATE_VECTOR_LENGTH+1 || vg_rand.index < 0 ){
+ vg_rand_seed( 4357 );
+ }
+ for( kk=0; kk<MT_STATE_VECTOR_LENGTH-MT_STATE_VECTOR_M; kk++ ){
+ y = (vg_rand.mt[kk] & MT_UPPER_MASK) |
+ (vg_rand.mt[kk+1] & MT_LOWER_MASK);
+ vg_rand.mt[kk] = vg_rand.mt[kk+MT_STATE_VECTOR_M] ^
+ (y >> 1) ^ mag[y & 0x1];
+ }
+ for( ; kk<MT_STATE_VECTOR_LENGTH-1; kk++ ){
+ y = (vg_rand.mt[kk] & MT_UPPER_MASK) |
+ (vg_rand.mt[kk+1] & MT_LOWER_MASK);
+ vg_rand.mt[kk] =
+ vg_rand.mt[ kk+(MT_STATE_VECTOR_M-MT_STATE_VECTOR_LENGTH)] ^
+ (y >> 1) ^ mag[y & 0x1];
+ }
+ y = (vg_rand.mt[MT_STATE_VECTOR_LENGTH-1] & MT_UPPER_MASK) |
+ (vg_rand.mt[0] & MT_LOWER_MASK);
+ vg_rand.mt[MT_STATE_VECTOR_LENGTH-1] =
+ vg_rand.mt[MT_STATE_VECTOR_M-1] ^ (y >> 1) ^ mag[y & 0x1];
+ vg_rand.index = 0;
+ }
+ y = vg_rand.mt[vg_rand.index++];
+ y ^= (y >> 11);
+ y ^= (y << 7) & MT_TEMPERING_MASK_B;
+ y ^= (y << 15) & MT_TEMPERING_MASK_C;
+ y ^= (y >> 18);
+ return y;
+}
+
+/*
+ * Generates a pseudo-randomly generated f64 in the range [0..1].
+ */
+static inline f64 vg_randf64(void)
+{
+ return (f64)vg_randu32()/(f64)0xffffffff;
+}
+
+static inline f64 vg_randf64_range( f64 min, f64 max )
+{
+ return vg_lerp( min, max, (f64)vg_randf64() );
+}
+
+static inline void vg_rand_dir( v3f dir )
+{
+ dir[0] = vg_randf64();
+ dir[1] = vg_randf64();
+ dir[2] = vg_randf64();
+
+ v3_muls( dir, 2.0f, dir );
+ v3_sub( dir, (v3f){1.0f,1.0f,1.0f}, dir );
+
+ v3_normalize( dir );
+}
+
+static inline void vg_rand_sphere( v3f co )
+{
+ vg_rand_dir(co);
+ v3_muls( co, cbrtf( vg_randf64() ), co );
+}
+