Add RNGs to the language

17 files changed, 554 insertions, 125 deletions

diff --git a/stdlib/arrays.c b/stdlib/arrays.c
index 552fb4cb..8814f31a 100644
--- a/stdlib/arrays.c
+++ b/stdlib/arrays.c
@@ -8,6 +8,7 @@
 #include "arrays.h"
 #include "metamethods.h"
 #include "optionals.h"
+#include "rng.h"
 #include "tables.h"
 #include "text.h"
 #include "util.h"
@@ -249,51 +250,34 @@ public Array_t Array$sorted(Array_t arr, Closure_t comparison, int64_t padded_it
     return arr;
 }
 
-static uint64_t random_range(Closure_t rng, uint64_t upper_bound)
-{
-    if (upper_bound < 2)
-        return 0;
-
-    // This approach is taken from arc4random_uniform()
-    uint64_t min = -upper_bound % upper_bound;
-    uint64_t r;
-    for (;;) {
-        r = ((uint64_t(*)(void*))rng.fn)(rng.userdata);
-        if (r >= min)
-            break;
-    }
-
-    return r % upper_bound;
-}
-
-public void Array$shuffle(Array_t *arr, Closure_t rng, int64_t padded_item_size)
+public void Array$shuffle(Array_t *arr, RNG_t rng, int64_t padded_item_size)
 {
     if (arr->data_refcount != 0 || (int64_t)arr->stride != padded_item_size)
         Array$compact(arr, padded_item_size);
 
     char tmp[padded_item_size];
     for (int64_t i = arr->length-1; i > 1; i--) {
-        int64_t j = (int64_t)random_range(rng, (uint64_t)(i+1));
+        int64_t j = RNG$int64(rng, 0, i);
         memcpy(tmp, arr->data + i*padded_item_size, (size_t)padded_item_size);
         memcpy((void*)arr->data + i*padded_item_size, arr->data + j*padded_item_size, (size_t)padded_item_size);
         memcpy((void*)arr->data + j*padded_item_size, tmp, (size_t)padded_item_size);
     }
 }
 
-public Array_t Array$shuffled(Array_t arr, Closure_t rng, int64_t padded_item_size)
+public Array_t Array$shuffled(Array_t arr, RNG_t rng, int64_t padded_item_size)
 {
     Array$compact(&arr, padded_item_size);
     Array$shuffle(&arr, rng, padded_item_size);
     return arr;
 }
 
-public void *Array$random(Array_t arr, Closure_t rng)
+public void *Array$random(Array_t arr, RNG_t rng)
 {
     if (arr.length == 0)
         return NULL; // fail("Cannot get a random item from an empty array!");
 
-    uint64_t index = random_range(rng, (uint64_t)arr.length);
-    return arr.data + arr.stride*(int64_t)index;
+    int64_t index = RNG$int64(rng, 0, arr.length-1);
+    return arr.data + arr.stride*index;
 }
 
 public Table_t Array$counts(Array_t arr, const TypeInfo_t *type)
@@ -310,7 +294,7 @@ public Table_t Array$counts(Array_t arr, const TypeInfo_t *type)
     return counts;
 }
 
-public Array_t Array$sample(Array_t arr, Int_t int_n, Array_t weights, int64_t padded_item_size)
+public Array_t Array$sample(Array_t arr, Int_t int_n, Array_t weights, RNG_t rng, int64_t padded_item_size)
 {
     int64_t n = Int_to_Int64(int_n, false);
     if (n < 0)
@@ -329,7 +313,7 @@ public Array_t Array$sample(Array_t arr, Int_t int_n, Array_t weights, int64_t p
 
     if (weights.length < 0) {
         for (int64_t i = 0; i < n; i++) {
-            int64_t index = arc4random_uniform(arr.length);
+            int64_t index = RNG$int64(rng, 0, arr.length-1);
             memcpy(selected.data + i*padded_item_size, arr.data + arr.stride*index, (size_t)padded_item_size);
         }
         return selected;
@@ -393,7 +377,7 @@ public Array_t Array$sample(Array_t arr, Int_t int_n, Array_t weights, int64_t p
             aliases[i].alias = i;
 
     for (int64_t i = 0; i < n; i++) {
-        double r = drand48() * arr.length;
+        double r = RNG$num(rng, 0, arr.length);
         int64_t index = (int64_t)r;
         if ((r - (double)index) > aliases[index].odds)
             index = aliases[index].alias;
diff --git a/stdlib/arrays.h b/stdlib/arrays.h
index 5d452d30..251e9f92 100644
--- a/stdlib/arrays.h
+++ b/stdlib/arrays.h
@@ -70,11 +70,11 @@ Int_t Array$find(Array_t arr, void *item, const TypeInfo_t *type);
 Int_t Array$first(Array_t arr, Closure_t predicate);
 void Array$sort(Array_t *arr, Closure_t comparison, int64_t padded_item_size);
 Array_t Array$sorted(Array_t arr, Closure_t comparison, int64_t padded_item_size);
-void Array$shuffle(Array_t *arr, Closure_t rng, int64_t padded_item_size);
-Array_t Array$shuffled(Array_t arr, Closure_t rng, int64_t padded_item_size);
-void *Array$random(Array_t arr, Closure_t rng);
+void Array$shuffle(Array_t *arr, RNG_t rng, int64_t padded_item_size);
+Array_t Array$shuffled(Array_t arr, RNG_t rng, int64_t padded_item_size);
+void *Array$random(Array_t arr, RNG_t rng);
 #define Array$random_value(arr, rng, t) ({ Array_t _arr = arr; if (_arr.length == 0) fail("Cannot get a random value from an empty array!"); *(t*)Array$random(_arr, rng); })
-Array_t Array$sample(Array_t arr, Int_t n, Array_t weights, int64_t padded_item_size);
+Array_t Array$sample(Array_t arr, Int_t n, Array_t weights, RNG_t rng, int64_t padded_item_size);
 Table_t Array$counts(Array_t arr, const TypeInfo_t *type);
 void Array$clear(Array_t *array);
 void Array$compact(Array_t *arr, int64_t padded_item_size);
diff --git a/stdlib/bools.c b/stdlib/bools.c
index c815c733..bc1a6445 100644
--- a/stdlib/bools.c
+++ b/stdlib/bools.c
@@ -39,11 +39,6 @@ PUREFUNC public OptionalBool_t Bool$from_text(Text_t text)
     }
 }
 
-public Bool_t Bool$random(double p)
-{
-    return (drand48() < p); 
-}
-
 public const TypeInfo_t Bool$info = {
     .size=sizeof(bool),
     .align=__alignof__(bool),
diff --git a/stdlib/bools.h b/stdlib/bools.h
index 4b0b42aa..fbc40f8b 100644
--- a/stdlib/bools.h
+++ b/stdlib/bools.h
@@ -15,7 +15,6 @@
 
 PUREFUNC Text_t Bool$as_text(const bool *b, bool colorize, const TypeInfo_t *type);
 OptionalBool_t Bool$from_text(Text_t text);
-Bool_t Bool$random(double p);
 
 extern const TypeInfo_t Bool$info;
 
diff --git a/stdlib/bytes.c b/stdlib/bytes.c
index f7e90576..8d665790 100644
--- a/stdlib/bytes.c
+++ b/stdlib/bytes.c
@@ -17,17 +17,6 @@ PUREFUNC public Text_t Byte$as_text(const Byte_t *b, bool colorize, const TypeIn
     return Text$format(colorize ? "\x1b[35m%u[B]\x1b[m" : "%u[B]", *b);
 }
 
-public Byte_t Byte$random(Byte_t min, Byte_t max)
-{
-    if (min > max)
-        fail("Random minimum value (%u) is larger than the maximum value (%u)", min, max);
-    if (min == max)
-        return min;
-
-    uint32_t r = arc4random_uniform((uint32_t)max - (uint32_t)min + 1u);
-    return (Byte_t)(min + r);
-}
-
 public const TypeInfo_t Byte$info = {
     .size=sizeof(Byte_t),
     .align=__alignof__(Byte_t),
diff --git a/stdlib/bytes.h b/stdlib/bytes.h
index 62b266e7..2a92f658 100644
--- a/stdlib/bytes.h
+++ b/stdlib/bytes.h
@@ -12,7 +12,6 @@
 #define Byte(b) ((Byte_t)(b))
 
 PUREFUNC Text_t Byte$as_text(const Byte_t *b, bool colorize, const TypeInfo_t *type);
-Byte_t Byte$random(Byte_t min, Byte_t max);
 
 extern const Byte_t Byte$min;
 extern const Byte_t Byte$max;
diff --git a/stdlib/chacha.h b/stdlib/chacha.h
new file mode 100644
index 00000000..1bca008d
--- /dev/null
+++ b/stdlib/chacha.h
@@ -0,0 +1,223 @@
+/*
+chacha-merged.c version 20080118
+D. J. Bernstein
+Public domain.
+*/
+
+/* $OpenBSD: chacha_private.h,v 1.3 2022/02/28 21:56:29 dtucker Exp $ */
+/* Tomo: chacha.h,v 1.0 2024/11/03 Bruce Hill */
+
+typedef unsigned char u8;
+typedef unsigned int u32;
+
+typedef struct
+{
+  u32 input[16]; /* could be compressed */
+} chacha_ctx;
+
+#define U8C(v) (v##U)
+#define U32C(v) (v##U)
+
+#define U8V(v) ((u8)(v) & U8C(0xFF))
+#define U32V(v) ((u32)(v) & U32C(0xFFFFFFFF))
+
+#define ROTL32(v, n) \
+  (U32V((v) << (n)) | ((v) >> (32 - (n))))
+
+#define U8TO32_LITTLE(p) \
+  (((u32)((p)[0])      ) | \
+   ((u32)((p)[1]) <<  8) | \
+   ((u32)((p)[2]) << 16) | \
+   ((u32)((p)[3]) << 24))
+
+#define U32TO8_LITTLE(p, v) \
+  do { \
+    (p)[0] = U8V((v)      ); \
+    (p)[1] = U8V((v) >>  8); \
+    (p)[2] = U8V((v) >> 16); \
+    (p)[3] = U8V((v) >> 24); \
+  } while (0)
+
+#define ROTATE(v,c) (ROTL32(v,c))
+#define XOR(v,w) ((v) ^ (w))
+#define PLUS(v,w) (U32V((v) + (w)))
+#define PLUSONE(v) (PLUS((v),1))
+
+#define QUARTERROUND(a,b,c,d) \
+  a = PLUS(a,b); d = ROTATE(XOR(d,a),16); \
+  c = PLUS(c,d); b = ROTATE(XOR(b,c),12); \
+  a = PLUS(a,b); d = ROTATE(XOR(d,a), 8); \
+  c = PLUS(c,d); b = ROTATE(XOR(b,c), 7);
+
+static const char sigma[16] = "expand 32-byte k";
+static const char tau[16] = "expand 16-byte k";
+
+static void
+chacha_keysetup(chacha_ctx *x,const u8 *k,u32 kbits)
+{
+  const char *constants;
+
+  x->input[4] = U8TO32_LITTLE(k + 0);
+  x->input[5] = U8TO32_LITTLE(k + 4);
+  x->input[6] = U8TO32_LITTLE(k + 8);
+  x->input[7] = U8TO32_LITTLE(k + 12);
+  if (kbits == 256) { /* recommended */
+    k += 16;
+    constants = sigma;
+  } else { /* kbits == 128 */
+    constants = tau;
+  }
+  x->input[8] = U8TO32_LITTLE(k + 0);
+  x->input[9] = U8TO32_LITTLE(k + 4);
+  x->input[10] = U8TO32_LITTLE(k + 8);
+  x->input[11] = U8TO32_LITTLE(k + 12);
+  x->input[0] = U8TO32_LITTLE(constants + 0);
+  x->input[1] = U8TO32_LITTLE(constants + 4);
+  x->input[2] = U8TO32_LITTLE(constants + 8);
+  x->input[3] = U8TO32_LITTLE(constants + 12);
+}
+
+static void
+chacha_ivsetup(chacha_ctx *x,const u8 *iv)
+{
+  x->input[12] = 0;
+  x->input[13] = 0;
+  x->input[14] = U8TO32_LITTLE(iv + 0);
+  x->input[15] = U8TO32_LITTLE(iv + 4);
+}
+
+static void
+chacha_encrypt_bytes(chacha_ctx *x,const u8 *m,u8 *c,u32 bytes)
+{
+  u32 x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15;
+  u32 j0, j1, j2, j3, j4, j5, j6, j7, j8, j9, j10, j11, j12, j13, j14, j15;
+  u8 *ctarget = NULL;
+  u8 tmp[64];
+  u_int i;
+
+  if (!bytes) return;
+
+  j0 = x->input[0];
+  j1 = x->input[1];
+  j2 = x->input[2];
+  j3 = x->input[3];
+  j4 = x->input[4];
+  j5 = x->input[5];
+  j6 = x->input[6];
+  j7 = x->input[7];
+  j8 = x->input[8];
+  j9 = x->input[9];
+  j10 = x->input[10];
+  j11 = x->input[11];
+  j12 = x->input[12];
+  j13 = x->input[13];
+  j14 = x->input[14];
+  j15 = x->input[15];
+
+  for (;;) {
+    if (bytes < 64) {
+      for (i = 0;i < bytes;++i) tmp[i] = m[i];
+      m = tmp;
+      ctarget = c;
+      c = tmp;
+    }
+    x0 = j0;
+    x1 = j1;
+    x2 = j2;
+    x3 = j3;
+    x4 = j4;
+    x5 = j5;
+    x6 = j6;
+    x7 = j7;
+    x8 = j8;
+    x9 = j9;
+    x10 = j10;
+    x11 = j11;
+    x12 = j12;
+    x13 = j13;
+    x14 = j14;
+    x15 = j15;
+    for (i = 20;i > 0;i -= 2) {
+      QUARTERROUND( x0, x4, x8,x12)
+      QUARTERROUND( x1, x5, x9,x13)
+      QUARTERROUND( x2, x6,x10,x14)
+      QUARTERROUND( x3, x7,x11,x15)
+      QUARTERROUND( x0, x5,x10,x15)
+      QUARTERROUND( x1, x6,x11,x12)
+      QUARTERROUND( x2, x7, x8,x13)
+      QUARTERROUND( x3, x4, x9,x14)
+    }
+    x0 = PLUS(x0,j0);
+    x1 = PLUS(x1,j1);
+    x2 = PLUS(x2,j2);
+    x3 = PLUS(x3,j3);
+    x4 = PLUS(x4,j4);
+    x5 = PLUS(x5,j5);
+    x6 = PLUS(x6,j6);
+    x7 = PLUS(x7,j7);
+    x8 = PLUS(x8,j8);
+    x9 = PLUS(x9,j9);
+    x10 = PLUS(x10,j10);
+    x11 = PLUS(x11,j11);
+    x12 = PLUS(x12,j12);
+    x13 = PLUS(x13,j13);
+    x14 = PLUS(x14,j14);
+    x15 = PLUS(x15,j15);
+
+#ifndef KEYSTREAM_ONLY
+    x0 = XOR(x0,U8TO32_LITTLE(m + 0));
+    x1 = XOR(x1,U8TO32_LITTLE(m + 4));
+    x2 = XOR(x2,U8TO32_LITTLE(m + 8));
+    x3 = XOR(x3,U8TO32_LITTLE(m + 12));
+    x4 = XOR(x4,U8TO32_LITTLE(m + 16));
+    x5 = XOR(x5,U8TO32_LITTLE(m + 20));
+    x6 = XOR(x6,U8TO32_LITTLE(m + 24));
+    x7 = XOR(x7,U8TO32_LITTLE(m + 28));
+    x8 = XOR(x8,U8TO32_LITTLE(m + 32));
+    x9 = XOR(x9,U8TO32_LITTLE(m + 36));
+    x10 = XOR(x10,U8TO32_LITTLE(m + 40));
+    x11 = XOR(x11,U8TO32_LITTLE(m + 44));
+    x12 = XOR(x12,U8TO32_LITTLE(m + 48));
+    x13 = XOR(x13,U8TO32_LITTLE(m + 52));
+    x14 = XOR(x14,U8TO32_LITTLE(m + 56));
+    x15 = XOR(x15,U8TO32_LITTLE(m + 60));
+#endif
+
+    j12 = PLUSONE(j12);
+    if (!j12) {
+      j13 = PLUSONE(j13);
+      /* stopping at 2^70 bytes per nonce is user's responsibility */
+    }
+
+    U32TO8_LITTLE(c + 0,x0);
+    U32TO8_LITTLE(c + 4,x1);
+    U32TO8_LITTLE(c + 8,x2);
+    U32TO8_LITTLE(c + 12,x3);
+    U32TO8_LITTLE(c + 16,x4);
+    U32TO8_LITTLE(c + 20,x5);
+    U32TO8_LITTLE(c + 24,x6);
+    U32TO8_LITTLE(c + 28,x7);
+    U32TO8_LITTLE(c + 32,x8);
+    U32TO8_LITTLE(c + 36,x9);
+    U32TO8_LITTLE(c + 40,x10);
+    U32TO8_LITTLE(c + 44,x11);
+    U32TO8_LITTLE(c + 48,x12);
+    U32TO8_LITTLE(c + 52,x13);
+    U32TO8_LITTLE(c + 56,x14);
+    U32TO8_LITTLE(c + 60,x15);
+
+    if (bytes <= 64) {
+      if (bytes < 64) {
+        for (i = 0;i < bytes;++i) ctarget[i] = c[i];
+      }
+      x->input[12] = j12;
+      x->input[13] = j13;
+      return;
+    }
+    bytes -= 64;
+    c += 64;
+#ifndef KEYSTREAM_ONLY
+    m += 64;
+#endif
+  }
+}
diff --git a/stdlib/datatypes.h b/stdlib/datatypes.h
index 8d342d2e..f033a21a 100644
--- a/stdlib/datatypes.h
+++ b/stdlib/datatypes.h
@@ -93,4 +93,6 @@ typedef struct Text_s {
 typedef struct timeval DateTime_t;
 #define OptionalDateTime_t DateTime_t
 
+typedef struct RNGState_t* RNG_t;
+
 // vim: ts=4 sw=0 et cino=L2,l1,(0,W4,m1,\:0
diff --git a/stdlib/integers.c b/stdlib/integers.c
index 6b0a9dd2..f6140da4 100644
--- a/stdlib/integers.c
+++ b/stdlib/integers.c
@@ -14,13 +14,6 @@
 #include "text.h"
 #include "types.h"
 
-static _Thread_local gmp_randstate_t Int_rng = {};
-
-public void Int$init_random(long seed) {
-    gmp_randinit_default(Int_rng);
-    gmp_randseed_ui(Int_rng, (unsigned long)seed);
-}
-
 public Text_t Int$value_as_text(Int_t i) {
     if (__builtin_expect(i.small & 1, 1)) {
         return Text$format("%ld", (i.small)>>2);
@@ -316,31 +309,6 @@ public Int_t Int$sqrt(Int_t i)
     return Int$from_mpz(result);
 }
 
-public Int_t Int$random(Int_t min, Int_t max) {
-    int32_t cmp = Int$compare_value(min, max);
-    if (cmp > 0) {
-        Text_t min_text = Int$as_text(&min, false, &Int$info), max_text = Int$as_text(&max, false, &Int$info);
-        fail("Random minimum value (%k) is larger than the maximum value (%k)",
-             &min_text, &max_text);
-    }
-    if (cmp == 0) return min;
-
-    mpz_t range_size;
-    mpz_init_set_int(range_size, max);
-    if (min.small & 1) {
-        mpz_t min_mpz;
-        mpz_init_set_si(min_mpz, min.small >> 2);
-        mpz_sub(range_size, range_size, min_mpz);
-    } else {
-        mpz_sub(range_size, range_size, *min.big);
-    }
-
-    mpz_t r;
-    mpz_init(r);
-    mpz_urandomm(r, Int_rng, range_size);
-    return Int$plus(min, Int$from_mpz(r));
-}
-
 public PUREFUNC Range_t Int$to(Int_t from, Int_t to) {
     return (Range_t){from, to, Int$compare_value(to, from) >= 0 ? (Int_t){.small=(1<<2)|1} : (Int_t){.small=(-1>>2)|1}};
 }
@@ -440,25 +408,6 @@ public const TypeInfo_t Int$info = {
         } \
         return bit_array; \
     } \
-    public c_type KindOfInt ## $full_random(void) { \
-        c_type r; \
-        arc4random_buf(&r, sizeof(r)); \
-        return r; \
-    } \
-    public c_type KindOfInt ## $random(c_type min, c_type max) { \
-        if (min > max) fail("Random minimum value (%ld) is larger than the maximum value (%ld)", min, max); \
-        if (min == max) return min; \
-        if (min == min_val && max == max_val) \
-            return KindOfInt ## $full_random(); \
-        uint64_t range = (uint64_t)max - (uint64_t)min + 1; \
-        uint64_t min_r = -range % range; \
-        uint64_t r; \
-        for (;;) { \
-            arc4random_buf(&r, sizeof(r)); \
-            if (r >= min_r) break; \
-        } \
-        return (c_type)((uint64_t)min + (r % range)); \
-    } \
     public to_attr Range_t KindOfInt ## $to(c_type from, c_type to) { \
         return (Range_t){Int64_to_Int(from), Int64_to_Int(to), to >= from ? (Int_t){.small=(1<<2)&1} : (Int_t){.small=(1<<2)&1}}; \
     } \
diff --git a/stdlib/integers.h b/stdlib/integers.h
index 04699162..b441e2be 100644
--- a/stdlib/integers.h
+++ b/stdlib/integers.h
@@ -34,8 +34,6 @@
     Text_t type_name ## $hex(c_type i, Int_t digits, bool uppercase, bool prefix); \
     Text_t type_name ## $octal(c_type i, Int_t digits, bool prefix); \
     Array_t type_name ## $bits(c_type x); \
-    c_type type_name ## $random(c_type min, c_type max); \
-    c_type type_name ## $full_random(void); \
     to_attr Range_t type_name ## $to(c_type from, c_type to); \
     PUREFUNC Optional ## type_name ## _t type_name ## $from_text(Text_t text); \
     MACROLIKE PUREFUNC c_type type_name ## $clamped(c_type x, c_type min, c_type max) { \
@@ -103,8 +101,6 @@ PUREFUNC bool Int$equal_value(const Int_t x, const Int_t y);
 Text_t Int$format(Int_t i, Int_t digits);
 Text_t Int$hex(Int_t i, Int_t digits, bool uppercase, bool prefix);
 Text_t Int$octal(Int_t i, Int_t digits, bool prefix);
-void Int$init_random(long seed);
-Int_t Int$random(Int_t min, Int_t max);
 PUREFUNC Range_t Int$to(Int_t from, Int_t to);
 OptionalInt_t Int$from_str(const char *str);
 OptionalInt_t Int$from_text(Text_t text);
@@ -127,7 +123,7 @@ Int_t Int$sqrt(Int_t i);
 } while (0)
 
 #define I(i) ((int64_t)(i) == (int32_t)(i) ? ((Int_t){.small=(int64_t)((uint64_t)(i)<<2)|1}) : Int64_to_Int(i))
-#define I_small(i) ((Int_t){.small=((uint64_t)(i)<<2)|1})
+#define I_small(i) ((Int_t){.small=(int64_t)((uint64_t)(i)<<2)|1})
 #define I_is_zero(i) ((i).small == 1)
 
 Int_t Int$slow_plus(Int_t x, Int_t y);
diff --git a/stdlib/nums.c b/stdlib/nums.c
index 1a8fec21..b8de553e 100644
--- a/stdlib/nums.c
+++ b/stdlib/nums.c
@@ -57,10 +57,6 @@ public CONSTFUNC double Num$mod(double num, double modulus) {
     return (result < 0) != (modulus < 0) ? result + modulus : result; 
 }
 
-public double Num$random(void) { 
-    return drand48(); 
-}
-
 public CONSTFUNC double Num$mix(double amount, double x, double y) { 
     return (1.0-amount)*x + amount*y;
 }
@@ -138,10 +134,6 @@ public CONSTFUNC float Num32$mod(float num, float modulus) {
     return (result < 0) != (modulus < 0) ? result + modulus : result; 
 }
 
-public float Num32$random(void) { 
-    return (float)drand48(); 
-}
-
 public CONSTFUNC float Num32$mix(float amount, float x, float y) { 
     return (1.0f-amount)*x + amount*y;
 }
diff --git a/stdlib/nums.h b/stdlib/nums.h
index fc6d804f..a07a2710 100644
--- a/stdlib/nums.h
+++ b/stdlib/nums.h
@@ -27,7 +27,6 @@ CONSTFUNC bool Num$isinf(double n);
 CONSTFUNC bool Num$finite(double n);
 CONSTFUNC bool Num$isnan(double n);
 double Num$nan(Text_t tag);
-double Num$random(void);
 CONSTFUNC double Num$mix(double amount, double x, double y);
 OptionalNum_t Num$from_text(Text_t text);
 MACROLIKE CONSTFUNC double Num$clamped(double x, double low, double high) {
@@ -45,7 +44,6 @@ float Num32$mod(float num, float modulus);
 CONSTFUNC bool Num32$isinf(float n);
 CONSTFUNC bool Num32$finite(float n);
 CONSTFUNC bool Num32$isnan(float n);
-float Num32$random(void);
 CONSTFUNC float Num32$mix(float amount, float x, float y);
 OptionalNum32_t Num32$from_text(Text_t text);
 float Num32$nan(Text_t tag);
diff --git a/stdlib/rng.c b/stdlib/rng.c
new file mode 100644
index 00000000..c69a2771
--- /dev/null
+++ b/stdlib/rng.c
@@ -0,0 +1,265 @@
+// Random Number Generator (RNG) implementation based on ChaCha
+
+#include <ctype.h>
+#include <err.h>
+#include <gc.h>
+#include <gmp.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <sys/param.h>
+
+#include "arrays.h"
+#include "datatypes.h"
+#include "rng.h"
+#include "text.h"
+#include "util.h"
+
+#include "chacha.h"
+
+public _Thread_local RNG_t default_rng;
+
+struct RNGState_t {
+    chacha_ctx chacha;
+    size_t unused_bytes;
+    uint8_t buf[16*64];
+};
+
+PUREFUNC static Text_t RNG$as_text(const RNG_t *rng, bool colorize, const TypeInfo_t *type)
+{
+    (void)type;
+    if (!rng) return Text("RNG");
+    return Text$format(colorize ? "\x1b[34;1mRNG(%p)\x1b[m" : "RNG(%p)", *rng);
+}
+
+#define KEYSZ 32
+#define IVSZ 8
+
+public void RNG$set_seed(RNG_t rng, Array_t seed)
+{
+    uint8_t seed_bytes[KEYSZ + IVSZ] = {};
+    for (int64_t i = 0; i < (int64_t)sizeof(seed_bytes); i++)
+        seed_bytes[i] = i < seed.length ? *(uint8_t*)(seed.data + i*seed.stride) : 0;
+
+    rng->unused_bytes = 0;
+    chacha_keysetup(&rng->chacha, seed_bytes, KEYSZ/8);
+    chacha_ivsetup(&rng->chacha, seed_bytes + KEYSZ);
+}
+
+public RNG_t RNG$copy(RNG_t rng)
+{
+    RNG_t copy = GC_MALLOC_ATOMIC(sizeof(struct RNGState_t));
+    *copy = *rng;
+    return copy;
+}
+
+public RNG_t RNG$new(Array_t seed)
+{
+    RNG_t rng = GC_MALLOC_ATOMIC(sizeof(struct RNGState_t));
+    RNG$set_seed(rng, seed);
+    return rng;
+}
+
+static void rekey(RNG_t rng)
+{
+    // Fill the buffer with the keystream
+    chacha_encrypt_bytes(&rng->chacha, rng->buf, rng->buf, sizeof(rng->buf));
+    // Immediately reinitialize for backtracking resistance
+    chacha_keysetup(&rng->chacha, rng->buf, KEYSZ/8);
+    chacha_ivsetup(&rng->chacha, rng->buf + KEYSZ);
+    memset(rng->buf, 0, KEYSZ + IVSZ);
+    rng->unused_bytes = sizeof(rng->buf) - KEYSZ - IVSZ;
+}
+
+static void random_bytes(RNG_t rng, uint8_t *dest, size_t needed)
+{
+    while (needed > 0) {
+        if (rng->unused_bytes > 0) {
+            size_t to_get = MIN(needed, rng->unused_bytes);
+            uint8_t *keystream = rng->buf + sizeof(rng->buf) - rng->unused_bytes;
+            memcpy(dest, keystream, to_get);
+            memset(keystream, 0, to_get);
+            dest += to_get;
+            needed -= to_get;
+            rng->unused_bytes -= to_get;
+        }
+        if (rng->unused_bytes == 0)
+            rekey(rng);
+    }
+}
+
+public Bool_t RNG$bool(RNG_t rng, Num_t p)
+{
+    if (p == 0.5) {
+        uint8_t b;
+        random_bytes(rng, &b, sizeof(b));
+        return b & 1;
+    } else {
+        return RNG$num(rng, 0.0, 1.0) < p;
+    }
+}
+
+public Int_t RNG$int(RNG_t rng, Int_t min, Int_t max)
+{
+    if (__builtin_expect(((min.small & max.small) & 1) != 0, 1)) {
+        int32_t r = RNG$int32(rng, (int32_t)(min.small >> 2), (int32_t)(max.small >> 2));
+        return I_small(r);
+    }
+
+    int32_t cmp = Int$compare_value(min, max);
+    if (cmp > 0) {
+        Text_t min_text = Int$as_text(&min, false, &Int$info), max_text = Int$as_text(&max, false, &Int$info);
+        fail("Random minimum value (%k) is larger than the maximum value (%k)",
+             &min_text, &max_text);
+    }
+    if (cmp == 0) return min;
+
+    mpz_t range_size;
+    mpz_init_set_int(range_size, max);
+    if (min.small & 1) {
+        mpz_t min_mpz;
+        mpz_init_set_si(min_mpz, min.small >> 2);
+        mpz_sub(range_size, range_size, min_mpz);
+    } else {
+        mpz_sub(range_size, range_size, *min.big);
+    }
+
+    gmp_randstate_t gmp_rng;
+    gmp_randinit_default(gmp_rng);
+    gmp_randseed_ui(gmp_rng, (unsigned long)RNG$int64(rng, INT64_MIN, INT64_MAX));
+
+    mpz_t r;
+    mpz_init(r);
+    mpz_urandomm(r, gmp_rng, range_size);
+
+    gmp_randclear(gmp_rng);
+    return Int$plus(min, Int$from_mpz(r));
+}
+
+public Int64_t RNG$int64(RNG_t rng, Int64_t min, Int64_t max)
+{
+    if (min > max) fail("Random minimum value (%ld) is larger than the maximum value (%ld)", min, max);
+    if (min == max) return min;
+    if (min == INT64_MIN && max == INT64_MAX) {
+        int64_t r;
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        return r;
+    }
+    uint64_t range = (uint64_t)max - (uint64_t)min + 1;
+    uint64_t min_r = -range % range;
+    uint64_t r;
+    for (;;) {
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        if (r >= min_r) break;
+    }
+    return (int64_t)((uint64_t)min + (r % range));
+}
+
+public Int32_t RNG$int32(RNG_t rng, Int32_t min, Int32_t max)
+{
+    if (min > max) fail("Random minimum value (%d) is larger than the maximum value (%d)", min, max);
+    if (min == max) return min;
+    if (min == INT32_MIN && max == INT32_MAX) {
+        int32_t r;
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        return r;
+    }
+    uint32_t range = (uint32_t)max - (uint32_t)min + 1;
+    uint32_t min_r = -range % range;
+    uint32_t r;
+    for (;;) {
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        if (r >= min_r) break;
+    }
+    return (int32_t)((uint32_t)min + (r % range));
+}
+
+public Int16_t RNG$int16(RNG_t rng, Int16_t min, Int16_t max)
+{
+    if (min > max) fail("Random minimum value (%d) is larger than the maximum value (%d)", min, max);
+    if (min == max) return min;
+    if (min == INT16_MIN && max == INT16_MAX) {
+        int16_t r;
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        return r;
+    }
+    uint16_t range = (uint16_t)max - (uint16_t)min + 1;
+    uint16_t min_r = -range % range;
+    uint16_t r;
+    for (;;) {
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        if (r >= min_r) break;
+    }
+    return (int16_t)((uint16_t)min + (r % range));
+}
+
+public Int8_t RNG$int8(RNG_t rng, Int8_t min, Int8_t max)
+{
+    if (min > max) fail("Random minimum value (%d) is larger than the maximum value (%d)", min, max);
+    if (min == max) return min;
+    if (min == INT8_MIN && max == INT8_MAX) {
+        int8_t r;
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        return r;
+    }
+    uint8_t range = (uint8_t)max - (uint8_t)min + 1;
+    uint8_t min_r = -range % range;
+    uint8_t r;
+    for (;;) {
+        random_bytes(rng, (uint8_t*)&r, sizeof(r));
+        if (r >= min_r) break;
+    }
+    return (int8_t)((uint8_t)min + (r % range));
+}
+
+public Num_t RNG$num(RNG_t rng, Num_t min, Num_t max)
+{
+    if (min > max) fail("Random minimum value (%g) is larger than the maximum value (%g)", min, max);
+    if (min == max) return min;
+
+    union {
+        Num_t num;
+        uint64_t bits;
+    } r, one = {.num=1.0};
+    random_bytes(rng, (void*)&r, sizeof(r));
+
+    // Set r.num to 1.<random-bits>
+    r.bits &= ~(0xFFFULL << 52);
+    r.bits |= (one.bits & (0xFFFULL << 52));
+
+    r.num -= 1.0;
+
+    if (min == 0.0 && max == 1.0)
+        return r.num;
+
+    return (1.0-r.num)*min + r.num*max;
+}
+
+public Num32_t RNG$num32(RNG_t rng, Num32_t min, Num32_t max)
+{
+    return (Num32_t)RNG$num(rng, (Num_t)min, (Num_t)max);
+}
+
+public Byte_t RNG$byte(RNG_t rng)
+{
+    Byte_t b;
+    random_bytes(rng, &b, sizeof(b));
+    return b;
+}
+
+public Array_t RNG$bytes(RNG_t rng, Int_t count)
+{
+    int64_t n = Int_to_Int64(count, false);
+    Byte_t *r = GC_MALLOC_ATOMIC(sizeof(Byte_t[n]));
+    random_bytes(rng, r, sizeof(Byte_t[n]));
+    return (Array_t){.data=r, .length=n, .stride=1, .atomic=1};
+}
+
+public const TypeInfo_t RNG$info = {
+    .size=sizeof(void*),
+    .align=__alignof__(void*),
+    .tag=CustomInfo,
+    .CustomInfo={.as_text=(void*)RNG$as_text},
+};
+
+// vim: ts=4 sw=0 et cino=L2,l1,(0,W4,m1,\:0
diff --git a/stdlib/rng.h b/stdlib/rng.h
new file mode 100644
index 00000000..5bc4794f
--- /dev/null
+++ b/stdlib/rng.h
@@ -0,0 +1,31 @@
+#pragma once
+
+// Random Number Generator (RNG) functions/type info
+
+#include <stdbool.h>
+#include <stdint.h>
+
+#include "datatypes.h"
+#include "types.h"
+#include "bools.h"
+#include "bytes.h"
+#include "util.h"
+
+RNG_t RNG$new(Array_t seed);
+void RNG$set_seed(RNG_t rng, Array_t seed);
+RNG_t RNG$copy(RNG_t rng);
+Bool_t RNG$bool(RNG_t rng, Num_t p);
+Int_t RNG$int(RNG_t rng, Int_t min, Int_t max);
+Int64_t RNG$int64(RNG_t rng, Int64_t min, Int64_t max);
+Int32_t RNG$int32(RNG_t rng, Int32_t min, Int32_t max);
+Int16_t RNG$int16(RNG_t rng, Int16_t min, Int16_t max);
+Int8_t RNG$int8(RNG_t rng, Int8_t min, Int8_t max);
+Byte_t RNG$byte(RNG_t rng);
+Array_t RNG$bytes(RNG_t rng, Int_t count);
+Num_t RNG$num(RNG_t rng, Num_t min, Num_t max);
+Num32_t RNG$num32(RNG_t rng, Num32_t min, Num32_t max);
+
+extern const TypeInfo_t RNG$info;
+extern _Thread_local RNG_t default_rng;
+
+// vim: ts=4 sw=0 et cino=L2,l1,(0,W4,m1,\:0
diff --git a/stdlib/stdlib.c b/stdlib/stdlib.c
index fb37fcfe..f2c6897c 100644
--- a/stdlib/stdlib.c
+++ b/stdlib/stdlib.c
@@ -2,6 +2,7 @@
 
 #include <errno.h>
 #include <execinfo.h>
+#include <fcntl.h>
 #include <gc.h>
 #include <stdbool.h>
 #include <stdint.h>
@@ -17,6 +18,7 @@
 #include "metamethods.h"
 #include "patterns.h"
 #include "paths.h"
+#include "rng.h"
 #include "siphash.h"
 #include "stdlib.h"
 #include "tables.h"
@@ -30,14 +32,15 @@ public void tomo_init(void)
    GC_INIT();
    USE_COLOR = getenv("COLOR") ? strcmp(getenv("COLOR"), "1") == 0 : isatty(STDOUT_FILENO);
    getrandom(TOMO_HASH_KEY, sizeof(TOMO_HASH_KEY), 0);
-   unsigned int seed;
-   getrandom(&seed, sizeof(seed), 0);
-   srand(seed);
-   srand48(seed);
-
-   long long_seed;
-   getrandom(&long_seed, sizeof(long_seed), 0);
-   Int$init_random(long_seed);
+
+   int rng_fd = open("/dev/urandom", O_RDONLY);
+   if (rng_fd < 0)
+       fail("Couldn't read from /dev/urandom");
+   uint8_t *random_bytes = GC_MALLOC_ATOMIC(40);
+   if (read(rng_fd, (void*)random_bytes, 40) < 40)
+       fail("Couldn't read from /dev/urandom");
+   Array_t rng_seed = {.length=40, .data=random_bytes, .stride=1, .atomic=1};
+   default_rng = RNG$new(rng_seed);
 
    if (register_printf_specifier('k', printf_text, printf_text_size))
        errx(1, "Couldn't set printf specifier");
diff --git a/stdlib/threads.c b/stdlib/threads.c
index beb6771e..3ae2980b 100644
--- a/stdlib/threads.c
+++ b/stdlib/threads.c
@@ -2,6 +2,7 @@
 
 #include <ctype.h>
 #include <err.h>
+#include <fcntl.h>
 #include <gc.h>
 #include <math.h>
 #include <stdbool.h>
@@ -13,6 +14,7 @@
 
 #include "arrays.h"
 #include "datatypes.h"
+#include "rng.h"
 #include "text.h"
 #include "threads.h"
 #include "types.h"
@@ -20,9 +22,10 @@
 
 static void *run_thread(Closure_t *closure)
 {
-    long seed;
-    getrandom(&seed, sizeof(seed), 0);
-    Int$init_random(seed);
+    uint8_t *random_bytes = GC_MALLOC_ATOMIC(40);
+    getrandom(random_bytes, 40, 0);
+    Array_t rng_seed = {.length=40, .data=random_bytes, .stride=1, .atomic=1};
+    default_rng = RNG$new(rng_seed);
     ((void(*)(void*))closure->fn)(closure->userdata);
     return NULL;
 }
diff --git a/stdlib/tomo.h b/stdlib/tomo.h
index 64a9979b..8b378c0a 100644
--- a/stdlib/tomo.h
+++ b/stdlib/tomo.h
@@ -24,6 +24,7 @@
 #include "patterns.h"
 #include "pointers.h"
 #include "ranges.h"
+#include "rng.h"
 #include "shell.h"
 #include "siphash.h"
 #include "tables.h"