diff --git a/src/header/stdlib/lcg48.rs b/src/header/stdlib/lcg48.rs index cba8188ca0f30e550d81343dde1c88b865220d4d..a425110b04cb5e6247e3c293f2ca72940332e81b 100644 --- a/src/header/stdlib/lcg48.rs +++ b/src/header/stdlib/lcg48.rs @@ -1,43 +1,86 @@ -//! Helper functions for pseudorandom number generation using LCG, see http://pubs.opengroup.org/onlinepubs/7908799/xsh/drand48.html +//! Helper functions for pseudorandom number generation using LCG, see https://pubs.opengroup.org/onlinepubs/9699919799.2018edition/functions/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; +/* 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 STASHED_XI (see below). */ +pub static mut DEFAULT_XI: [c_ushort; 3] = [0; 3]; -/* 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; +// Used by seed48() (returns a pointer to this array). +pub static mut STASHED_XI: [c_ushort; 3] = [0; 3]; + +/* 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; + +pub static mut A: u64 = A_DEFAULT_VALUE; +pub static mut C: u16 = C_DEFAULT_VALUE; + +/// Used by `srand48()` and `seed48()`. +pub unsafe fn reset_a_and_c() { + A = A_DEFAULT_VALUE; + C = C_DEFAULT_VALUE; +} + +/// Build a 48-bit integer from a size-3 array of unsigned short. +/// +/// Takes a pointer argument due to the inappropriate C function +/// signatures generated from Rust's sized arrays, see +/// https://github.com/eqrion/cbindgen/issues/171 +pub unsafe fn uint48_from_ushort_arr3(arr_ptr: *const c_ushort) -> u64 { + let arr = [*arr_ptr.offset(0), *arr_ptr.offset(1), *arr_ptr.offset(2)]; + + /* 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) +} + +/// Set a size-3 array of unsigned short from a 48-bit integer. +pub unsafe fn set_ushort_arr3_from_uint48(arr_ptr: *mut c_ushort, value: u64) { + *arr_ptr.offset(0) = c_ushort::from(value as u16); + *arr_ptr.offset(1) = c_ushort::from((value >> 16) as u16); + *arr_ptr.offset(2) = c_ushort::from((value >> 32) as u16); +} + +/// 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. The input argument must be a size-3 array. +pub unsafe fn generator_step(xi_arr_ptr: *mut c_ushort) -> u64 { + let old_xi: u64 = uint48_from_ushort_arr3(xi_arr_ptr); -/// 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; + * 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_xi: u64 = A.wrapping_mul(old_xi).wrapping_add(u64::from(C)) & 0xffff_ffff_ffff; + + set_ushort_arr3_from_uint48(xi_arr_ptr, new_xi); + new_xi } -/// Get a C `double` from a 48-bit integer (for `drand48()` and `erand48()`). -pub fn x_to_float64(x: u64) -> c_double { +/// Get a C `double` from a 48-bit integer (for `drand48()` and +/// `erand48()`). +pub fn float64_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.0f64 } -/// Get the high 31 bits of a 48-bit integer (for `lrand48()` and `nrand48()`). -pub fn x_to_uint31(x: u64) -> c_long { +/// Get the high 31 bits of a 48-bit integer (for `lrand48()` and +/// `nrand48()`). +pub fn uint31_from_x(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 { +/// Get the high 32 bits, signed, of a 48-bit integer (for `mrand48()` +/// and `jrand48()`). +pub fn int32_from_x(x: u64) -> c_long { // Cast via i32 to ensure we get the sign correct - (x >> 16) as i32 as c_long + c_long::from((x >> 16) as i32) } diff --git a/src/header/stdlib/mod.rs b/src/header/stdlib/mod.rs index b1170693b6bf85749d61da9923ecac17bce9ef24..380bca304290fc447555f5d0691e1ff8b885d853 100644 --- a/src/header/stdlib/mod.rs +++ b/src/header/stdlib/mod.rs @@ -232,8 +232,8 @@ pub extern "C" fn div(numer: c_int, denom: c_int) -> div_t { #[no_mangle] pub unsafe extern "C" fn drand48() -> c_double { - lcg48::generator_step(); - lcg48::x_to_float64(lcg48::XI) + let new_xi = lcg48::generator_step(lcg48::DEFAULT_XI.as_mut_ptr()); + lcg48::float64_from_x(new_xi) } // #[no_mangle] @@ -246,9 +246,10 @@ pub extern "C" fn ecvt( unimplemented!(); } -// #[no_mangle] -pub extern "C" fn erand(xsubi: [c_ushort; 3]) -> c_double { - unimplemented!(); +#[no_mangle] +pub unsafe extern "C" fn erand48(xsubi: *mut c_ushort) -> c_double { + let new_xi = lcg48::generator_step(xsubi); + lcg48::float64_from_x(new_xi) } #[no_mangle] @@ -361,9 +362,10 @@ pub extern "C" fn initstate(seec: c_uint, state: *mut c_char, size: size_t) -> * unimplemented!(); } -// #[no_mangle] -pub extern "C" fn jrand48(xsubi: [c_ushort; 3]) -> c_long { - unimplemented!(); +#[no_mangle] +pub unsafe extern "C" fn jrand48(xsubi: *mut c_ushort) -> c_long { + let new_xi = lcg48::generator_step(xsubi); + lcg48::int32_from_x(new_xi) } // #[no_mangle] @@ -376,9 +378,16 @@ pub extern "C" fn labs(i: c_long) -> c_long { i.abs() } -// #[no_mangle] -pub extern "C" fn lcong48(param: [c_ushort; 7]) { - unimplemented!(); +#[no_mangle] +pub unsafe extern "C" fn lcong48(param: *mut c_ushort) { + // Input should be a size-7 array. + + /* Go through this ptr -> u64 -> ptr conversion to ensure we only + * get the lower 16 bits of each element. */ + let new_xi = lcg48::uint48_from_ushort_arr3(param.offset(0)); + lcg48::set_ushort_arr3_from_uint48(lcg48::DEFAULT_XI.as_mut_ptr(), new_xi); + lcg48::A = lcg48::uint48_from_ushort_arr3(param.offset(3)); + lcg48::C = *param.offset(6) as u16; // c_ushort may be more than 16 bits } #[repr(C)] @@ -416,8 +425,8 @@ pub extern "C" fn lldiv(numer: c_longlong, denom: c_longlong) -> lldiv_t { #[no_mangle] pub unsafe extern "C" fn lrand48() -> c_long { - lcg48::generator_step(); - lcg48::x_to_uint31(lcg48::XI) + let new_xi = lcg48::generator_step(lcg48::DEFAULT_XI.as_mut_ptr()); + lcg48::uint31_from_x(new_xi) } #[no_mangle] @@ -573,13 +582,14 @@ pub extern "C" fn mkstemps(name: *mut c_char, suffix_len: c_int) -> c_int { #[no_mangle] pub unsafe extern "C" fn mrand48() -> c_long { - lcg48::generator_step(); - lcg48::x_to_int32(lcg48::XI) + let new_xi = lcg48::generator_step(lcg48::DEFAULT_XI.as_mut_ptr()); + lcg48::int32_from_x(new_xi) } -// #[no_mangle] -pub extern "C" fn nrand48(xsubi: [c_ushort; 3]) -> c_long { - unimplemented!(); +#[no_mangle] +pub unsafe extern "C" fn nrand48(xsubi: *mut c_ushort) -> c_long { + let new_xi = lcg48::generator_step(xsubi); + lcg48::uint31_from_x(new_xi) } #[no_mangle] @@ -719,9 +729,16 @@ pub unsafe extern "C" fn realpath(pathname: *const c_char, resolved: *mut c_char ptr } -// #[no_mangle] -pub extern "C" fn seed48(seed16v: [c_ushort; 3]) -> c_ushort { - unimplemented!(); +#[no_mangle] +pub unsafe extern "C" fn seed48(seed16v: *mut c_ushort) -> *mut c_ushort { + lcg48::reset_a_and_c(); + + lcg48::STASHED_XI = lcg48::DEFAULT_XI; + + let new_xi = lcg48::uint48_from_ushort_arr3(seed16v); + lcg48::set_ushort_arr3_from_uint48(lcg48::DEFAULT_XI.as_mut_ptr(), new_xi); + + lcg48::STASHED_XI.as_mut_ptr() } #[no_mangle] @@ -808,10 +825,13 @@ pub unsafe extern "C" fn srand(seed: c_uint) { #[no_mangle] pub unsafe extern "C" fn srand48(seedval: c_long) { + lcg48::reset_a_and_c(); + /* 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; + let new_xi = (((seedval as u32) as u64) << 16) | 0x330e_u64; + lcg48::set_ushort_arr3_from_uint48(lcg48::DEFAULT_XI.as_mut_ptr(), new_xi); } // #[no_mangle] diff --git a/tests/expected/stdlib/lcg48.stdout b/tests/expected/stdlib/lcg48.stdout index 14522ff60bc904344273f0d18ce229d04b20ca00..fa871030baa911d47a02c434e66f92db34c033f6 100644 --- a/tests/expected/stdlib/lcg48.stdout +++ b/tests/expected/stdlib/lcg48.stdout @@ -1,4 +1,14 @@ 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 +drand48 (seeded with srand48): 0.750266 0.607593 0.567593 0.799563 0.984984 0.205670 0.625922 0.794426 0.369416 0.854100 +lrand48 (seeded with srand48): 1611183183 1304796356 1218897288 1717049088 2115236938 441672110 1344158015 1706017430 793314380 1834165927 +mrand48 (seeded with srand48): -1072600929 -1685374584 -1857172720 -860869119 -64493420 883344220 -1606651266 -882932436 1586628760 -626635441 +erand48: 0.210555 0.014158 0.111353 0.658369 0.103767 0.180385 0.945033 0.745768 0.290272 0.111716 +nrand48: 514983590 590935818 1480794144 1496813112 2133865028 1816766485 2095074020 126058208 909762120 14734916 +jrand48: -1066398599 903693914 -1922375113 -2090140830 1218074962 1662411059 -722435322 764426686 -874142666 -1454656015 +seed48_return: [40c0, 4d4f, daa6] +lrand48 (seeded with seed48): 386486647 2049879217 706208537 1265096744 1586830881 1884641178 1566935266 1256810805 875501172 1670187935 +xsubi restored froom seed48 return value: [40c0, 4d4f, daa6] +nrand48 (from xsubi value returned by seed48): 1654466549 509433115 2007628085 1022767111 1935848687 967444157 1967758021 686622820 989863078 1688030308 +lrand48 (with parameters from lcong48): 2015972364 2009368981 971134301 317520085 149004773 538917235 242519436 1066970146 991527304 1588277058 +lrand48 (seeded with srand48 after lcong48 call): 1611183183 1304796356 1218897288 1717049088 2115236938 441672110 1344158015 1706017430 793314380 1834165927 +lrand48 (seeded with seed48 after lcong48 call): 386486647 2049879217 706208537 1265096744 1586830881 1884641178 1566935266 1256810805 875501172 1670187935 diff --git a/tests/stdlib/lcg48.c b/tests/stdlib/lcg48.c index da8c76a46ab55ca529348ab7ba03d81b4c830247..a1685f27fd207ec9f19615bd89fbb95bdfa85bc2 100644 --- a/tests/stdlib/lcg48.c +++ b/tests/stdlib/lcg48.c @@ -7,7 +7,12 @@ int main(void) { long x_l, x_m; double x_d; long seedval = 0xcafebeef; + unsigned short seed[3] = {0xfedc, 0xba98, 0x7654}; + unsigned short xsubi[3] = {0xabcd, 0xef42, 0x5678}; + unsigned short lcong48_params[7] = {0x0123, 0x4567, 0x89ab, 0xcdef, 0x4242, 0xf000, 0xbaaa}; + unsigned short xsubi_from_seed48[3] = {0, 0, 0}; + /* Test uninitialized behavior */ printf("lrand48 (uninitialized):"); for (int i = 0; i < 10; i++) { @@ -16,8 +21,10 @@ int main(void) { } printf("\n"); + /* Test different output types with same seed, builtin X_i and + * default multiplier and addend */ srand48(seedval); - printf("drand48:"); + printf("drand48 (seeded with srand48):"); for (int i = 0; i < 10; i++) { x_d = drand48(); @@ -26,7 +33,7 @@ int main(void) { printf("\n"); srand48(seedval); - printf("lrand48:"); + printf("lrand48 (seeded with srand48):"); for (int i = 0; i < 10; i++) { x_l = lrand48(); @@ -35,11 +42,96 @@ int main(void) { printf("\n"); srand48(seedval); - printf("mrand48:"); + printf("mrand48 (seeded with srand48):"); for (int i = 0; i < 10; i++) { x_m = mrand48(); printf(" %ld", x_m); } printf("\n"); + + /* Test corresponding functions taking user-supplied X_i, with + * default multiplier and addend */ + printf("erand48:"); + for (int i = 0; i < 10; i++) + { + x_d = erand48(xsubi); + printf(" %lf", x_d); + } + printf("\n"); + + printf("nrand48:"); + for (int i = 0; i < 10; i++) + { + x_l = nrand48(xsubi); + printf(" %ld", x_l); + } + printf("\n"); + + printf("jrand48:"); + for (int i = 0; i < 10; i++) + { + x_l = jrand48(xsubi); + printf(" %ld", x_l); + } + printf("\n"); + + /* Test seed48() "stashing" behavior. */ + unsigned short *seed48_return = seed48(seed); + printf("seed48_return: [%x, %x, %x]\n", + seed48_return[0], seed48_return[1], seed48_return[2]); + + /* Test seeding behavior of seed48() */ + printf("lrand48 (seeded with seed48):"); + for (int i = 0; i < 10; i++) + { + x_l = lrand48(); + printf(" %ld", x_l); + } + printf("\n"); + + /* Test restore from seed48()'s "stashed" value */ + xsubi_from_seed48[0] = seed48_return[0]; + xsubi_from_seed48[1] = seed48_return[1]; + xsubi_from_seed48[2] = seed48_return[2]; + printf("xsubi restored froom seed48 return value: [%x, %x, %x]\n", + xsubi_from_seed48[0], xsubi_from_seed48[1], xsubi_from_seed48[2]); + printf("nrand48 (from xsubi value returned by seed48):"); + for (int i = 0; i < 10; i++) + { + x_l = nrand48(xsubi_from_seed48); + printf(" %ld", x_l); + } + printf("\n"); + + /* Test behavior with all-user-defined parameters */ + lcong48(lcong48_params); + printf("lrand48 (with parameters from lcong48):"); + for (int i = 0; i < 10; i++) + { + x_l = lrand48(); + printf(" %ld", x_l); + } + printf("\n"); + + /* Test multiplier- and addend-restoring behavior of srand48() */ + srand48(seedval); + printf("lrand48 (seeded with srand48 after lcong48 call):"); + for (int i = 0; i < 10; i++) + { + x_l = lrand48(); + printf(" %ld", x_l); + } + printf("\n"); + + /* Test multiplier- and addend-restoring behavior of seed48() */ + lcong48(lcong48_params); + seed48(seed); + printf("lrand48 (seeded with seed48 after lcong48 call):"); + for (int i = 0; i < 10; i++) + { + x_l = lrand48(); + printf(" %ld", x_l); + } + printf("\n"); }