846 lines
25 KiB
C
846 lines
25 KiB
C
/**
|
||
* @license Apache-2.0
|
||
*
|
||
* Copyright (c) 2018 The Stdlib Authors.
|
||
*
|
||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||
* you may not use this file except in compliance with the License.
|
||
* You may obtain a copy of the License at
|
||
*
|
||
* http://www.apache.org/licenses/LICENSE-2.0
|
||
*
|
||
* Unless required by applicable law or agreed to in writing, software
|
||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||
* See the License for the specific language governing permissions and
|
||
* limitations under the License.
|
||
*
|
||
*
|
||
* ## Notice
|
||
*
|
||
* The original C code and copyright notice are from the [source implementation]{@link http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c}. The implementation has been modified according to the styles and conventions of this project.
|
||
*
|
||
* ```text
|
||
* Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
|
||
* All rights reserved.
|
||
*
|
||
* Redistribution and use in source and binary forms, with or without
|
||
* modification, are permitted provided that the following conditions
|
||
* are met:
|
||
*
|
||
* 1. Redistributions of source code must retain the above copyright
|
||
* notice, this list of conditions and the following disclaimer.
|
||
*
|
||
* 2. Redistributions in binary form must reproduce the above copyright
|
||
* notice, this list of conditions and the following disclaimer in the
|
||
* documentation and/or other materials provided with the distribution.
|
||
*
|
||
* 3. The names of its contributors may not be used to endorse or promote
|
||
* products derived from this software without specific prior written
|
||
* permission.
|
||
*
|
||
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
||
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
||
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
||
* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
|
||
* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
||
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
||
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
||
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
|
||
* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
||
* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
||
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
||
* ```
|
||
*/
|
||
|
||
// Note: keep project includes in alphabetical order...
|
||
#include <stdlib.h>
|
||
#include <stdint.h>
|
||
#include <string.h>
|
||
#include "stdlib/random/base.h"
|
||
#include "stdlib/random/base/mt19937.h"
|
||
|
||
// Forward declarations:
|
||
static inline int8_t next( struct BasePRNGObject *obj, uint64_t *out );
|
||
static inline int8_t normalized( struct BasePRNGObject *obj, double *out );
|
||
static inline void mt19937_free( struct BasePRNGObject *obj );
|
||
|
||
static inline void create_state( uint32_t *state, const int32_t N, const uint32_t s );
|
||
static inline void init_state( uint32_t *state, const int32_t N, const uint32_t *seed, const int64_t M );
|
||
static inline void twist( uint32_t *state, const int32_t N );
|
||
|
||
// Define the size of the state array (see refs):
|
||
static const int32_t N = 624;
|
||
|
||
// Define a (magic) constant used for indexing into the state array:
|
||
static const int32_t M = 397;
|
||
|
||
// Define an initial state (magic) constant: 19650218 => 00000001001010111101011010101010
|
||
static const uint32_t SEED_ARRAY_INIT_STATE = 19650218;
|
||
|
||
// Define a mask for the most significant `w-r` bits, where `w` is the word size (32 bits) and `r` is the separation point of one word (see refs): 2147483648 => 0x80000000 => 10000000000000000000000000000000
|
||
static const uint32_t UPPER_MASK = 0x80000000;
|
||
|
||
// Define a mask for the least significant `r` bits (see refs): 2147483647 => 0x7fffffff => 01111111111111111111111111111111
|
||
static const uint32_t LOWER_MASK = 0x7fffffff;
|
||
|
||
// Define a multiplier (see Knuth TAOCP Vol2. 3rd Ed. P.106): 1812433253 => 01101100000001111000100101100101
|
||
static const uint32_t KNUTH_MULTIPLIER = 1812433253;
|
||
|
||
// Define a (magic) multiplier: 1664525 => 00000000000110010110011000001101
|
||
static const uint32_t MAGIC_MULTIPLIER_1 = 1664525;
|
||
|
||
// Define a (magic) multiplier: 1566083941 => 01011101010110001000101101100101
|
||
static const uint32_t MAGIC_MULTIPLIER_2 = 1566083941;
|
||
|
||
// Define a tempering coefficient: 2636928640 => 0x9d2c5680 => 10011101001011000101011010000000
|
||
static const uint32_t TEMPERING_COEFFICIENT_1 = 0x9d2c5680;
|
||
|
||
// Define a tempering coefficient: 4022730752 => 0xefc60000 => 11101111110001100000000000000000
|
||
static const uint32_t TEMPERING_COEFFICIENT_2 = 0xefc60000;
|
||
|
||
// Define a constant vector `a` (see refs): 2567483615 => 0x9908b0df => 10011001000010001011000011011111
|
||
static const uint32_t MATRIX_A = 0x9908b0df;
|
||
|
||
// MAG01[x] = x * MATRIX_A; for x = {0,1}
|
||
static const uint32_t MAG01[2] = { 0x0, MATRIX_A };
|
||
|
||
// Define the maximum unsigned 32-bit integer: 4294967295 => 11111111111111111111111111111111
|
||
static const uint32_t MAX_UINT32 = 4294967295;
|
||
|
||
// Define the maximum "safe" double-precision floating-point integer:
|
||
static const double FLOAT64_MAX_SAFE_INTEGER = 9007199254740991.0;
|
||
|
||
// Define a normalization constant when generating double-precision floating-point numbers: 2^53 => 9007199254740992
|
||
static const double FLOAT64_NORMALIZATION_CONSTANT = 1.0 / ( FLOAT64_MAX_SAFE_INTEGER+1.0 );
|
||
|
||
// 2^26: 67108864
|
||
static const double TWO_26 = 67108864.0;
|
||
|
||
// 2^32: 2147483648 => 0x80000000 => 10000000000000000000000000000000
|
||
static const uint32_t TWO_32 = 0x80000000;
|
||
|
||
// 1 (as a 32-bit unsigned integer): 1 => 0x1 => 00000000000000000000000000000001
|
||
static const uint32_t ONE = 0x1;
|
||
|
||
// Define the maximum normalized pseudorandom double-precision floating-point number: ( (((2^32-1)>>>5)*2^26)+( (2^32-1)>>>6) ) / 2^53
|
||
static const double MAX_NORMALIZED = FLOAT64_MAX_SAFE_INTEGER * FLOAT64_NORMALIZATION_CONSTANT;
|
||
|
||
/**
|
||
* MT19937 PRNG.
|
||
*
|
||
* @private
|
||
*/
|
||
static const struct BasePRNG mt19937_prng = {
|
||
"mt19937", // name
|
||
(uint64_t)1, // min
|
||
(uint64_t)MAX_UINT32, // max: (2^{32}-1)
|
||
0.0, // min (normalized)
|
||
MAX_NORMALIZED, // max (normalized): (2^{53}-1)/2^{53}
|
||
sizeof( stdlib_base_random_mt19937_state_t ), // state_size
|
||
&next, // next()
|
||
&normalized, // normalized()
|
||
&mt19937_free // free()
|
||
};
|
||
|
||
/**
|
||
* Returns a pseudorandom integer.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The function returns `-1` if unable to generate a pseudorandom integer and `0` otherwise.
|
||
*
|
||
* ## References
|
||
*
|
||
* - Matsumoto, Makoto, and Takuji Nishimura. 1998. "Mersenne Twister: A 623-dimensionally Equidistributed Uniform Pseudo-random Number Generator." _ACM Transactions on Modeling and Computer Simulation_ 8 (1). New York, NY, USA: ACM: 3–30. doi:[10.1145/272991.272995][@matsumoto:1998a].
|
||
*
|
||
* @private
|
||
* @param obj PRNG object
|
||
* @param out output address
|
||
* @return status code
|
||
*/
|
||
static inline int8_t next( struct BasePRNGObject *obj, uint64_t *out ) {
|
||
stdlib_base_random_mt19937_state_t *so;
|
||
uint32_t *state;
|
||
uint32_t r;
|
||
int32_t i;
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return -1;
|
||
}
|
||
// Retrieve the state object:
|
||
so = (stdlib_base_random_mt19937_state_t *)( obj->state );
|
||
|
||
// Retrieve the current state and current state index:
|
||
state = so->state;
|
||
i = so->i;
|
||
|
||
// Determine if we need to update our internal state array:
|
||
if ( i >= N ) {
|
||
twist( so->state, N );
|
||
i = 0;
|
||
}
|
||
// Get the next state value:
|
||
r = so->state[ i ];
|
||
|
||
// Update the state index:
|
||
so->i = i + 1;
|
||
|
||
// Tempering transform to compensate for the reduced dimensionality of equidistribution:
|
||
r ^= r >> 11;
|
||
r ^= ( r << 7 ) & TEMPERING_COEFFICIENT_1;
|
||
r ^= ( r << 15 ) & TEMPERING_COEFFICIENT_2;
|
||
r ^= r >> 18;
|
||
|
||
// Set the output value:
|
||
*out = (uint64_t)( r );
|
||
|
||
return 0;
|
||
}
|
||
|
||
/**
|
||
* Returns a pseudorandom double-precision floating-point number on the interval `[0,1)`.
|
||
*
|
||
* ## Method
|
||
*
|
||
* - When generating normalized double-precision floating-point numbers, we first generate two pseudorandom integers \\( x \\) and \\( y \\) on the interval \\( [1,2^{32}-1) \\) for a combined \\( 64 \\) random bits.
|
||
*
|
||
* - We would like \\( 53 \\) random bits to generate a 53-bit precision integer and, thus, want to discard \\( 11 \\) of the generated bits.
|
||
*
|
||
* - We do so by discarding \\( 5 \\) bits from \\( x \\) and \\( 6 \\) bits from \\( y \\).
|
||
*
|
||
* - Accordingly, \\( x \\) contains \\( 27 \\) random bits, which are subsequently shifted left \\( 26 \\) bits (multiplied by \\( 2^{26} \\), and \\( y \\) contains \\( 26 \\) random bits to fill in the lower \\( 26 \\) bits. When summed, they combine to comprise \\( 53 \\) random bits of a double-precision floating-point integer.
|
||
*
|
||
* - As an example, suppose, for the sake of argument, the 32-bit PRNG generates the maximum unsigned 32-bit integer \\( 2^{32}-1 \\) twice in a row. Then,
|
||
*
|
||
* ```c
|
||
* uint32_t x = 4294967295 >> 5; // 00000111111111111111111111111111
|
||
* uint32_t y = 4294967295 >> 6; // 00000011111111111111111111111111
|
||
* ```
|
||
*
|
||
* Multiplying \\( x \\) by \\( 2^{26} \\) returns \\( 9007199187632128 \\), which, in binary, is
|
||
*
|
||
* ```binarystring
|
||
* 0 10000110011 11111111111111111111 11111100000000000000000000000000
|
||
* ```
|
||
*
|
||
* Adding \\( y \\) yields \\( 9007199254740991 \\) (the maximum "safe" double-precision floating-point integer value), which, in binary, is
|
||
*
|
||
* ```binarystring
|
||
* 0 10000110011 11111111111111111111 11111111111111111111111111111111
|
||
* ```
|
||
*
|
||
* - Similarly, suppose the 32-bit PRNG generates the following values
|
||
*
|
||
* ```c
|
||
* uint32_t x = 1 >> 5; // 0 => 00000000000000000000000000000000
|
||
* uint32_t y = 64 >> 6; // 1 => 00000000000000000000000000000001
|
||
* ```
|
||
*
|
||
* Multiplying \\( x \\) by \\( 2^{26} \\) returns \\( 0 \\), which, in binary, is
|
||
*
|
||
* ```binarystring
|
||
* 0 00000000000 00000000000000000000 00000000000000000000000000000000
|
||
* ```
|
||
*
|
||
* Adding \\( y \\) yields \\( 1 \\), which, in binary, is
|
||
*
|
||
* ```binarystring
|
||
* 0 01111111111 00000000000000000000 00000000000000000000000000000000
|
||
* ```
|
||
*
|
||
* - As different combinations of \\( x \\) and \\( y \\) are generated, different combinations of double-precision floating-point exponent and significand bits will be toggled, thus generating pseudorandom double-precision floating-point numbers.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The function returns `-1` if unable to generate a pseudorandom number and `0` otherwise.
|
||
*
|
||
* ## References
|
||
*
|
||
* - Harase, Shin. 2017. "Conversion of Mersenne Twister to double-precision floating-point numbers." _ArXiv_ abs/1708.06018 (September). <https://arxiv.org/abs/1708.06018>.
|
||
*
|
||
* @private
|
||
* @param obj PRNG object
|
||
* @param out output address
|
||
* @return status code
|
||
*/
|
||
static inline int8_t normalized( struct BasePRNGObject *obj, double *out ) {
|
||
uint64_t x;
|
||
uint64_t y;
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return -1;
|
||
}
|
||
next( obj, &x);
|
||
x >>= 5;
|
||
|
||
next( obj, &y );
|
||
y >>= 6;
|
||
|
||
// Note: casting `x` and `y` to doubles here is fine, as neither will ever exceed the maximum "safe" double-precision floating-point number:
|
||
*out = ( ((double)x*TWO_26)+(double)y ) * FLOAT64_NORMALIZATION_CONSTANT;
|
||
|
||
return 0;
|
||
}
|
||
|
||
/**
|
||
* Frees a PRNG's allocated memory.
|
||
*
|
||
* @private
|
||
* @param obj PRNG object
|
||
*/
|
||
static inline void mt19937_free( struct BasePRNGObject *obj ) {
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return;
|
||
}
|
||
stdlib_base_random_mt19937_state_t *state = (stdlib_base_random_mt19937_state_t *)( obj->state );
|
||
free( state->seed );
|
||
free( obj->state );
|
||
free( obj );
|
||
}
|
||
|
||
/**
|
||
* Returns an initial PRNG state.
|
||
*
|
||
* @private
|
||
* @param state pointer to state array
|
||
* @param N state array length
|
||
* @param s seed
|
||
*/
|
||
static inline void create_state( uint32_t *state, const int32_t N, const uint32_t s ) {
|
||
int32_t i;
|
||
|
||
// Set the first element of the state array to the provided seed:
|
||
state[ 0 ] = s;
|
||
|
||
// Initialize the remaining state array elements:
|
||
for ( i = 1; i < N; i++ ) {
|
||
state[ i ] = (KNUTH_MULTIPLIER * (state[i-1] ^ (state[i-1]>>30)) + i);
|
||
}
|
||
}
|
||
|
||
/**
|
||
* Initializes a PRNG state array according to a seed array.
|
||
*
|
||
* @private
|
||
* @param state pointer to state array
|
||
* @param N state array length
|
||
* @param seed pointer to seed array
|
||
* @param M seed array length
|
||
*/
|
||
static inline void init_state( uint32_t *state, const int32_t N, const uint32_t *seed, const int64_t M ) {
|
||
int32_t i;
|
||
int64_t j;
|
||
int64_t k;
|
||
|
||
i = 1;
|
||
j = 0;
|
||
k = ( N > M ) ? N : M;
|
||
for ( ; k > 0; k-- ) {
|
||
state[ i ] = (state[i]^((state[i-1]^(state[i-1]>>30))*MAGIC_MULTIPLIER_1)) + seed[j] + j;
|
||
i += 1;
|
||
j += 1;
|
||
if ( i >= N ) {
|
||
state[ 0 ] = state[ N-1 ];
|
||
i = 1;
|
||
}
|
||
if ( j >= M ) {
|
||
j = 0;
|
||
}
|
||
}
|
||
for ( k = N-1; k > 0; k-- ) {
|
||
state[ i ] = (state[i]^((state[i-1]^(state[i-1]>>30))*MAGIC_MULTIPLIER_2)) - i;
|
||
i += 1;
|
||
if ( i >= N ) {
|
||
state[ 0 ] = state[ N-1 ];
|
||
i = 1;
|
||
}
|
||
}
|
||
// Ensure a non-zero initial state array:
|
||
state[ 0 ] = TWO_32; // MSB (most significant bit) is 1
|
||
}
|
||
|
||
/**
|
||
* Updates a PRNG's internal state by generating the next `N` words.
|
||
*
|
||
* @private
|
||
* @param state pointer to a PRNG's internal state array
|
||
* @param N state array length
|
||
*/
|
||
static inline void twist( uint32_t *state, const int32_t N ) {
|
||
uint32_t w;
|
||
int32_t i;
|
||
int32_t j;
|
||
int32_t k;
|
||
|
||
k = N - M;
|
||
for ( i = 0; i < k; i++ ) {
|
||
w = ( state[i]&UPPER_MASK ) | ( state[i+1]&LOWER_MASK );
|
||
state[ i ] = state[ i+M ] ^ ( w>>1 ) ^ MAG01[ w&ONE ];
|
||
}
|
||
j = N - 1;
|
||
for ( ; i < j; i++ ) {
|
||
w = ( state[i]&UPPER_MASK ) | ( state[i+1]&LOWER_MASK );
|
||
state[ i ] = state[ i-k ] ^ ( w>>1 ) ^ MAG01[ w&ONE ];
|
||
}
|
||
w = ( state[j]&UPPER_MASK ) | ( state[0]&LOWER_MASK );
|
||
state[ j ] = state[ M-1 ] ^ ( w>>1 ) ^ MAG01[ w&ONE ];
|
||
}
|
||
|
||
/**
|
||
* Returns a pointer to a dynamically allocated PRNG.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The user is responsible for freeing the allocated memory.
|
||
*
|
||
* @param seed pointer to a seed array
|
||
* @param len seed array length
|
||
* @return pointer to a dynamically allocated PRNG or, if unable to allocate memory, a null pointer
|
||
*
|
||
* @example
|
||
* #include <stdlib.h>
|
||
* #include <stdio.h>
|
||
* #include <stdint.h>
|
||
* #include "stdlib/random/base.h"
|
||
* #include "stdlib/random/base/mt19937.h"
|
||
*
|
||
* // Define a PRNG seed:
|
||
* uint32_t seed[] = { 12345 };
|
||
*
|
||
* // Create a PRNG:
|
||
* struct BasePRNGObject *obj = stdlib_base_random_mt19937_allocate( seed, 1 );
|
||
* if ( obj == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* uint64_t r;
|
||
* int8_t status = obj->prng->next( obj, &r );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Unexpected result.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* // Free allocated memory:
|
||
* stdlib_base_random_mt19937_free( obj );
|
||
*/
|
||
struct BasePRNGObject * stdlib_base_random_mt19937_allocate( const uint32_t *seed, const int64_t len ) {
|
||
stdlib_base_random_mt19937_state_t *state;
|
||
struct BasePRNGObject *obj;
|
||
uint32_t *iseed;
|
||
|
||
obj = malloc( sizeof( struct BasePRNGObject ) );
|
||
if ( obj == NULL ) {
|
||
return NULL;
|
||
}
|
||
state = malloc( sizeof( stdlib_base_random_mt19937_state_t ) );
|
||
if ( state == NULL ) {
|
||
free( obj ); // prevent memory leaks
|
||
return NULL;
|
||
}
|
||
obj->prng = &mt19937_prng;
|
||
obj->state = state;
|
||
|
||
// Create an internal copy of the provided seed array to prevent the inability to reproduce PRNG values based on the PRNG's stated seed due to external state mutation:
|
||
iseed = (uint32_t *)malloc( sizeof( uint32_t )*len );
|
||
if ( iseed == NULL ) {
|
||
free( obj ); // prevent memory leaks
|
||
free( state );
|
||
return NULL;
|
||
}
|
||
memcpy( iseed, seed, sizeof( uint32_t )*len );
|
||
state->seed = iseed;
|
||
state->seed_length = len;
|
||
|
||
// Initialize the PRNG state:
|
||
create_state( state->state, N, SEED_ARRAY_INIT_STATE );
|
||
init_state( state->state, N, iseed, len );
|
||
|
||
// Set the state index which determines when we need to update the PRNG's internal state:
|
||
state->i = N;
|
||
|
||
return obj;
|
||
}
|
||
|
||
/**
|
||
* Frees a PRNG's allocated memory.
|
||
*
|
||
* @param obj PRNG object
|
||
*
|
||
* @example
|
||
* #include <stdlib.h>
|
||
* #include <stdio.h>
|
||
* #include <stdint.h>
|
||
* #include "stdlib/random/base.h"
|
||
* #include "stdlib/random/base/mt19937.h"
|
||
*
|
||
* // Define a PRNG seed:
|
||
* uint32_t seed[] = { 12345 };
|
||
*
|
||
* // Create a PRNG:
|
||
* struct BasePRNGObject *obj = stdlib_base_random_mt19937_allocate( seed, 1 );
|
||
* if ( obj == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* uint64_t r;
|
||
* int8_t status = obj->prng->next( obj, &r );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Unexpected result.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* // Free allocated memory:
|
||
* stdlib_base_random_mt19937_free( obj );
|
||
*/
|
||
void stdlib_base_random_mt19937_free( struct BasePRNGObject *obj ) {
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return;
|
||
}
|
||
obj->prng->free( obj );
|
||
}
|
||
|
||
/**
|
||
* Returns the PRNG seed length.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The function returns `-1` if unable to resolve a PRNG seed length and `0` otherwise.
|
||
*
|
||
* @param obj PRNG object
|
||
* @param len output address
|
||
* @return status code
|
||
*
|
||
* @example
|
||
* #include <stdlib.h>
|
||
* #include <stdio.h>
|
||
* #include <stdint.h>
|
||
* #include "stdlib/random/base.h"
|
||
* #include "stdlib/random/base/mt19937.h"
|
||
*
|
||
* // Define a PRNG seed:
|
||
* uint32_t seed1[] = { 12345 };
|
||
*
|
||
* // Create a PRNG:
|
||
* struct BasePRNGObject *obj = stdlib_base_random_mt19937_allocate( seed1, 1 );
|
||
* if ( obj == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* // Get the seed length:
|
||
* int64_t seed_length;
|
||
* int8_t status = stdlib_base_random_mt19937_seed_length( obj, &seed_length );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Error encountered when attempting to retrieve the PRNG seed length.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // Get the PRNG seed:
|
||
* uint32_t *seed2 = (uint32_t *)malloc( sizeof( uint32_t )*seed_length );
|
||
* if ( seed2 == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
* status = stdlib_base_random_mt19937_seed( obj, seed2 );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Error encountered when attempting to retrieve the PRNG seed.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // Use the seed to, e.g., create another PRNG which will generate the same sequence...
|
||
*
|
||
* // Free allocated memory:
|
||
* stdlib_base_random_mt19937_free( obj );
|
||
* free( seed2 );
|
||
*/
|
||
int8_t stdlib_base_random_mt19937_seed_length( const struct BasePRNGObject *obj, int64_t *len ) {
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return -1;
|
||
}
|
||
// Retrieve the state object:
|
||
const stdlib_base_random_mt19937_state_t *state = (stdlib_base_random_mt19937_state_t *)( obj->state );
|
||
|
||
// Set the seed length:
|
||
*len = (int64_t)( state->seed_length );
|
||
|
||
return 0;
|
||
}
|
||
|
||
/**
|
||
* Returns a PRNG seed.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The function returns `-1` if unable to resolve a PRNG seed and `0` otherwise.
|
||
*
|
||
* @param obj PRNG object
|
||
* @param seed output address
|
||
* @return status code
|
||
*
|
||
* @example
|
||
* #include <stdlib.h>
|
||
* #include <stdio.h>
|
||
* #include <stdint.h>
|
||
* #include "stdlib/random/base.h"
|
||
* #include "stdlib/random/base/mt19937.h"
|
||
*
|
||
* // Define a PRNG seed:
|
||
* uint32_t seed1[] = { 12345 };
|
||
*
|
||
* // Create a PRNG:
|
||
* struct BasePRNGObject *obj = stdlib_base_random_mt19937_allocate( seed1, 1 );
|
||
* if ( obj == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* // Get the seed length:
|
||
* int64_t seed_length;
|
||
* int8_t status = stdlib_base_random_mt19937_seed_length( obj, &seed_length );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Error encountered when attempting to retrieve the PRNG seed length.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // Get the PRNG seed:
|
||
* uint32_t *seed2 = (uint32_t *)malloc( sizeof( uint32_t )*seed_length );
|
||
* if ( seed2 == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
* status = stdlib_base_random_mt19937_seed( obj, seed2 );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Error encountered when attempting to retrieve the PRNG seed.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // Use the seed to, e.g., create another PRNG which will generate the same sequence...
|
||
*
|
||
* // Free allocated memory:
|
||
* stdlib_base_random_mt19937_free( obj );
|
||
* free( seed2 );
|
||
*/
|
||
int8_t stdlib_base_random_mt19937_seed( const struct BasePRNGObject *obj, uint32_t *seed ) {
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return -1;
|
||
}
|
||
// Retrieve the state object:
|
||
const stdlib_base_random_mt19937_state_t *state = (stdlib_base_random_mt19937_state_t *)( obj->state );
|
||
|
||
// Copy the seed array:
|
||
memcpy( seed, state->seed, sizeof( uint32_t )*(state->seed_length) );
|
||
|
||
return 0;
|
||
}
|
||
|
||
/**
|
||
* Returns a **copy** of the current PRNG state.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The user is responsible for freeing the allocated memory.
|
||
*
|
||
* @param obj PRNG object
|
||
* @return pointer to a copy of the PRNG's internal state or, if unable to allocate memory, a null pointer
|
||
*
|
||
* @example
|
||
* #include <stdlib.h>
|
||
* #include <stdio.h>
|
||
* #include <stdint.h>
|
||
* #include "stdlib/random/base.h"
|
||
* #include "stdlib/random/base/mt19937.h"
|
||
*
|
||
* // Define a PRNG seed:
|
||
* uint32_t seed[] = { 12345 };
|
||
*
|
||
* // Create a PRNG:
|
||
* struct BasePRNGObject *obj = stdlib_base_random_mt19937_allocate( seed, 1 );
|
||
* if ( obj == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* stdlib_base_random_mt19937_state_t *state = (stdlib_base_random_mt19937_state_t *)stdlib_base_random_mt19937_state( obj );
|
||
* if ( state == NULL ) {
|
||
* fprintf( stderr, "Unable to retrieve PRNG state.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // Use the captured state to, e.g., sync another PRNG or to reset a PRNG to a particular state in order to "replay" generated values at a later point in time...
|
||
*
|
||
* // Free allocated memory:
|
||
* stdlib_base_random_mt19937_free( obj );
|
||
*
|
||
* free( state->seed );
|
||
* free( state );
|
||
*/
|
||
void * stdlib_base_random_mt19937_state( const struct BasePRNGObject *obj ) {
|
||
stdlib_base_random_mt19937_state_t *state;
|
||
stdlib_base_random_mt19937_state_t *so;
|
||
uint64_t nbytes;
|
||
uint32_t *seed;
|
||
if ( obj == NULL || obj->prng != &mt19937_prng ) {
|
||
return NULL;
|
||
}
|
||
state = (stdlib_base_random_mt19937_state_t *)malloc( obj->prng->state_size );
|
||
if ( state == NULL ) {
|
||
return NULL;
|
||
}
|
||
nbytes = sizeof( uint32_t ) * (state->seed_length);
|
||
seed = (uint32_t *)malloc( nbytes );
|
||
if ( seed == NULL ) {
|
||
free( state ); // prevent memory leaks
|
||
return NULL;
|
||
}
|
||
// Copy the state:
|
||
memcpy( state, obj->state, obj->prng->state_size );
|
||
|
||
// Explicitly copy the seed array to prevent external mutation:
|
||
so = (stdlib_base_random_mt19937_state_t *)( obj->state );
|
||
memcpy( seed, so->seed, nbytes );
|
||
state->seed = seed;
|
||
|
||
return (void *)state;
|
||
}
|
||
|
||
/**
|
||
* Sets the PRNG state.
|
||
*
|
||
* ## Notes
|
||
*
|
||
* - The function returns `-1` if unable to set a PRNG state and `0` otherwise.
|
||
*
|
||
* @param obj PRNG object
|
||
* @param state state
|
||
* @return status code
|
||
*
|
||
* @example
|
||
* #include <stdlib.h>
|
||
* #include <stdio.h>
|
||
* #include <stdint.h>
|
||
* #include "stdlib/random/base.h"
|
||
* #include "stdlib/random/base/mt19937.h"
|
||
*
|
||
* // Define a PRNG seed:
|
||
* uint32_t seed[] = { 12345 };
|
||
*
|
||
* // Create a PRNG:
|
||
* struct BasePRNGObject *obj = stdlib_base_random_mt19937_allocate( seed, 1 );
|
||
* if ( obj == NULL ) {
|
||
* fprintf( stderr, "Error allocating memory.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* uint64_t r;
|
||
* int8_t status = obj->prng->next( obj, &r );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Unexpected result.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* // Retrieve the current PRNG state...
|
||
* stdlib_base_random_mt19937_state_t *state = (stdlib_base_random_mt19937_state_t *)stdlib_base_random_mt19937_state( obj );
|
||
* if ( state == NULL ) {
|
||
* fprintf( stderr, "Error encountered when attempting to retrieve PRNG state.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* // Reset the PRNG to a previous state...
|
||
* status = stdlib_base_random_mt19937_set( obj, (void *)state );
|
||
* if ( status != 0 ) {
|
||
* fprintf( stderr, "Error encountered when attempting to set PRNG state.\n" );
|
||
* exit( 1 );
|
||
* }
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* status = obj->prng->next( obj, &r );
|
||
*
|
||
* // ...
|
||
*
|
||
* // Free allocated memory:
|
||
* stdlib_base_random_mt19937_free( obj );
|
||
*
|
||
* free( state->seed );
|
||
* free( state );
|
||
*/
|
||
int8_t stdlib_base_random_mt19937_set( struct BasePRNGObject *obj, const void *state ) {
|
||
stdlib_base_random_mt19937_state_t *vstate;
|
||
uint64_t nbytes;
|
||
uint32_t *seed;
|
||
if ( obj == NULL || state == NULL || obj->prng != &mt19937_prng ) {
|
||
return -1;
|
||
}
|
||
// Copy the provided seed array:
|
||
vstate = ( stdlib_base_random_mt19937_state_t *)state;
|
||
nbytes = sizeof( uint32_t ) * ( vstate->seed_length );
|
||
seed = (uint32_t *)malloc( nbytes );
|
||
if ( seed == NULL ) {
|
||
return -1;
|
||
}
|
||
memcpy( seed, vstate->seed, nbytes );
|
||
|
||
// Retrieve the current PRNG state:
|
||
vstate = ( stdlib_base_random_mt19937_state_t *)( obj->state );
|
||
|
||
// Free the memory held by the current seed array:
|
||
free( vstate->seed );
|
||
|
||
// Overwrite the current state with the provided state:
|
||
memcpy( vstate, state, obj->prng->state_size );
|
||
|
||
// Update the seed array pointer:
|
||
vstate->seed = seed;
|
||
|
||
return 0;
|
||
}
|