Improve package math/rand

core/math/math.odin

@@ -61,9 +61,15 @@ foreign _ {
 	log_f32 :: proc(x: f32) -> f32 ---;
 	@(link_name="llvm.log.f64")
 	log_f64 :: proc(x: f64) -> f64 ---;
+
+	@(link_name="llvm.exp.f32")
+	exp_f32 :: proc(x: f32) -> f32 ---;
+	@(link_name="llvm.exp.f64")
+	exp_f64 :: proc(x: f64) -> f64 ---;
 }
 
 log :: proc{log_f32, log_f64};
+exp :: proc{exp_f32, exp_f64};
 
 tan_f32 :: proc "c" (θ: f32) -> f32 { return sin(θ)/cos(θ); }
 tan_f64 :: proc "c" (θ: f64) -> f64 { return sin(θ)/cos(θ); }

core/math/rand/normal.odin

@@ -0,0 +1,146 @@
+package rand
+
+import "core:math"
+
+//
+//    Normal distribution
+//
+//    "The Ziggurat Method for Generating Random Variables"
+//    Authors: George Marsaglia, Wai Wan Tsang
+//    Submitted: 2000-04-15. Published: 2000-10-02.
+//    https://www.jstatsoft.org/index.php/jss/article/view/v005i08/ziggurat.pdf [pdf]
+//    https://www.jstatsoft.org/article/view/v005i08 [web page]
+//
+
+// norm_float64 returns a normally distributed f64 in the range -max(f64) through +max(f64) inclusive,
+// with a standard normal distribution with a mean of 0 and standard deviation of 1.
+//
+// sample = norm_float64() * std_dev + mean
+//
+norm_float64 :: proc(r: ^Rand = global_rand_ptr) -> f64 {
+	rn :: 3.442619855899;
+
+	@(static)
+	kn := [128]u32{
+		0x76ad2212, 0x00000000, 0x600f1b53, 0x6ce447a6, 0x725b46a2,
+		0x7560051d, 0x774921eb, 0x789a25bd, 0x799045c3, 0x7a4bce5d,
+		0x7adf629f, 0x7b5682a6, 0x7bb8a8c6, 0x7c0ae722, 0x7c50cce7,
+		0x7c8cec5b, 0x7cc12cd6, 0x7ceefed2, 0x7d177e0b, 0x7d3b8883,
+		0x7d5bce6c, 0x7d78dd64, 0x7d932886, 0x7dab0e57, 0x7dc0dd30,
+		0x7dd4d688, 0x7de73185, 0x7df81cea, 0x7e07c0a3, 0x7e163efa,
+		0x7e23b587, 0x7e303dfd, 0x7e3beec2, 0x7e46db77, 0x7e51155d,
+		0x7e5aabb3, 0x7e63abf7, 0x7e6c222c, 0x7e741906, 0x7e7b9a18,
+		0x7e82adfa, 0x7e895c63, 0x7e8fac4b, 0x7e95a3fb, 0x7e9b4924,
+		0x7ea0a0ef, 0x7ea5b00d, 0x7eaa7ac3, 0x7eaf04f3, 0x7eb3522a,
+		0x7eb765a5, 0x7ebb4259, 0x7ebeeafd, 0x7ec2620a, 0x7ec5a9c4,
+		0x7ec8c441, 0x7ecbb365, 0x7ece78ed, 0x7ed11671, 0x7ed38d62,
+		0x7ed5df12, 0x7ed80cb4, 0x7eda175c, 0x7edc0005, 0x7eddc78e,
+		0x7edf6ebf, 0x7ee0f647, 0x7ee25ebe, 0x7ee3a8a9, 0x7ee4d473,
+		0x7ee5e276, 0x7ee6d2f5, 0x7ee7a620, 0x7ee85c10, 0x7ee8f4cd,
+		0x7ee97047, 0x7ee9ce59, 0x7eea0eca, 0x7eea3147, 0x7eea3568,
+		0x7eea1aab, 0x7ee9e071, 0x7ee98602, 0x7ee90a88, 0x7ee86d08,
+		0x7ee7ac6a, 0x7ee6c769, 0x7ee5bc9c, 0x7ee48a67, 0x7ee32efc,
+		0x7ee1a857, 0x7edff42f, 0x7ede0ffa, 0x7edbf8d9, 0x7ed9ab94,
+		0x7ed7248d, 0x7ed45fae, 0x7ed1585c, 0x7ece095f, 0x7eca6ccb,
+		0x7ec67be2, 0x7ec22eee, 0x7ebd7d1a, 0x7eb85c35, 0x7eb2c075,
+		0x7eac9c20, 0x7ea5df27, 0x7e9e769f, 0x7e964c16, 0x7e8d44ba,
+		0x7e834033, 0x7e781728, 0x7e6b9933, 0x7e5d8a1a, 0x7e4d9ded,
+		0x7e3b737a, 0x7e268c2f, 0x7e0e3ff5, 0x7df1aa5d, 0x7dcf8c72,
+		0x7da61a1e, 0x7d72a0fb, 0x7d30e097, 0x7cd9b4ab, 0x7c600f1a,
+		0x7ba90bdc, 0x7a722176, 0x77d664e5,
+	};
+
+	@(static)
+	wn := [128]f32{
+		1.7290405e-09, 1.2680929e-10, 1.6897518e-10, 1.9862688e-10,
+		2.2232431e-10, 2.4244937e-10, 2.601613e-10,  2.7611988e-10,
+		2.9073963e-10, 3.042997e-10,  3.1699796e-10, 3.289802e-10,
+		3.4035738e-10, 3.5121603e-10, 3.616251e-10,  3.7164058e-10,
+		3.8130857e-10, 3.9066758e-10, 3.9975012e-10, 4.08584e-10,
+		4.1719309e-10, 4.2559822e-10, 4.338176e-10,  4.418672e-10,
+		4.497613e-10,  4.5751258e-10, 4.651324e-10,  4.7263105e-10,
+		4.8001775e-10, 4.87301e-10,   4.944885e-10,  5.015873e-10,
+		5.0860405e-10, 5.155446e-10,  5.2241467e-10, 5.2921934e-10,
+		5.359635e-10,  5.426517e-10,  5.4928817e-10, 5.5587696e-10,
+		5.624219e-10,  5.6892646e-10, 5.753941e-10,  5.818282e-10,
+		5.882317e-10,  5.946077e-10,  6.00959e-10,   6.072884e-10,
+		6.135985e-10,  6.19892e-10,   6.2617134e-10, 6.3243905e-10,
+		6.386974e-10,  6.449488e-10,  6.511956e-10,  6.5744005e-10,
+		6.6368433e-10, 6.699307e-10,  6.7618144e-10, 6.824387e-10,
+		6.8870465e-10, 6.949815e-10,  7.012715e-10,  7.075768e-10,
+		7.1389966e-10, 7.202424e-10,  7.266073e-10,  7.329966e-10,
+		7.394128e-10,  7.4585826e-10, 7.5233547e-10, 7.58847e-10,
+		7.653954e-10,  7.719835e-10,  7.7861395e-10, 7.852897e-10,
+		7.920138e-10,  7.987892e-10,  8.0561924e-10, 8.125073e-10,
+		8.194569e-10,  8.2647167e-10, 8.3355556e-10, 8.407127e-10,
+		8.479473e-10,  8.55264e-10,   8.6266755e-10, 8.7016316e-10,
+		8.777562e-10,  8.8545243e-10, 8.932582e-10,  9.0117996e-10,
+		9.09225e-10,   9.174008e-10,  9.2571584e-10, 9.341788e-10,
+		9.427997e-10,  9.515889e-10,  9.605579e-10,  9.697193e-10,
+		9.790869e-10,  9.88676e-10,   9.985036e-10,  1.0085882e-09,
+		1.0189509e-09, 1.0296151e-09, 1.0406069e-09, 1.0519566e-09,
+		1.063698e-09,  1.0758702e-09, 1.0885183e-09, 1.1016947e-09,
+		1.1154611e-09, 1.1298902e-09, 1.1450696e-09, 1.1611052e-09,
+		1.1781276e-09, 1.1962995e-09, 1.2158287e-09, 1.2369856e-09,
+		1.2601323e-09, 1.2857697e-09, 1.3146202e-09, 1.347784e-09,
+		1.3870636e-09, 1.4357403e-09, 1.5008659e-09, 1.6030948e-09,
+	};
+
+	@(static)
+	fn := [128]f32{
+		1.00000000,  0.9635997,   0.9362827,   0.9130436,   0.89228165,
+		0.87324303,  0.8555006,   0.8387836,   0.8229072,   0.8077383,
+		0.793177,    0.7791461,   0.7655842,   0.7524416,   0.73967725,
+		0.7272569,   0.7151515,   0.7033361,   0.69178915,  0.68049186,
+		0.6694277,   0.658582,    0.6479418,   0.63749546,  0.6272325,
+		0.6171434,   0.6072195,   0.5974532,   0.58783704,  0.5783647,
+		0.56903,     0.5598274,   0.5507518,   0.54179835,  0.5329627,
+		0.52424055,  0.5156282,   0.50712204,  0.49871865,  0.49041483,
+		0.48220766,  0.4740943,   0.46607214,  0.4581387,   0.45029163,
+		0.44252872,  0.43484783,  0.427247,    0.41972435,  0.41227803,
+		0.40490642,  0.39760786,  0.3903808,   0.3832238,   0.37613547,
+		0.36911446,  0.3621595,   0.35526937,  0.34844297,  0.34167916,
+		0.33497685,  0.3283351,   0.3217529,   0.3152294,   0.30876362,
+		0.30235484,  0.29600215,  0.28970486,  0.2834622,   0.2772735,
+		0.27113807,  0.2650553,   0.25902456,  0.2530453,   0.24711695,
+		0.241239,    0.23541094,  0.22963232,  0.2239027,   0.21822165,
+		0.21258877,  0.20700371,  0.20146611,  0.19597565,  0.19053204,
+		0.18513499,  0.17978427,  0.17447963,  0.1692209,   0.16400786,
+		0.15884037,  0.15371831,  0.14864157,  0.14361008,  0.13862377,
+		0.13368265,  0.12878671,  0.12393598,  0.119130544, 0.11437051,
+		0.10965602,  0.104987256, 0.10036444,  0.095787846, 0.0912578,
+		0.08677467,  0.0823389,   0.077950984, 0.073611505, 0.06932112,
+		0.06508058,  0.06089077,  0.056752663, 0.0526674,   0.048636295,
+		0.044660863, 0.040742867, 0.03688439,  0.033087887, 0.029356318,
+		0.025693292, 0.022103304, 0.018592102, 0.015167298, 0.011839478,
+		0.008624485, 0.005548995, 0.0026696292,
+	};
+
+
+	for {
+		j := i32(uint32(r));
+		i := j & 0x7f;
+		x := f64(j) * f64(wn[i]);
+		if u32(abs(j)) < kn[i] {
+			// 99% of the time this will be hit
+			return x;
+		}
+
+		if i == 0 {
+			for {
+				x = -math.log(float64(r)) * (1.0/ rn);
+				y := -math.log(float64(r));
+				if y+y >= x*x {
+					break;
+				}
+			}
+			return j > 0 ? rn + x : -rn - x;
+		}
+		if fn[i]+f32(float64(r))*(fn[i-1]-fn[i]) < f32(math.exp(-0.5*x*x)) {
+			return x;
+		}
+	}
+
+	return 0; // NOTE(bill): Will never be hit but this is here for sanity's sake
+}
+

core/math/rand/rand.odin

@@ -5,7 +5,25 @@ Rand :: struct {
 	inc:   u64,
 }
 
-init :: proc(r: ^Rand, seed: u64 = 8675309) {
+
+@(private, static)
+_GLOBAL_SEED_DATA := 1234567890;
+@(private, static)
+global_rand := create(u64(uintptr(&_GLOBAL_SEED_DATA)));
+@(private, static)
+global_rand_ptr := &global_rand;
+
+set_global_seed :: proc(seed: u64) {
+	init(global_rand_ptr, seed);
+}
+
+create :: proc(seed: u64) -> Rand {
+	r: Rand;
+	init(&r, seed);
+	return r;
+}
+
+init :: proc(r: ^Rand, seed: u64) {
 	r.state = 0;
 	r.inc = (seed << 1) | 1;
 	_random(r);
@@ -21,18 +39,18 @@ _random :: proc(r: ^Rand) -> u32 {
 	return (xor_shifted >> rot) | (xor_shifted << ((-rot) & 31));
 }
 
-uint32 :: proc(r: ^Rand) -> u32 { return _random(r); }
+uint32 :: proc(r: ^Rand = global_rand_ptr) -> u32 { return _random(r); }
 
-uint64 :: proc(r: ^Rand) -> u64 {
+uint64 :: proc(r: ^Rand = global_rand_ptr) -> u64 {
 	a := u64(_random(r));
 	b := u64(_random(r));
 	return (a<<32) | b;
 }
 
-int31 :: proc(r: ^Rand) -> i32 { return i32(uint32(r) << 1 >> 1); }
-int63 :: proc(r: ^Rand) -> i64 { return i64(uint64(r) << 1 >> 1); }
+int31 :: proc(r: ^Rand = global_rand_ptr) -> i32 { return i32(uint32(r) << 1 >> 1); }
+int63 :: proc(r: ^Rand = global_rand_ptr) -> i64 { return i64(uint64(r) << 1 >> 1); }
 
-int31_max :: proc(r: ^Rand, n: i32) -> i32 {
+int31_max :: proc(n: i32, r: ^Rand = global_rand_ptr) -> i32 {
 	if n <= 0 do panic("Invalid argument to int31_max");
 	if n&(n-1) == 0 {
 		return int31(r) & (n-1);
@@ -45,7 +63,7 @@ int31_max :: proc(r: ^Rand, n: i32) -> i32 {
 	return v % n;
 }
 
-int63_max :: proc(r: ^Rand, n: i64) -> i64 {
+int63_max :: proc(n: i64, r: ^Rand = global_rand_ptr) -> i64 {
 	if n <= 0 do panic("Invalid argument to int63_max");
 	if n&(n-1) == 0 {
 		return int63(r) & (n-1);
@@ -58,5 +76,52 @@ int63_max :: proc(r: ^Rand, n: i64) -> i64 {
 	return v % n;
 }
 
-float64 :: proc(r: ^Rand) -> f64 { return f64(int63_max(r, 1<<53)) / (1 << 53); }
-float32 :: proc(r: ^Rand) -> f32 { return f32(float64(r)); }
+int_max :: proc(n: int, r: ^Rand = global_rand_ptr) -> int {
+	if n <= 0 do panic("Invalid argument to int_max");
+	when size_of(int) == 4 {
+		return int(int31_max(i32(n), r));
+	} else {
+		return int(int63_max(i64(n), r));
+	}
+}
+
+float64 :: proc(r: ^Rand = global_rand_ptr) -> f64 { return f64(int63_max(1<<53, r)) / (1 << 53); }
+float32 :: proc(r: ^Rand = global_rand_ptr) -> f32 { return f32(float64(r)); }
+
+
+read :: proc(p: []byte, r: ^Rand = global_rand_ptr) -> (n: int) {
+	pos := i8(0);
+	val := i64(0);
+	for n = 0; n < len(p); n += 1 {
+		if pos == 0 {
+			val = int63(r);
+			pos = 7;
+		}
+		p[n] = byte(val);
+		val >>= 8;
+		pos -= 1;
+	}
+	return;
+}
+
+// perm returns a slice of n ints in a pseudo-random permutation of integers in the range [0, n)
+perm :: proc(n: int, r: ^Rand = global_rand_ptr) -> []int {
+	m := make([]int, n);
+	for i := 0; i < n; i += 1 {
+		j := int_max(i+1);
+		m[i] = m[j];
+		m[j] = i;
+	}
+	return m;
+}
+
+
+shuffle :: proc(array: $T/[]$E, r: ^Rand = global_rand_ptr) {
+	n := i64(len(array));
+	if n < 2 do return;
+
+	for i := i64(0); i < n; i += 1 {
+		j := int63_max(n, r);
+		array[i], array[j] = array[j], array[i];
+	}
+}