Add gen := context.random_generator parameters to rand calls

core/math/rand/distributions.odin

@@ -8,12 +8,12 @@ float32_uniform :: float32_range
 // Triangular Distribution
 // See: http://wikipedia.org/wiki/Triangular_distribution
 @(require_results)
-float64_triangular :: proc(lo, hi: f64, mode: Maybe(f64)) -> f64 {
+float64_triangular :: proc(lo, hi: f64, mode: Maybe(f64), gen := context.random_generator) -> f64 {
 	if hi-lo == 0 {
 		return lo
 	}
 	lo, hi := lo, hi
-	u := float64()
+	u := float64(gen)
 	c := f64(0.5) if mode == nil else clamp((mode.?-lo) / (hi-lo), 0, 1)
 	if u > c {
 		u = 1-u
@@ -26,12 +26,12 @@ float64_triangular :: proc(lo, hi: f64, mode: Maybe(f64)) -> f64 {
 // Triangular Distribution
 // See: http://wikipedia.org/wiki/Triangular_distribution
 @(require_results)
-float32_triangular :: proc(lo, hi: f32, mode: Maybe(f32)) -> f32 {
+float32_triangular :: proc(lo, hi: f32, mode: Maybe(f32), gen := context.random_generator) -> f32 {
 	if hi-lo == 0 {
 		return lo
 	}
 	lo, hi := lo, hi
-	u := float32()
+	u := float32(gen)
 	c := f32(0.5) if mode == nil else clamp((mode.?-lo) / (hi-lo), 0, 1)
 	if u > c {
 		u = 1-u
@@ -44,25 +44,25 @@ float32_triangular :: proc(lo, hi: f32, mode: Maybe(f32)) -> f32 {
 
 // Normal/Gaussian Distribution
 @(require_results)
-float64_normal :: proc(mean, stddev: f64) -> f64 {
+float64_normal :: proc(mean, stddev: f64, gen := context.random_generator) -> f64 {
-	return norm_float64() * stddev + mean
+	return norm_float64(gen) * stddev + mean
 }
 // Normal/Gaussian Distribution
 @(require_results)
-float32_normal :: proc(mean, stddev: f32) -> f32 {
+float32_normal :: proc(mean, stddev: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_normal(f64(mean), f64(stddev)))
+	return f32(float64_normal(f64(mean), f64(stddev), gen))
 }
 
 
 // Log Normal Distribution
 @(require_results)
-float64_log_normal :: proc(mean, stddev: f64) -> f64 {
+float64_log_normal :: proc(mean, stddev: f64, gen := context.random_generator) -> f64 {
-	return math.exp(float64_normal(mean, stddev))
+	return math.exp(float64_normal(mean, stddev, gen))
 }
 // Log Normal Distribution
 @(require_results)
-float32_log_normal :: proc(mean, stddev: f32) -> f32 {
+float32_log_normal :: proc(mean, stddev: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_log_normal(f64(mean), f64(stddev)))
+	return f32(float64_log_normal(f64(mean), f64(stddev), gen))
 }
 
 
@@ -72,8 +72,8 @@ float32_log_normal :: proc(mean, stddev: f32) -> f32 {
 //     0 to positive infinity if lambda >  0
 //     negative infinity to 0 if lambda <= 0
 @(require_results)
-float64_exponential :: proc(lambda: f64) -> f64 {
+float64_exponential :: proc(lambda: f64, gen := context.random_generator) -> f64 {
-	return - math.ln(1 - float64()) / lambda
+	return - math.ln(1 - float64(gen)) / lambda
 }
 // Exponential Distribution
 // `lambda` is 1.0/(desired mean). It should be non-zero.
@@ -81,8 +81,8 @@ float64_exponential :: proc(lambda: f64) -> f64 {
 //     0 to positive infinity if lambda >  0
 //     negative infinity to 0 if lambda <= 0
 @(require_results)
-float32_exponential :: proc(lambda: f32) -> f32 {
+float32_exponential :: proc(lambda: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_exponential(f64(lambda)))
+	return f32(float64_exponential(f64(lambda), gen))
 }
 
 
@@ -96,7 +96,7 @@ float32_exponential :: proc(lambda: f32) -> f32 {
 //
 // mean is alpha*beta, variance is math.pow(alpha*beta, 2)
 @(require_results)
-float64_gamma :: proc(alpha, beta: f64) -> f64 {
+float64_gamma :: proc(alpha, beta: f64, gen := context.random_generator) -> f64 {
 	if alpha <= 0 || beta <= 0 {
 		panic(#procedure + ": alpha and beta must be > 0.0")
 	}
@@ -112,11 +112,11 @@ float64_gamma :: proc(alpha, beta: f64) -> f64 {
 		bbb := alpha - LOG4
 		ccc := alpha + ainv
 		for {
-			u1 := float64()
+			u1 := float64(gen)
 			if !(1e-7 < u1 && u1 < 0.9999999) {
 				continue
 			}
-			u2 := 1 - float64()
+			u2 := 1 - float64(gen)
 			v := math.ln(u1 / (1 - u1)) / ainv
 			x := alpha * math.exp(v)
 			z := u1 * u1 * u2
@@ -127,12 +127,12 @@ float64_gamma :: proc(alpha, beta: f64) -> f64 {
 		}
 	case alpha == 1:
 		// float64_exponential(1/beta)
-		return -math.ln(1 - float64()) * beta
+		return -math.ln(1 - float64(gen)) * beta
 	case:
 		// ALGORITHM GS of Statistical Computing - Kennedy & Gentle
 		x: f64
 		for {
-			u := float64()
+			u := float64(gen)
 			b := (math.e + alpha) / math.e
 			p := b * u
 			if p <= 1 {
@@ -140,7 +140,7 @@ float64_gamma :: proc(alpha, beta: f64) -> f64 {
 			} else {
 				x = -math.ln((b - p) / alpha)
 			}
-			u1 := float64()
+			u1 := float64(gen)
 			if p > 1 {
 				if u1 <= math.pow(x, alpha-1) {
 					break
@@ -162,8 +162,8 @@ float64_gamma :: proc(alpha, beta: f64) -> f64 {
 //
 // mean is alpha*beta, variance is math.pow(alpha*beta, 2)
 @(require_results)
-float32_gamma :: proc(alpha, beta: f32) -> f32 {
+float32_gamma :: proc(alpha, beta: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_gamma(f64(alpha), f64(beta)))
+	return f32(float64_gamma(f64(alpha), f64(beta), gen))
 }
 
 
@@ -173,14 +173,14 @@ float32_gamma :: proc(alpha, beta: f32) -> f32 {
 //
 // Return values range between 0 and 1
 @(require_results)
-float64_beta :: proc(alpha, beta: f64) -> f64 {
+float64_beta :: proc(alpha, beta: f64, gen := context.random_generator) -> f64 {
 	if alpha <= 0 || beta <= 0 {
 		panic(#procedure + ": alpha and beta must be > 0.0")
 	}
 	// Knuth Vol 2 Ed 3 pg 134 "the beta distribution"
-	y := float64_gamma(alpha, 1.0)
+	y := float64_gamma(alpha, 1.0, gen)
 	if y != 0 {
-		return y / (y + float64_gamma(beta, 1.0))
+		return y / (y + float64_gamma(beta, 1.0, gen))
 	}
 	return 0
 }
@@ -190,35 +190,35 @@ float64_beta :: proc(alpha, beta: f64) -> f64 {
 //
 // Return values range between 0 and 1
 @(require_results)
-float32_beta :: proc(alpha, beta: f32) -> f32 {
+float32_beta :: proc(alpha, beta: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_beta(f64(alpha), f64(beta)))
+	return f32(float64_beta(f64(alpha), f64(beta), gen))
 }
 
 
 // Pareto distribution, `alpha` is the shape parameter.
 // https://wikipedia.org/wiki/Pareto_distribution
 @(require_results)
-float64_pareto :: proc(alpha: f64) -> f64 {
+float64_pareto :: proc(alpha: f64, gen := context.random_generator) -> f64 {
-	return math.pow(1 - float64(), -1.0 / alpha)
+	return math.pow(1 - float64(gen), -1.0 / alpha)
 }
 // Pareto distribution, `alpha` is the shape parameter.
 // https://wikipedia.org/wiki/Pareto_distribution
 @(require_results)
-float32_pareto :: proc(alpha, beta: f32) -> f32 {
+float32_pareto :: proc(alpha, beta: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_pareto(f64(alpha)))
+	return f32(float64_pareto(f64(alpha), gen))
 }
 
 
 // Weibull distribution, `alpha` is the scale parameter, `beta` is the shape parameter.
 @(require_results)
-float64_weibull :: proc(alpha, beta: f64) -> f64 {
+float64_weibull :: proc(alpha, beta: f64, gen := context.random_generator) -> f64 {
-	u := 1 - float64()
+	u := 1 - float64(gen)
 	return alpha * math.pow(-math.ln(u), 1.0/beta)
 }
 // Weibull distribution, `alpha` is the scale parameter, `beta` is the shape parameter.
 @(require_results)
-float32_weibull :: proc(alpha, beta: f32) -> f32 {
+float32_weibull :: proc(alpha, beta: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_weibull(f64(alpha), f64(beta)))
+	return f32(float64_weibull(f64(alpha), f64(beta), gen))
 }
 
 
@@ -227,23 +227,23 @@ float32_weibull :: proc(alpha, beta: f32) -> f32 {
 // `kappa` is the concentration parameter which must be >= 0
 // When `kappa` is zero, the Distribution is a uniform Distribution over the range 0 to 2pi
 @(require_results)
-float64_von_mises :: proc(mean_angle, kappa: f64) -> f64 {
+float64_von_mises :: proc(mean_angle, kappa: f64, gen := context.random_generator) -> f64 {
 	// Fisher, N.I., "Statistical Analysis of Circular Data", Cambridge University Press, 1993.
 
 	mu := mean_angle
 	if kappa <= 1e-6 {
-		return math.TAU * float64()
+		return math.TAU * float64(gen)
 	}
 
 	s := 0.5 / kappa
 	t := s + math.sqrt(1 + s*s)
 	z: f64
 	for {
-		u1 := float64()
+		u1 := float64(gen)
 		z = math.cos(math.TAU * 0.5 * u1)
 
 		d := z / (t + z)
-		u2 := float64()
+		u2 := float64(gen)
 		if u2 < 1 - d*d || u2 <= (1-d)*math.exp(d) {
 			break
 		}
@@ -251,7 +251,7 @@ float64_von_mises :: proc(mean_angle, kappa: f64) -> f64 {
 
 	q := 1.0 / t
 	f := (q + z) / (1 + q*z)
-	u3 := float64()
+	u3 := float64(gen)
 	if u3 > 0.5 {
 		return math.mod(mu + math.acos(f), math.TAU)
 	} else {
@@ -263,57 +263,57 @@ float64_von_mises :: proc(mean_angle, kappa: f64) -> f64 {
 // `kappa` is the concentration parameter which must be >= 0
 // When `kappa` is zero, the Distribution is a uniform Distribution over the range 0 to 2pi
 @(require_results)
-float32_von_mises :: proc(mean_angle, kappa: f32) -> f32 {
+float32_von_mises :: proc(mean_angle, kappa: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_von_mises(f64(mean_angle), f64(kappa)))
+	return f32(float64_von_mises(f64(mean_angle), f64(kappa), gen))
 }
 
 
 // Cauchy-Lorentz Distribution
 // `x_0` is the location, `gamma` is the scale where `gamma` > 0
 @(require_results)
-float64_cauchy_lorentz :: proc(x_0, gamma: f64) -> f64 {
+float64_cauchy_lorentz :: proc(x_0, gamma: f64, gen := context.random_generator) -> f64 {
 	assert(gamma > 0)
 
 	// Calculated from the inverse CDF
 
-	return math.tan(math.PI * (float64() - 0.5))*gamma + x_0
+	return math.tan(math.PI * (float64(gen) - 0.5))*gamma + x_0
 }
 // Cauchy-Lorentz Distribution
 // `x_0` is the location, `gamma` is the scale where `gamma` > 0
 @(require_results)
-float32_cauchy_lorentz :: proc(x_0, gamma: f32) -> f32 {
+float32_cauchy_lorentz :: proc(x_0, gamma: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_cauchy_lorentz(f64(x_0), f64(gamma)))
+	return f32(float64_cauchy_lorentz(f64(x_0), f64(gamma), gen))
 }
 
 
 // Log Cauchy-Lorentz Distribution
 // `x_0` is the location, `gamma` is the scale where `gamma` > 0
 @(require_results)
-float64_log_cauchy_lorentz :: proc(x_0, gamma: f64) -> f64 {
+float64_log_cauchy_lorentz :: proc(x_0, gamma: f64, gen := context.random_generator) -> f64 {
 	assert(gamma > 0)
-	return math.exp(math.tan(math.PI * (float64() - 0.5))*gamma + x_0)
+	return math.exp(math.tan(math.PI * (float64(gen) - 0.5))*gamma + x_0)
 }
 // Log Cauchy-Lorentz Distribution
 // `x_0` is the location, `gamma` is the scale where `gamma` > 0
 @(require_results)
-float32_log_cauchy_lorentz :: proc(x_0, gamma: f32) -> f32 {
+float32_log_cauchy_lorentz :: proc(x_0, gamma: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_log_cauchy_lorentz(f64(x_0), f64(gamma)))
+	return f32(float64_log_cauchy_lorentz(f64(x_0), f64(gamma), gen))
 }
 
 
 // Laplace Distribution
 // `b` is the scale where `b` > 0
 @(require_results)
-float64_laplace :: proc(mean, b: f64) -> f64 {
+float64_laplace :: proc(mean, b: f64, gen := context.random_generator) -> f64 {
 	assert(b > 0)
-	p := float64()-0.5
+	p := float64(gen)-0.5
 	return -math.sign(p)*math.ln(1 - 2*abs(p))*b + mean
 }
 // Laplace Distribution
 // `b` is the scale where `b` > 0
 @(require_results)
-float32_laplace :: proc(mean, b: f32) -> f32 {
+float32_laplace :: proc(mean, b: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_laplace(f64(mean), f64(b)))
+	return f32(float64_laplace(f64(mean), f64(b), gen))
 }
 
 
@@ -321,18 +321,18 @@ float32_laplace :: proc(mean, b: f32) -> f32 {
 // `eta` is the shape, `b` is the scale
 // Both `eta` and `b` must be > 0
 @(require_results)
-float64_gompertz :: proc(eta, b: f64) -> f64 {
+float64_gompertz :: proc(eta, b: f64, gen := context.random_generator) -> f64 {
 	if eta <= 0 || b <= 0 {
 		panic(#procedure + ": eta and b must be > 0.0")
 	}
 
-	p := float64()
+	p := float64(gen)
 	return math.ln(1 - math.ln(1 - p)/eta)/b
 }
 // Gompertz Distribution
 // `eta` is the shape, `b` is the scale
 // Both `eta` and `b` must be > 0
 @(require_results)
-float32_gompertz :: proc(eta, b: f32) -> f32 {
+float32_gompertz :: proc(eta, b: f32, gen := context.random_generator) -> f32 {
-	return f32(float64_gompertz(f64(eta), f64(b)))
+	return f32(float64_gompertz(f64(eta), f64(b), gen))
 }

core/math/rand/exp.odin

@@ -16,7 +16,7 @@ import "core:math"
 //    https://www.jstatsoft.org/article/view/v005i08 [web page]
 //
 @(require_results)
-exp_float64 :: proc() -> f64 {
+exp_float64 :: proc(gen := context.random_generator) -> f64 {
 	re :: 7.69711747013104972
 
 	@(static, rodata)
@@ -199,16 +199,16 @@ exp_float64 :: proc() -> f64 {
 	}
 
 	for {
-		j := uint32()
+		j := uint32(gen)
 		i := j & 0xFF
 		x := f64(j) * f64(we[i])
 		if j < ke[i] {
 			return x
 		}
 		if i == 0 {
-			return re - math.ln(float64())
+			return re - math.ln(float64(gen))
 		}
-		if fe[i]+f32(float64())*(fe[i-1]-fe[i]) < f32(math.exp(-x)) {
+		if fe[i]+f32(float64(gen))*(fe[i-1]-fe[i]) < f32(math.exp(-x)) {
 			return x
 		}
 	}

core/math/rand/normal.odin

@@ -18,7 +18,7 @@ import "core:math"
 //    https://www.jstatsoft.org/article/view/v005i08 [web page]
 //
 @(require_results)
-norm_float64 :: proc() -> f64 {
+norm_float64 :: proc(gen := context.random_generator) -> f64 {
 	rn :: 3.442619855899
 
 	@(static, rodata)
@@ -116,7 +116,7 @@ norm_float64 :: proc() -> f64 {
 	}
 
 	for {
-		j := i32(uint32())
+		j := i32(uint32(gen))
 		i := j & 0x7f
 		x := f64(j) * f64(wn[i])
 		if u32(abs(j)) < kn[i] {
@@ -126,15 +126,15 @@ norm_float64 :: proc() -> f64 {
 
 		if i == 0 {
 			for {
-				x = -math.ln(float64()) * (1.0/ rn)
+				x = -math.ln(float64(gen)) * (1.0/ rn)
-				y := -math.ln(float64())
+				y := -math.ln(float64(gen))
 				if y+y >= x*x {
 					break
 				}
 			}
 			return j > 0 ? rn + x : -rn - x
 		}
-		if fn[i]+f32(float64())*(fn[i-1]-fn[i]) < f32(math.exp(-0.5*x*x)) {
+		if fn[i]+f32(float64(gen))*(fn[i-1]-fn[i]) < f32(math.exp(-0.5*x*x)) {
 			return x
 		}
 	}

core/math/rand/rand.odin

@@ -62,15 +62,15 @@ Possible Output:
 
 	10
 */
-reset :: proc(seed: u64) {
+reset :: proc(seed: u64, gen := context.random_generator) {
-	runtime.random_generator_reset_u64(context.random_generator, seed)
+	runtime.random_generator_reset_u64(gen, seed)
 }
 
 
 @(private)
-_random_u64 :: proc() -> (res: u64) {
+_random_u64 :: proc(gen := context.random_generator) -> (res: u64) {
-	ok := runtime.random_generator_read_ptr(context.random_generator, &res, size_of(res))
+	ok := runtime.random_generator_read_ptr(gen, &res, size_of(res))
-	assert(ok, "uninitialized context.random_generator")
+	assert(ok, "uninitialized gen/context.random_generator")
 	return
 }
 
@@ -95,7 +95,7 @@ Possible Output:
 
 */
 @(require_results)
-uint32 :: proc() -> (val: u32) { return u32(_random_u64()) }
+uint32 :: proc(gen := context.random_generator) -> (val: u32) { return u32(_random_u64(gen)) }
 
 /*
 Generates a random 64 bit value using the provided random number generator. If no generator is provided the global random number generator will be used.
@@ -118,7 +118,7 @@ Possible Output:
 
 */
 @(require_results)
-uint64 :: proc() -> (val: u64) { return _random_u64() }
+uint64 :: proc(gen := context.random_generator) -> (val: u64) { return _random_u64(gen) }
 
 /*
 Generates a random 128 bit value using the provided random number generator. If no generator is provided the global random number generator will be used.
@@ -141,9 +141,9 @@ Possible Output:
 
 */
 @(require_results)
-uint128 :: proc() -> (val: u128) {
+uint128 :: proc(gen := context.random_generator) -> (val: u128) {
-	a := u128(_random_u64())
+	a := u128(_random_u64(gen))
-	b := u128(_random_u64())
+	b := u128(_random_u64(gen))
 	return (a<<64) | b
 }
 
@@ -168,7 +168,7 @@ Possible Output:
 	389
 
 */
-@(require_results) int31  :: proc() -> (val: i32)  { return i32(uint32() << 1 >> 1) }
+@(require_results) int31  :: proc(gen := context.random_generator) -> (val: i32)  { return i32(uint32(gen) << 1 >> 1) }
 
 /*
 Generates a random 63 bit value using the provided random number generator. If no generator is provided the global random number generator will be used.  
@@ -191,7 +191,7 @@ Possible Output:
 	389
 
 */
-@(require_results) int63  :: proc() -> (val: i64)  { return i64(uint64() << 1 >> 1) }
+@(require_results) int63  :: proc(gen := context.random_generator) -> (val: i64)  { return i64(uint64(gen) << 1 >> 1) }
 
 /*
 Generates a random 127 bit value using the provided random number generator. If no generator is provided the global random number generator will be used.  
@@ -214,7 +214,7 @@ Possible Output:
 	389
 
 */
-@(require_results) int127 :: proc() -> (val: i128) { return i128(uint128() << 1 >> 1) }
+@(require_results) int127 :: proc(gen := context.random_generator) -> (val: i128) { return i128(uint128(gen) << 1 >> 1) }
 
 /*
 Generates a random 31 bit value in the range `[0, n)` using the provided random number generator. If no generator is provided the global random number generator will be used.
@@ -242,17 +242,17 @@ Possible Output:
 
 */
 @(require_results)
-int31_max :: proc(n: i32) -> (val: i32) {
+int31_max :: proc(n: i32, gen := context.random_generator) -> (val: i32) {
 	if n <= 0 {
 		panic("Invalid argument to int31_max")
 	}
 	if n&(n-1) == 0 {
-		return int31() & (n-1)
+		return int31(gen) & (n-1)
 	}
 	max := i32((1<<31) - 1 - (1<<31)%u32(n))
-	v := int31()
+	v := int31(gen)
 	for v > max {
-		v = int31()
+		v = int31(gen)
 	}
 	return v % n
 }
@@ -283,17 +283,17 @@ Possible Output:
 
 */
 @(require_results)
-int63_max :: proc(n: i64) -> (val: i64) {
+int63_max :: proc(n: i64, gen := context.random_generator) -> (val: i64) {
 	if n <= 0 {
 		panic("Invalid argument to int63_max")
 	}
 	if n&(n-1) == 0 {
-		return int63() & (n-1)
+		return int63(gen) & (n-1)
 	}
 	max := i64((1<<63) - 1 - (1<<63)%u64(n))
-	v := int63()
+	v := int63(gen)
 	for v > max {
-		v = int63()
+		v = int63(gen)
 	}
 	return v % n
 }
@@ -324,17 +324,17 @@ Possible Output:
 
 */
 @(require_results)
-int127_max :: proc(n: i128) -> (val: i128) {
+int127_max :: proc(n: i128, gen := context.random_generator) -> (val: i128) {
 	if n <= 0 {
 		panic("Invalid argument to int127_max")
 	}
 	if n&(n-1) == 0 {
-		return int127() & (n-1)
+		return int127(gen) & (n-1)
 	}
 	max := i128((1<<127) - 1 - (1<<127)%u128(n))
-	v := int127()
+	v := int127(gen)
 	for v > max {
-		v = int127()
+		v = int127(gen)
 	}
 	return v % n
 }
@@ -365,14 +365,14 @@ Possible Output:
 
 */
 @(require_results)
-int_max :: proc(n: int) -> (val: int) {
+int_max :: proc(n: int, gen := context.random_generator) -> (val: int) {
 	if n <= 0 {
 		panic("Invalid argument to int_max")
 	}
 	when size_of(int) == 4 {
-		return int(int31_max(i32(n)))
+		return int(int31_max(i32(n), gen))
 	} else {
-		return int(int63_max(i64(n)))
+		return int(int63_max(i64(n), gen))
 	}
 }
 
@@ -396,7 +396,7 @@ Possible Output:
 	0.511
 
 */
-@(require_results) float64 :: proc() -> (val: f64) { return f64(int63_max(1<<53)) / (1 << 53) }
+@(require_results) float64 :: proc(gen := context.random_generator) -> (val: f64) { return f64(int63_max(1<<53, gen)) / (1 << 53) }
 
 /*
 Generates a random single floating point value in the range `[0, 1)` using the provided random number generator. If no generator is provided the global random number generator will be used.
@@ -418,7 +418,7 @@ Possible Output:
 	0.511
 
 */
-@(require_results) float32 :: proc() -> (val: f32) { return f32(int31_max(1<<24)) / (1 << 24) }
+@(require_results) float32 :: proc(gen := context.random_generator) -> (val: f32) { return f32(int31_max(1<<24, gen)) / (1 << 24) }
 
 /*
 Generates a random double floating point value in the range `[low, high)` using the provided random number generator. If no generator is provided the global random number generator will be used.
@@ -446,9 +446,9 @@ Possible Output:
 	673.130
 
 */
-@(require_results) float64_range :: proc(low, high: f64) -> (val: f64) {
+@(require_results) float64_range :: proc(low, high: f64, gen := context.random_generator) -> (val: f64) {
 	assert(low <= high, "low must be lower than or equal to high")
-	val = (high-low)*float64() + low
+	val = (high-low)*float64(gen) + low
 	if val >= high {
 		val = max(low, high * (1 - math.F64_EPSILON))
 	}
@@ -481,9 +481,9 @@ Possible Output:
 	673.130
 
 */
-@(require_results) float32_range :: proc(low, high: f32) -> (val: f32) {
+@(require_results) float32_range :: proc(low, high: f32, gen := context.random_generator) -> (val: f32) {
 	assert(low <= high, "low must be lower than or equal to high")
-	val = (high-low)*float32() + low
+	val = (high-low)*float32(gen) + low
 	if val >= high {
 		val = max(low, high * (1 - math.F32_EPSILON))
 	}
@@ -518,12 +518,12 @@ Possible Output:
 
 */
 @(require_results)
-read :: proc(p: []byte) -> (n: int) {
+read :: proc(p: []byte, gen := context.random_generator) -> (n: int) {
 	pos := i8(0)
 	val := i64(0)
 	for n = 0; n < len(p); n += 1 {
 		if pos == 0 {
-			val = int63()
+			val = int63(gen)
 			pos = 7
 		}
 		p[n] = byte(val)
@@ -566,10 +566,10 @@ Possible Output:
 
 */
 @(require_results)
-perm :: proc(n: int, allocator := context.allocator) -> (res: []int, err: mem.Allocator_Error) #optional_allocator_error {
+perm :: proc(n: int, allocator := context.allocator, gen := context.random_generator) -> (res: []int, err: mem.Allocator_Error) #optional_allocator_error {
 	m := make([]int, n, allocator) or_return
 	for i := 0; i < n; i += 1 {
-		j := int_max(i+1)
+		j := int_max(i+1, gen)
 		m[i] = m[j]
 		m[j] = i
 	}
@@ -599,14 +599,14 @@ Possible Output:
 	[2, 4, 3, 1]
 
 */
-shuffle :: proc(array: $T/[]$E) {
+shuffle :: proc(array: $T/[]$E, gen := context.random_generator) {
 	n := i64(len(array))
 	if n < 2 {
 		return
 	}
 
 	for i := i64(n - 1); i > 0; i -= 1 {
-		j := int63_max(i + 1)
+		j := int63_max(i + 1, gen)
 		array[i], array[j] = array[j], array[i]
 	}
 }
@@ -641,17 +641,17 @@ Possible Output:
 
 */
 @(require_results)
-choice :: proc(array: $T/[]$E) -> (res: E) {
+choice :: proc(array: $T/[]$E, gen := context.random_generator) -> (res: E) {
 	n := i64(len(array))
 	if n < 1 {
 		return E{}
 	}
-	return array[int63_max(n)]
+	return array[int63_max(n, gen)]
 }
 
 
 @(require_results)
-choice_enum :: proc($T: typeid) -> T
+choice_enum :: proc($T: typeid, gen := context.random_generator) -> T
 	where
 		intrinsics.type_is_enum(T),
 		size_of(T) <= 8,
@@ -659,11 +659,11 @@ choice_enum :: proc($T: typeid) -> T
 {
 	when intrinsics.type_is_unsigned(intrinsics.type_core_type(T)) &&
 	     u64(max(T)) > u64(max(i64)) {
-		i := uint64() % u64(len(T))
+		i := uint64(gen) % u64(len(T))
 		i += u64(min(T))
 		return T(i)
 	} else {
-		i := int63_max(i64(len(T)))
+		i := int63_max(i64(len(T)), gen)
 		i += i64(min(T))
 		return T(i)
 	}