Skip to content
Snippets Groups Projects
Commit 3b12a400 authored by Jeremy Soller's avatar Jeremy Soller
Browse files

Merge branch 'lcg48_refactor' into 'master'

lcg48 refactor

See merge request !243
parents cc4c3a5c 7fdd450e
No related branches found
No related tags found
No related merge requests found
...@@ -5,11 +5,11 @@ use crate::platform::types::*; ...@@ -5,11 +5,11 @@ use crate::platform::types::*;
/* The default element buffer for the linear congruential generator's /* The default element buffer for the linear congruential generator's
* sequence. Implemented using a c_ushort array for consistency between * sequence. Implemented using a c_ushort array for consistency between
* the drand48()/lrand48()/mrand48() and erand48()/nrand48()/jrand48() * the drand48()/lrand48()/mrand48() and erand48()/nrand48()/jrand48()
* functions, and with STASHED_XI (see below). */ * functions, and with SEED48_XSUBI (see below). */
pub static mut DEFAULT_XI: [c_ushort; 3] = [0; 3]; pub static mut DEFAULT_XSUBI: [c_ushort; 3] = [0; 3];
// Used by seed48() (returns a pointer to this array). // Used by seed48() (returns a pointer to this array).
pub static mut STASHED_XI: [c_ushort; 3] = [0; 3]; pub static mut SEED48_XSUBI: [c_ushort; 3] = [0; 3];
/* Multiplier and addend, which may be set through lcong48(). Default /* Multiplier and addend, which may be set through lcong48(). Default
* values as specified in POSIX. */ * values as specified in POSIX. */
...@@ -27,60 +27,64 @@ pub unsafe fn reset_a_and_c() { ...@@ -27,60 +27,64 @@ pub unsafe fn reset_a_and_c() {
/// Build a 48-bit integer from a size-3 array of unsigned short. /// Build a 48-bit integer from a size-3 array of unsigned short.
/// ///
/// Takes a pointer argument due to the inappropriate C function /// Pointers to c_ushort can be converted to &[c_ushort; 3] by taking
/// signatures generated from Rust's sized arrays, see /// &*(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 /// https://github.com/eqrion/cbindgen/issues/171
pub unsafe fn uint48_from_ushort_arr3(arr_ptr: *const c_ushort) -> u64 { pub fn u48_from_ushort_arr3(arr: &[c_ushort; 3]) -> 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 /* Cast via u16 to ensure we get only the lower 16 bits of each
* element, as specified by POSIX. */ * 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) 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. /// Make 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) { pub fn ushort_arr3_from_u48(value: u64) -> [c_ushort; 3] {
*arr_ptr.offset(0) = c_ushort::from(value as u16); [
*arr_ptr.offset(1) = c_ushort::from((value >> 16) as u16); c_ushort::from(value as u16),
*arr_ptr.offset(2) = c_ushort::from((value >> 32) as u16); c_ushort::from((value >> 16) as u16),
c_ushort::from((value >> 32) as u16),
]
} }
/// Advances the buffer from the input argument to the next element in /// Advances the buffer from the input argument to the next element in
/// the linear congruential generator's sequence. /// the linear congruential generator's sequence.
/// ///
/// Modifies the passed argument in-place and returns the new value as a /// Modifies the passed argument in-place and returns the new value as a
/// u64. The input argument must be a size-3 array. /// u64.
pub unsafe fn generator_step(xi_arr_ptr: *mut c_ushort) -> u64 { pub unsafe fn generator_step(xsubi: &mut [c_ushort; 3]) -> u64 {
let old_xi: u64 = uint48_from_ushort_arr3(xi_arr_ptr); let old_xsubi_value: u64 = u48_from_ushort_arr3(xsubi);
/* The recurrence relation of the linear congruential generator, /* The recurrence relation of the linear congruential generator,
* X_(n+1) = (a * X_n + c) % m, * X_(n+1) = (a * X_n + c) % m,
* with m = 2**48. The multiplication and addition can overflow a * 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, 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; let new_xsubi_value: u64 =
A.wrapping_mul(old_xsubi_value).wrapping_add(u64::from(C)) & 0xffff_ffff_ffff;
set_ushort_arr3_from_uint48(xi_arr_ptr, new_xi); *xsubi = ushort_arr3_from_u48(new_xsubi_value);
new_xi new_xsubi_value
} }
/// Get a C `double` from a 48-bit integer (for `drand48()` and /// Get a C `double` from a 48-bit integer (for `drand48()` and
/// `erand48()`). /// `erand48()`).
pub fn float64_from_x(x: u64) -> c_double { 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 /* 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 * 48 of the 52 significand bits. The value then lies in the range
* [1.0, 2.0), from which we simply subtract 1.0. */ * [1.0, 2.0), from which we simply subtract 1.0. */
f64::from_bits(0x3ff0_0000_0000_0000_u64 | (x << 4)) - 1.0f64 f64::from_bits(0x3ff0_0000_0000_0000_u64 | (x << 4)) - 1.0
} }
/// Get the high 31 bits of a 48-bit integer (for `lrand48()` and /// Get the high 31 bits of a 48-bit integer (for `lrand48()` and
/// `nrand48()`). /// `nrand48()`).
pub fn uint31_from_x(x: u64) -> c_long { pub fn u31_from_x(x: u64) -> c_long {
(x >> 17) as c_long (x >> 17) as c_long
} }
/// Get the high 32 bits, signed, of a 48-bit integer (for `mrand48()` /// Get the high 32 bits, signed, of a 48-bit integer (for `mrand48()`
/// and `jrand48()`). /// and `jrand48()`).
pub fn int32_from_x(x: u64) -> c_long { pub fn i32_from_x(x: u64) -> c_long {
// Cast via i32 to ensure we get the sign correct // Cast via i32 to ensure we get the sign correct
c_long::from((x >> 16) as i32) c_long::from((x >> 16) as i32)
} }
...@@ -244,8 +244,8 @@ pub extern "C" fn div(numer: c_int, denom: c_int) -> div_t { ...@@ -244,8 +244,8 @@ pub extern "C" fn div(numer: c_int, denom: c_int) -> div_t {
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn drand48() -> c_double { pub unsafe extern "C" fn drand48() -> c_double {
let new_xi = lcg48::generator_step(lcg48::DEFAULT_XI.as_mut_ptr()); let new_xsubi_value = lcg48::generator_step(&mut lcg48::DEFAULT_XSUBI);
lcg48::float64_from_x(new_xi) lcg48::f64_from_x(new_xsubi_value)
} }
// #[no_mangle] // #[no_mangle]
...@@ -260,8 +260,8 @@ pub extern "C" fn ecvt( ...@@ -260,8 +260,8 @@ pub extern "C" fn ecvt(
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn erand48(xsubi: *mut c_ushort) -> c_double { pub unsafe extern "C" fn erand48(xsubi: *mut c_ushort) -> c_double {
let new_xi = lcg48::generator_step(xsubi); let new_xsubi_value = lcg48::generator_step(&mut *(xsubi as *mut [c_ushort; 3]));
lcg48::float64_from_x(new_xi) lcg48::f64_from_x(new_xsubi_value)
} }
#[no_mangle] #[no_mangle]
...@@ -376,8 +376,8 @@ pub extern "C" fn initstate(seec: c_uint, state: *mut c_char, size: size_t) -> * ...@@ -376,8 +376,8 @@ pub extern "C" fn initstate(seec: c_uint, state: *mut c_char, size: size_t) -> *
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn jrand48(xsubi: *mut c_ushort) -> c_long { pub unsafe extern "C" fn jrand48(xsubi: *mut c_ushort) -> c_long {
let new_xi = lcg48::generator_step(xsubi); let new_xsubi_value = lcg48::generator_step(&mut *(xsubi as *mut [c_ushort; 3]));
lcg48::int32_from_x(new_xi) lcg48::i32_from_x(new_xsubi_value)
} }
#[no_mangle] #[no_mangle]
...@@ -426,14 +426,16 @@ pub extern "C" fn labs(i: c_long) -> c_long { ...@@ -426,14 +426,16 @@ pub extern "C" fn labs(i: c_long) -> c_long {
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn lcong48(param: *mut c_ushort) { pub unsafe extern "C" fn lcong48(param: *mut c_ushort) {
// Input should be a size-7 array. // Set DEFAULT_XSUBI buffer from elements 0-2
let xsubi_value = lcg48::u48_from_ushort_arr3(&*(param as *const [c_ushort; 3]));
lcg48::DEFAULT_XSUBI = lcg48::ushort_arr3_from_u48(xsubi_value);
/* Go through this ptr -> u64 -> ptr conversion to ensure we only // Set multiplier from elements 3-5
* get the lower 16 bits of each element. */ lcg48::A = lcg48::u48_from_ushort_arr3(&*(param.offset(3) as *const [c_ushort; 3]));
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); /* Set addend from element 6. Note that c_ushort may be more than 16
lcg48::A = lcg48::uint48_from_ushort_arr3(param.offset(3)); * bits, thus the cast. */
lcg48::C = *param.offset(6) as u16; // c_ushort may be more than 16 bits lcg48::C = *param.offset(6) as u16;
} }
#[repr(C)] #[repr(C)]
...@@ -471,8 +473,8 @@ pub extern "C" fn lldiv(numer: c_longlong, denom: c_longlong) -> lldiv_t { ...@@ -471,8 +473,8 @@ pub extern "C" fn lldiv(numer: c_longlong, denom: c_longlong) -> lldiv_t {
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn lrand48() -> c_long { pub unsafe extern "C" fn lrand48() -> c_long {
let new_xi = lcg48::generator_step(lcg48::DEFAULT_XI.as_mut_ptr()); let new_xsubi_value = lcg48::generator_step(&mut lcg48::DEFAULT_XSUBI);
lcg48::uint31_from_x(new_xi) lcg48::u31_from_x(new_xsubi_value)
} }
#[no_mangle] #[no_mangle]
...@@ -625,14 +627,14 @@ pub extern "C" fn mkstemps(name: *mut c_char, suffix_len: c_int) -> c_int { ...@@ -625,14 +627,14 @@ pub extern "C" fn mkstemps(name: *mut c_char, suffix_len: c_int) -> c_int {
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn mrand48() -> c_long { pub unsafe extern "C" fn mrand48() -> c_long {
let new_xi = lcg48::generator_step(lcg48::DEFAULT_XI.as_mut_ptr()); let new_xsubi_value = lcg48::generator_step(&mut lcg48::DEFAULT_XSUBI);
lcg48::int32_from_x(new_xi) lcg48::i32_from_x(new_xsubi_value)
} }
#[no_mangle] #[no_mangle]
pub unsafe extern "C" fn nrand48(xsubi: *mut c_ushort) -> c_long { pub unsafe extern "C" fn nrand48(xsubi: *mut c_ushort) -> c_long {
let new_xi = lcg48::generator_step(xsubi); let new_xsubi_value = lcg48::generator_step(&mut *(xsubi as *mut [c_ushort; 3]));
lcg48::uint31_from_x(new_xi) lcg48::u31_from_x(new_xsubi_value)
} }
#[no_mangle] #[no_mangle]
...@@ -776,12 +778,15 @@ pub unsafe extern "C" fn realpath(pathname: *const c_char, resolved: *mut c_char ...@@ -776,12 +778,15 @@ pub unsafe extern "C" fn realpath(pathname: *const c_char, resolved: *mut c_char
pub unsafe extern "C" fn seed48(seed16v: *mut c_ushort) -> *mut c_ushort { pub unsafe extern "C" fn seed48(seed16v: *mut c_ushort) -> *mut c_ushort {
lcg48::reset_a_and_c(); lcg48::reset_a_and_c();
lcg48::STASHED_XI = lcg48::DEFAULT_XI; // Stash current DEFAULT_XSUBI value in SEED48_XSUBI
lcg48::SEED48_XSUBI = lcg48::DEFAULT_XSUBI;
let new_xi = lcg48::uint48_from_ushort_arr3(seed16v); // Set DEFAULT_XSUBI from the argument provided
lcg48::set_ushort_arr3_from_uint48(lcg48::DEFAULT_XI.as_mut_ptr(), new_xi); let xsubi_value = lcg48::u48_from_ushort_arr3(&*(seed16v as *const [c_ushort; 3]));
lcg48::DEFAULT_XSUBI = lcg48::ushort_arr3_from_u48(xsubi_value);
lcg48::STASHED_XI.as_mut_ptr() // Return the stashed value
lcg48::SEED48_XSUBI.as_mut_ptr()
} }
#[no_mangle] #[no_mangle]
...@@ -873,8 +878,8 @@ pub unsafe extern "C" fn srand48(seedval: c_long) { ...@@ -873,8 +878,8 @@ 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 /* 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 * of the input argument, and the lower 16 bits to 0x330e, as
* specified in POSIX. */ * specified in POSIX. */
let new_xi = (((seedval as u32) as u64) << 16) | 0x330e_u64; let xsubi_value = (u64::from(seedval as u32) << 16) | 0x330e;
lcg48::set_ushort_arr3_from_uint48(lcg48::DEFAULT_XI.as_mut_ptr(), new_xi); lcg48::DEFAULT_XSUBI = lcg48::ushort_arr3_from_u48(xsubi_value);
} }
// #[no_mangle] // #[no_mangle]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment