diff --git a/src/header/stdlib/mod.rs b/src/header/stdlib/mod.rs
index fa73a22adcf1de113cfc32b60c2831ae2bc337b3..fc65fba06ae780f9e213b6f26acd2cab227cf944 100644
--- a/src/header/stdlib/mod.rs
+++ b/src/header/stdlib/mod.rs
@@ -263,9 +263,11 @@ pub extern "C" fn div(numer: c_int, denom: c_int) -> div_t {
 }
 
 #[no_mangle]
-pub unsafe extern "C" fn drand48() -> c_double {
-    let new_xsubi_value = rand48::generator_step(&mut rand48::DEFAULT_XSUBI);
-    rand48::f64_from_x(new_xsubi_value)
+pub extern "C" fn drand48() -> c_double {
+    let params = rand48::params_lock();
+    let mut xsubi = rand48::xsubi_lock();
+    *xsubi = params.step(*xsubi);
+    xsubi.get_f64()
 }
 
 // #[no_mangle]
@@ -280,8 +282,11 @@ pub extern "C" fn ecvt(
 
 #[no_mangle]
 pub unsafe extern "C" fn erand48(xsubi: *mut c_ushort) -> c_double {
-    let new_xsubi_value = rand48::generator_step(&mut *(xsubi as *mut [c_ushort; 3]));
-    rand48::f64_from_x(new_xsubi_value)
+    let params = rand48::params_lock();
+    let xsubi_mut: &mut [c_ushort; 3] = slice::from_raw_parts_mut(xsubi, 3).try_into().unwrap();
+    let new_xsubi_value = params.step(xsubi_mut.into());
+    *xsubi_mut = new_xsubi_value.into();
+    new_xsubi_value.get_f64()
 }
 
 #[no_mangle]
@@ -440,8 +445,11 @@ pub unsafe extern "C" fn initstate(seed: c_uint, state: *mut c_char, size: size_
 
 #[no_mangle]
 pub unsafe extern "C" fn jrand48(xsubi: *mut c_ushort) -> c_long {
-    let new_xsubi_value = rand48::generator_step(&mut *(xsubi as *mut [c_ushort; 3]));
-    rand48::i32_from_x(new_xsubi_value)
+    let params = rand48::params_lock();
+    let xsubi_mut: &mut [c_ushort; 3] = slice::from_raw_parts_mut(xsubi, 3).try_into().unwrap();
+    let new_xsubi_value = params.step(xsubi_mut.into());
+    *xsubi_mut = new_xsubi_value.into();
+    new_xsubi_value.get_i32()
 }
 
 #[no_mangle]
@@ -490,16 +498,17 @@ pub extern "C" fn labs(i: c_long) -> c_long {
 
 #[no_mangle]
 pub unsafe extern "C" fn lcong48(param: *mut c_ushort) {
-    // Set DEFAULT_XSUBI buffer from elements 0-2
-    let xsubi_value = rand48::u48_from_ushort_arr3(&*(param as *const [c_ushort; 3]));
-    rand48::DEFAULT_XSUBI = rand48::ushort_arr3_from_u48(xsubi_value);
+    let mut xsubi = rand48::xsubi_lock();
+    let mut params = rand48::params_lock();
 
-    // Set multiplier from elements 3-5
-    rand48::A = rand48::u48_from_ushort_arr3(&*(param.offset(3) as *const [c_ushort; 3]));
+    let param_slice = slice::from_raw_parts(param, 7);
 
-    /* Set addend from element 6. Note that c_ushort may be more than 16
-     * bits, thus the cast. */
-    rand48::C = *param.offset(6) as u16;
+    let xsubi_ref: &[c_ushort; 3] = param_slice[0..3].try_into().unwrap();
+    let a_ref: &[c_ushort; 3] = param_slice[3..6].try_into().unwrap();
+    let c = param_slice[6];
+
+    *xsubi = xsubi_ref.into();
+    params.set(a_ref, c);
 }
 
 #[repr(C)]
@@ -536,9 +545,11 @@ pub extern "C" fn lldiv(numer: c_longlong, denom: c_longlong) -> lldiv_t {
 }
 
 #[no_mangle]
-pub unsafe extern "C" fn lrand48() -> c_long {
-    let new_xsubi_value = rand48::generator_step(&mut rand48::DEFAULT_XSUBI);
-    rand48::u31_from_x(new_xsubi_value)
+pub extern "C" fn lrand48() -> c_long {
+    let params = rand48::params_lock();
+    let mut xsubi = rand48::xsubi_lock();
+    *xsubi = params.step(*xsubi);
+    xsubi.get_u31()
 }
 
 #[no_mangle]
@@ -694,15 +705,20 @@ pub unsafe extern "C" fn mkstemps(name: *mut c_char, suffix_len: c_int) -> c_int
 }
 
 #[no_mangle]
-pub unsafe extern "C" fn mrand48() -> c_long {
-    let new_xsubi_value = rand48::generator_step(&mut rand48::DEFAULT_XSUBI);
-    rand48::i32_from_x(new_xsubi_value)
+pub extern "C" fn mrand48() -> c_long {
+    let params = rand48::params_lock();
+    let mut xsubi = rand48::xsubi_lock();
+    *xsubi = params.step(*xsubi);
+    xsubi.get_i32()
 }
 
 #[no_mangle]
 pub unsafe extern "C" fn nrand48(xsubi: *mut c_ushort) -> c_long {
-    let new_xsubi_value = rand48::generator_step(&mut *(xsubi as *mut [c_ushort; 3]));
-    rand48::u31_from_x(new_xsubi_value)
+    let params = rand48::params_lock();
+    let xsubi_mut: &mut [c_ushort; 3] = slice::from_raw_parts_mut(xsubi, 3).try_into().unwrap();
+    let new_xsubi_value = params.step(xsubi_mut.into());
+    *xsubi_mut = new_xsubi_value.into();
+    new_xsubi_value.get_u31()
 }
 
 #[no_mangle]
@@ -978,17 +994,17 @@ pub unsafe extern "C" fn realpath(pathname: *const c_char, resolved: *mut c_char
 
 #[no_mangle]
 pub unsafe extern "C" fn seed48(seed16v: *mut c_ushort) -> *mut c_ushort {
-    rand48::reset_a_and_c();
+    static mut BUFFER: [c_ushort; 3] = [0; 3];
 
-    // Stash current DEFAULT_XSUBI value in SEED48_XSUBI
-    rand48::SEED48_XSUBI = rand48::DEFAULT_XSUBI;
+    let mut params = rand48::params_lock();
+    let mut xsubi = rand48::xsubi_lock();
 
-    // Set DEFAULT_XSUBI from the argument provided
-    let xsubi_value = rand48::u48_from_ushort_arr3(&*(seed16v as *const [c_ushort; 3]));
-    rand48::DEFAULT_XSUBI = rand48::ushort_arr3_from_u48(xsubi_value);
+    let seed16v_ref: &[c_ushort; 3] = slice::from_raw_parts(seed16v, 3).try_into().unwrap();
 
-    // Return the stashed value
-    rand48::SEED48_XSUBI.as_mut_ptr()
+    BUFFER = (*xsubi).into();
+    *xsubi = seed16v_ref.into();
+    params.reset();
+    BUFFER.as_mut_ptr()
 }
 
 unsafe fn copy_kv(
@@ -1066,14 +1082,17 @@ pub unsafe extern "C" fn srand(seed: c_uint) {
 }
 
 #[no_mangle]
-pub unsafe extern "C" fn srand48(seedval: c_long) {
-    rand48::reset_a_and_c();
+pub extern "C" fn srand48(seedval: c_long) {
+    let mut params = rand48::params_lock();
+    let mut xsubi = rand48::xsubi_lock();
 
+    params.reset();
     /* Set the high 32 bits of the 48-bit X_i value to the lower 32 bits
      * of the input argument, and the lower 16 bits to 0x330e, as
      * specified in POSIX. */
-    let xsubi_value = (u64::from(seedval as u32) << 16) | 0x330e;
-    rand48::DEFAULT_XSUBI = rand48::ushort_arr3_from_u48(xsubi_value);
+    *xsubi = ((u64::from(seedval as u32) << 16) | 0x330e)
+        .try_into()
+        .unwrap();
 }
 
 #[no_mangle]
diff --git a/src/header/stdlib/rand48.rs b/src/header/stdlib/rand48.rs
index 57e968ec5c370cdba76f108972c02c2796dbabf7..6d78d6c0cc2a0a4598a6a67340cd9cec879d4447 100644
--- a/src/header/stdlib/rand48.rs
+++ b/src/header/stdlib/rand48.rs
@@ -1,90 +1,137 @@
 //! Helper functions for pseudorandom number generation using LCG, see https://pubs.opengroup.org/onlinepubs/9699919799.2018edition/functions/drand48.html
 
-use crate::platform::types::*;
+use crate::{
+    platform::types::*,
+    sync::{Mutex, MutexGuard},
+};
 
-/* The default element buffer for the linear congruential generator's
- * sequence. Implemented using a c_ushort array for consistency between
- * the drand48()/lrand48()/mrand48() and erand48()/nrand48()/jrand48()
- * functions, and with SEED48_XSUBI (see below). */
-pub static mut DEFAULT_XSUBI: [c_ushort; 3] = [0; 3];
+/// A 48-bit integer, used for the 48-bit arithmetic in these functions.
+#[derive(Clone, Copy, Default)]
+#[repr(C)]
+pub struct U48(u64);
 
-// Used by seed48() (returns a pointer to this array).
-pub static mut SEED48_XSUBI: [c_ushort; 3] = [0; 3];
+impl From<&[c_ushort; 3]> for U48 {
+    fn from(value: &[c_ushort; 3]) -> Self {
+        /* Cast via u16 to ensure we get only the lower 16 bits of each
+         * element, as specified by POSIX. */
+        Self {
+            0: u64::from(value[0] as u16)
+                | (u64::from(value[1] as u16) << 16)
+                | (u64::from(value[2] as u16) << 32),
+        }
+    }
+}
 
-/* Multiplier and addend, which may be set through lcong48(). Default
- * values as specified in POSIX. */
-const A_DEFAULT_VALUE: u64 = 0x5deece66d;
-const C_DEFAULT_VALUE: u16 = 0xb;
+impl From<&mut [c_ushort; 3]> for U48 {
+    fn from(value: &mut [c_ushort; 3]) -> Self {
+        Self::from(&*value)
+    }
+}
 
-pub static mut A: u64 = A_DEFAULT_VALUE;
-pub static mut C: u16 = C_DEFAULT_VALUE;
+impl TryFrom<u64> for U48 {
+    type Error = u64;
 
-/// Used by `srand48()` and `seed48()`.
-pub unsafe fn reset_a_and_c() {
-    A = A_DEFAULT_VALUE;
-    C = C_DEFAULT_VALUE;
+    fn try_from(value: u64) -> Result<Self, u64> {
+        if value < 0x1_0000_0000_0000 {
+            Ok(Self { 0: value })
+        } else {
+            Err(value)
+        }
+    }
 }
 
-/// Build a 48-bit integer from a size-3 array of unsigned short.
-///
-/// Pointers to c_ushort can be converted to &[c_ushort; 3] by taking
-/// &*(YOUR_C_USHORT_POINTER as *const [c_ushort; 3])
-///
-/// See also this cbindgen issue for why the stdlib functions can't just
-/// have an xsubi: *mut [c_ushort; 3] parameter:
-/// https://github.com/eqrion/cbindgen/issues/171
-pub fn u48_from_ushort_arr3(arr: &[c_ushort; 3]) -> u64 {
-    /* Cast via u16 to ensure we get only the lower 16 bits of each
-     * element, as specified by POSIX. */
-    u64::from(arr[0] as u16) | (u64::from(arr[1] as u16) << 16) | (u64::from(arr[2] as u16) << 32)
+impl From<U48> for u64 {
+    fn from(value: U48) -> Self {
+        value.0
+    }
 }
 
-/// Make a size-3 array of unsigned short from a 48-bit integer.
-pub fn ushort_arr3_from_u48(value: u64) -> [c_ushort; 3] {
-    [
-        c_ushort::from(value as u16),
-        c_ushort::from((value >> 16) as u16),
-        c_ushort::from((value >> 32) as u16),
-    ]
+impl From<U48> for [c_ushort; 3] {
+    fn from(value: U48) -> Self {
+        [
+            // "as u16" in case c_ushort is larger than u16
+            (value.0 as u16).into(),
+            ((value.0 >> 16) as u16).into(),
+            ((value.0 >> 32) as u16).into(),
+        ]
+    }
 }
 
-/// Advances the buffer from the input argument to the next element in
-/// the linear congruential generator's sequence.
-///
-/// Modifies the passed argument in-place and returns the new value as a
-/// u64.
-pub unsafe fn generator_step(xsubi: &mut [c_ushort; 3]) -> u64 {
-    let old_xsubi_value: u64 = u48_from_ushort_arr3(xsubi);
-
-    /* The recurrence relation of the linear congruential generator,
-     * X_(n+1) = (a * X_n + c) % m,
-     * with m = 2**48. The multiplication and addition can overflow a
-     * u64, but we just let it wrap since we take mod 2**48 anyway. */
-    let new_xsubi_value: u64 =
-        A.wrapping_mul(old_xsubi_value).wrapping_add(u64::from(C)) & 0xffff_ffff_ffff;
-
-    *xsubi = ushort_arr3_from_u48(new_xsubi_value);
-    new_xsubi_value
+impl U48 {
+    /// Get a C `double` in the interval [0.0, 1.0) (for `drand48()` and `erand48()`).
+    pub fn get_f64(self) -> c_double {
+        /* We set the exponent to 0, and the 48-bit integer is copied into the high
+         * 48 of the 52 significand bits. The value then lies in the range
+         * [1.0, 2.0), from which we simply subtract 1.0. */
+        f64::from_bits(0x3ff0_0000_0000_0000_u64 | (self.0 << 4)) - 1.0
+    }
+
+    /// Get the high 31 bits (for `lrand48()` and `nrand48()`).
+    pub fn get_u31(self) -> c_long {
+        (self.0 >> 17).try_into().unwrap()
+    }
+
+    /// Get the high 32 bits, signed (for `mrand48()` and `jrand48()`).
+    pub fn get_i32(self) -> c_long {
+        // Cast via i32 to ensure we get the sign correct
+        ((self.0 >> 16) as i32).into()
+    }
 }
 
-/// Get a C `double` from a 48-bit integer (for `drand48()` and
-/// `erand48()`).
-pub fn f64_from_x(x: u64) -> c_double {
-    /* We set the exponent to 0, and the 48-bit integer is copied into the high
-     * 48 of the 52 significand bits. The value then lies in the range
-     * [1.0, 2.0), from which we simply subtract 1.0. */
-    f64::from_bits(0x3ff0_0000_0000_0000_u64 | (x << 4)) - 1.0
+/// The a and c parameters of an LCG.
+#[derive(Default)]
+#[repr(C)]
+pub struct Params {
+    pub a: U48,
+    pub c: u16,
 }
 
-/// Get the high 31 bits of a 48-bit integer (for `lrand48()` and
-/// `nrand48()`).
-pub fn u31_from_x(x: u64) -> c_long {
-    (x >> 17) as c_long
+impl Params {
+    pub const fn new() -> Self {
+        // Default values as specified in POSIX
+        Params {
+            a: U48(0x5deece66d),
+            c: 0xb,
+        }
+    }
+
+    pub fn reset(&mut self) {
+        *self = Self::new();
+    }
+
+    /// For use in lcong48().
+    pub fn set(&mut self, a: &[c_ushort; 3], c: c_ushort) {
+        self.a = a.into();
+        self.c = c as u16; // Per POSIX, discard higher bits in case unsigned short is larger than u16
+    }
+
+    pub fn step(&self, xsubi: U48) -> U48 {
+        /* The recurrence relation of the linear congruential generator,
+         * X_(n+1) = (a * X_n + c) % m,
+         * with m = 2**48. The multiplication and addition can overflow a u64, but
+         * we just let it wrap since we take mod 2**48 anyway. */
+        (u64::from(self.a)
+            .wrapping_mul(u64::from(xsubi))
+            .wrapping_add(u64::from(self.c))
+            & 0xffff_ffff_ffff)
+            .try_into()
+            .unwrap()
+    }
 }
 
-/// Get the high 32 bits, signed, of a 48-bit integer (for `mrand48()`
-/// and `jrand48()`).
-pub fn i32_from_x(x: u64) -> c_long {
-    // Cast via i32 to ensure we get the sign correct
-    c_long::from((x >> 16) as i32)
+// TODO: consider using rwlock instead of mutex for more fine-grained access
+/// Immediately get the global Params lock, or panic if unsuccessful.
+pub fn params_lock<'a>() -> MutexGuard<'a, Params> {
+    static PARAMS: Mutex<Params> = Mutex::<Params>::new(Params::new());
+
+    PARAMS
+        .try_lock()
+        .expect("unable to acquire LCG parameter lock")
+}
+
+/// Immediately get the global X_i lock, or panic if unsuccessful.
+pub fn xsubi_lock<'a>() -> MutexGuard<'a, U48> {
+    static XSUBI: Mutex<U48> = Mutex::<U48>::new(U48(0));
+
+    XSUBI.try_lock().expect("unable to acquire LCG X_i lock")
 }
diff --git a/tests/test_helpers.h b/tests/test_helpers.h
index 2f9ed80cbae1ff1c6d4a1b9ad220ecb2fe7ea3be..6f5b3f01d3df99996cb58f772db8a29264e92a8b 100644
--- a/tests/test_helpers.h
+++ b/tests/test_helpers.h
@@ -98,7 +98,12 @@
         _exit(code); \
     } while(0)
 
-#define random_bool() (lrand48() % 2 == 0)
+// Duplicate of lrand48() logic but suitable for multithreaded use
+int random_bool() {
+    _Thread_local static uint64_t xsubi = 0;
+    xsubi = 0x5deece66d * xsubi + 0xb;
+    return (xsubi >> 17) % 2 == 0;
+}
 
 // Quick helper for checking desired errno status.
 // Use as macro: CHECK_AND_PRINT_ERRNO(<desired errno, e.g. EINVAL>);