diff --git a/TODO.md b/TODO.md
index 31eab53cff2d7c6970213f6789ef84bf61823bed..1ee2132cd54cdfb02ace914be0f3987d5c38417b 100644
--- a/TODO.md
+++ b/TODO.md
@@ -6,5 +6,3 @@
 
 * Define a combinator [smap] where the user function [f] has access to the
   size of the data.
-
-* Fix [random_bigint] so that it works above 2^30.
diff --git a/src/IFSeq.ml b/src/IFSeq.ml
index 14a4936a2d7cfe683a36e1600fa417e775cd9e7f..ea721411d8ecdd916adce1d52b9d8cafb1a559e5 100644
--- a/src/IFSeq.ml
+++ b/src/IFSeq.ml
@@ -8,27 +8,18 @@ include IFSeqSyn.Make(Z)
 let bigsum ss =
   List.fold_left sum zero ss
 
-(* For some reason, [Z.random] does not exist. *)
-(* For some reason, [Random.int] stops working at [2^30]. *)
-(* This is troublesome. *)
-
-let random_bigint (n : Z.t) : Z.t =
-  let open Z in
-  if n < one lsl 30 then
-    Z.of_int (Random.int (Z.to_int n))
-  else
-    failwith "Can't sample over more than 2^30 elements."
-      (* TEMPORARY really unsatisfactory! *)
-
 (* Extract a randomly chosen sample of [m] elements out of a sequence [s]. *)
 
 (* We do not protect against repetitions, as they are unlikely when [s] is
    long. *)
 
+(* For some reason, [Z.random] does not exist, and [Random.int] stops
+   working at [2^30]. [Bigint.random] works around these problems. *)
+
 let rec sample (m : int) (s : 'a seq) (k : 'a Seq.t) : 'a Seq.t =
   if m > 0 then
     fun () ->
-      let i = random_bigint (length s) in
+      let i = Bigint.random (length s) in
       let x = get s i in
       Seq.Cons (x, sample (m - 1) s k)
   else
diff --git a/src/bigint.ml b/src/bigint.ml
new file mode 100644
index 0000000000000000000000000000000000000000..52abbe8f3647d2613389f5579ada7dbdaa27b6ab
--- /dev/null
+++ b/src/bigint.ml
@@ -0,0 +1,62 @@
+(* Uniform random generation of large integers. Copied and adapted from Jane
+   Street's bignum library. *)
+
+(* [random state range] chooses a [depth] and generates random values using
+   [Random.State.bits state], called [1 lsl depth] times and concatenated. The
+   preliminary result [n] therefore satisfies [0 <= n < 1 lsl (30 lsl depth)].
+
+   In order for the random choice to be uniform between [0] and [range-1],
+   there must exist [k > 0] such that [n < k * range <= 1 lsl (30 lsl depth)].
+   If so, [n % range] is returned. Otherwise the random choice process is
+   repeated from scratch.
+
+   The [depth] value is chosen so that repeating is uncommon (1 in 1000 or
+   less). *)
+
+let bits_at_depth (depth : int) : int =
+  30 lsl depth
+
+let range_at_depth (depth : int) : Z.t =
+  Z.(one lsl (bits_at_depth depth))
+
+let rec choose_bit_depth_for_range_from range depth =
+  if Z.geq (range_at_depth depth) range then depth
+  else choose_bit_depth_for_range_from range (depth + 1)
+
+let choose_bit_depth_for_range (range : Z.t) : int =
+  choose_bit_depth_for_range_from range 0
+
+let rec random_bigint_at_depth (state : Random.State.t) depth : Z.t =
+  if depth = 0 then
+    Z.of_int (Random.State.bits state)
+  else
+    let depth = depth - 1 in
+    let prefix = random_bigint_at_depth state depth in
+    let suffix = random_bigint_at_depth state depth in
+    Z.(prefix lsl (bits_at_depth depth) lor suffix)
+
+let random_value_is_uniform_in_range range depth n =
+  let k = Z.(range_at_depth depth / range) in
+  Z.lt n Z.(k * range)
+
+let rec large_random_at_depth state range depth =
+  let result = random_bigint_at_depth state depth in
+  if random_value_is_uniform_in_range range depth result
+  then Z.(result mod range)
+  else large_random_at_depth state range depth
+
+let large_random state range =
+  let tolerance_factor = Z.of_int 1000 in
+  let depth = choose_bit_depth_for_range Z.(range * tolerance_factor) in
+  large_random_at_depth state range depth
+
+let random state range =
+  if Z.leq range Z.zero then
+    failwith (Printf.sprintf "Bigint.random: argument %s <= 0" (Z.to_string range))
+  else if Z.lt range Z.(one lsl 30) then
+    Z.of_int (Random.State.int state (Z.to_int range))
+  else
+    large_random state range
+
+let random range =
+  random (Random.get_state()) range
diff --git a/src/bigint.mli b/src/bigint.mli
new file mode 100644
index 0000000000000000000000000000000000000000..b90a2e60887d62713f57127857aa1fd92a67c407
--- /dev/null
+++ b/src/bigint.mli
@@ -0,0 +1,3 @@
+(* Uniform random generation of large integers. *)
+
+val random: Z.t -> Z.t