Implement LCG pseudorandom number functions

src/header/stdlib/lcg48.rs

@@ -0,0 +1,43 @@
+//! Helper functions for pseudorandom number generation using LCG, see http://pubs.opengroup.org/onlinepubs/7908799/xsh/drand48.html
+
+use platform::types::*;
+
+/* The current element of the linear congruential generator's sequence. Any
+ * function that sets this variable must ensure that only the lower 48 bits get
+ * set. */
+pub static mut XI: u64 = 0;
+
+/* Multiplier and addend, which may be set through lcong48(). Default values as
+ * specified in POSIX. */
+pub static mut A: u64 = 0x5deece66d;
+pub static mut C: u16 = 0xb;
+
+/// Advances the linear congruential generator to the next element in its
+/// sequence.
+pub unsafe fn generator_step() {
+    /* 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. */
+    XI = A.wrapping_mul(XI).wrapping_add(u64::from(C)) & 0xffff_ffff_ffff;
+}
+
+/// Get a C `double` from a 48-bit integer (for `drand48()` and `erand48()`).
+pub fn x_to_float64(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.0f64
+}
+
+/// Get the high 31 bits of a 48-bit integer (for `lrand48()` and `nrand48()`).
+pub fn x_to_uint31(x: u64) -> c_long {
+    (x >> 17) as c_long
+}
+
+/// Get the high 32 bits, signed, of a 48-bit integer (for `mrand48()` and
+/// `jrand48()`).
+pub fn x_to_int32(x: u64) -> c_long {
+    // Cast via i32 to ensure we get the sign correct
+    (x >> 16) as i32 as c_long
+}

src/header/stdlib/mod.rs

@@ -20,6 +20,7 @@ use platform;
 use platform::types::*;
 use platform::{Pal, Sys};
 
+mod lcg48;
 mod sort;
 
 pub const EXIT_FAILURE: c_int = 1;
@@ -227,9 +228,10 @@ pub extern "C" fn div(numer: c_int, denom: c_int) -> div_t {
     }
 }
 
-// #[no_mangle]
-pub extern "C" fn drand48() -> c_double {
-    unimplemented!();
+#[no_mangle]
+pub unsafe extern "C" fn drand48() -> c_double {
+    lcg48::generator_step();
+    lcg48::x_to_float64(lcg48::XI)
 }
 
 // #[no_mangle]
@@ -410,9 +412,10 @@ pub extern "C" fn lldiv(numer: c_longlong, denom: c_longlong) -> lldiv_t {
     }
 }
 
-// #[no_mangle]
-pub extern "C" fn lrand48() -> c_long {
-    unimplemented!();
+#[no_mangle]
+pub unsafe extern "C" fn lrand48() -> c_long {
+    lcg48::generator_step();
+    lcg48::x_to_uint31(lcg48::XI)
 }
 
 #[no_mangle]
@@ -566,9 +569,10 @@ pub extern "C" fn mkstemps(name: *mut c_char, suffix_len: c_int) -> c_int {
     mkostemps(name, suffix_len, 0)
 }
 
-// #[no_mangle]
-pub extern "C" fn mrand48() -> c_long {
-    unimplemented!();
+#[no_mangle]
+pub unsafe extern "C" fn mrand48() -> c_long {
+    lcg48::generator_step();
+    lcg48::x_to_int32(lcg48::XI)
 }
 
 // #[no_mangle]
@@ -775,9 +779,12 @@ pub unsafe extern "C" fn srand(seed: c_uint) {
     RNG = Some(XorShiftRng::from_seed([seed as u8; 16]));
 }
 
-// #[no_mangle]
-pub extern "C" fn srand48(seed: c_long) {
-    unimplemented!();
+#[no_mangle]
+pub unsafe extern "C" fn srand48(seedval: c_long) {
+    /* 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. */
+    lcg48::XI = (((seedval & 0xffff_ffff) as u64) << 16) | 0x330e_u64;
 }
 
 // #[no_mangle]