From 45860e92565c115e395a96bf6042117c719c5b4b Mon Sep 17 00:00:00 2001 From: Peter Limkilde Svendsen <peter.limkilde@gmail.com> Date: Sun, 12 May 2019 14:50:18 +0000 Subject: [PATCH] Implement LCG pseudorandom number functions --- src/header/stdlib/lcg48.rs | 43 ++++++++++++++++++++++++++++ src/header/stdlib/mod.rs | 31 ++++++++++++-------- tests/Makefile | 1 + tests/expected/stdlib/lcg48.stderr | 0 tests/expected/stdlib/lcg48.stdout | 4 +++ tests/stdlib/lcg48.c | 45 ++++++++++++++++++++++++++++++ 6 files changed, 112 insertions(+), 12 deletions(-) create mode 100644 src/header/stdlib/lcg48.rs create mode 100644 tests/expected/stdlib/lcg48.stderr create mode 100644 tests/expected/stdlib/lcg48.stdout create mode 100644 tests/stdlib/lcg48.c diff --git a/src/header/stdlib/lcg48.rs b/src/header/stdlib/lcg48.rs new file mode 100644 index 000000000..cba8188ca --- /dev/null +++ b/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 +} diff --git a/src/header/stdlib/mod.rs b/src/header/stdlib/mod.rs index 57ed49aa9..fb7457824 100644 --- a/src/header/stdlib/mod.rs +++ b/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] diff --git a/tests/Makefile b/tests/Makefile index c051b92e5..d4e0e9a40 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -41,6 +41,7 @@ EXPECT_NAMES=\ stdlib/atoi \ stdlib/div \ stdlib/env \ + stdlib/lcg48 \ stdlib/mkostemps \ stdlib/rand \ stdlib/strtod \ diff --git a/tests/expected/stdlib/lcg48.stderr b/tests/expected/stdlib/lcg48.stderr new file mode 100644 index 000000000..e69de29bb diff --git a/tests/expected/stdlib/lcg48.stdout b/tests/expected/stdlib/lcg48.stdout new file mode 100644 index 000000000..14522ff60 --- /dev/null +++ b/tests/expected/stdlib/lcg48.stdout @@ -0,0 +1,4 @@ +lrand48 (uninitialized): 0 2116118 89401895 379337186 782977366 196130996 198207689 1046291021 1131187612 975888346 +drand48: 0.750266 0.607593 0.567593 0.799563 0.984984 0.205670 0.625922 0.794426 0.369416 0.854100 +lrand48: 1611183183 1304796356 1218897288 1717049088 2115236938 441672110 1344158015 1706017430 793314380 1834165927 +mrand48: -1072600929 -1685374584 -1857172720 -860869119 -64493420 883344220 -1606651266 -882932436 1586628760 -626635441 diff --git a/tests/stdlib/lcg48.c b/tests/stdlib/lcg48.c new file mode 100644 index 000000000..da8c76a46 --- /dev/null +++ b/tests/stdlib/lcg48.c @@ -0,0 +1,45 @@ +#include <stdlib.h> +#include <stdio.h> + +#include "test_helpers.h" + +int main(void) { + long x_l, x_m; + double x_d; + long seedval = 0xcafebeef; + + printf("lrand48 (uninitialized):"); + for (int i = 0; i < 10; i++) + { + x_l = lrand48(); + printf(" %ld", x_l); + } + printf("\n"); + + srand48(seedval); + printf("drand48:"); + for (int i = 0; i < 10; i++) + { + x_d = drand48(); + printf(" %lf", x_d); + } + printf("\n"); + + srand48(seedval); + printf("lrand48:"); + for (int i = 0; i < 10; i++) + { + x_l = lrand48(); + printf(" %ld", x_l); + } + printf("\n"); + + srand48(seedval); + printf("mrand48:"); + for (int i = 0; i < 10; i++) + { + x_m = mrand48(); + printf(" %ld", x_m); + } + printf("\n"); +} -- GitLab