From 227dca23ff166604e4bbf4b540b3c8c935286371 Mon Sep 17 00:00:00 2001 From: Reini Urban Date: Fri, 16 Oct 2020 09:56:15 +0200 Subject: [PATCH 1/4] add gsl_rng_jsf and gsl_rng_jsf64 --- doc/rng.rst | 30 ++++++++++ rng/Makefile.am | 11 ++-- rng/benchmark.c | 2 + rng/gsl_rng.h | 2 + rng/jsf.c | 144 ++++++++++++++++++++++++++++++++++++++++++++++++ rng/test.c | 7 ++- rng/types.c | 2 + 7 files changed, 192 insertions(+), 6 deletions(-) create mode 100644 rng/jsf.c diff --git doc/rng.rst doc/rng.rst index 5f21b3f7..6aa731cc 100644 --- doc/rng.rst +++ doc/rng.rst @@ -1064,6 +1064,31 @@ significant bits. The seed specifies the initial value, :math:`x_1`. +.. index:: JSF random number generator + +.. var:: gsl_rng_jsf + + This is Bob Jenkin's small pseudo-random number generator (2007) + from https://burtleburtle.net/bob/rand/smallprng.html + + The state needs to contain 4 random 32-bit words, with those + forbidden fixpoints: + + { 0x00000000, 0x00000000, 0x00000000, 0x00000000 }, + { 0x77777777, 0x55555555, 0x11111111, 0x44444444 }, + { 0x5591F2E3, 0x69EBA6CD, 0x2A171E3D, 0x3FD48890 }, + { 0x47CB8D56, 0xAE9B35A7, 0x5C78F4A8, 0x522240FF } + { 0x71AAC8F9, 0x66B4F5D3, 0x1E950B8F, 0x481FEA44 }, + { 0xAB23E5C6, 0xD3D74D9A, 0x542E3C7A, 0x7FA91120 } + + The given set function avoids forbidden states for all given seeds. + + The period of this generator is about :math:`2^{126}`, with a shortest + cycle of :math:`2^{94}` and it uses 4 words of state per generator. + The 32bit variant passes all DIEHARDER tests, the 64bit variant + `gsl_rng_jsf64` fails some. + + Performance =========== @@ -1141,6 +1166,11 @@ available online, * DIEHARD source code, G. Marsaglia, http://stat.fsu.edu/pub/diehard/ +The source code for the DIEHARDER random number generator tests is +available online, + +* DIEHARDER source code, Robert G. Brown, https://webhome.phy.duke.edu/~rgb/General/dieharder.php + A comprehensive set of random number generator tests is available from NIST, diff --git rng/Makefile.am rng/Makefile.am index b0c6b8cc..71ad56f4 100644 --- rng/Makefile.am +++ rng/Makefile.am @@ -1,10 +1,11 @@ noinst_LTLIBRARIES = libgslrng.la +noinst_PROGRAMS = benchmark rng_dump pkginclude_HEADERS = gsl_rng.h AM_CPPFLAGS = -I$(top_srcdir) -libgslrng_la_SOURCES = borosh13.c cmrg.c coveyou.c default.c file.c fishman18.c fishman20.c fishman2x.c gfsr4.c knuthran2.c knuthran.c knuthran2002.c lecuyer21.c minstd.c mrg.c mt.c r250.c ran0.c ran1.c ran2.c ran3.c rand48.c rand.c random.c randu.c ranf.c ranlux.c ranlxd.c ranlxs.c ranmar.c rng.c slatec.c taus.c taus113.c transputer.c tt.c types.c uni32.c uni.c vax.c waterman14.c zuf.c inline.c +libgslrng_la_SOURCES = borosh13.c cmrg.c coveyou.c default.c file.c fishman18.c fishman20.c fishman2x.c gfsr4.c knuthran2.c knuthran.c knuthran2002.c lecuyer21.c minstd.c mrg.c mt.c jsf.c r250.c ran0.c ran1.c ran2.c ran3.c rand48.c rand.c random.c randu.c ranf.c ranlux.c ranlxd.c ranlxs.c ranmar.c rng.c slatec.c taus.c taus113.c transputer.c tt.c types.c uni32.c uni.c vax.c waterman14.c zuf.c inline.c CLEANFILES = test.dat @@ -16,7 +17,7 @@ test_LDADD = libgslrng.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la . TESTS = $(check_PROGRAMS) check_PROGRAMS = test -# benchmark_SOURCES = benchmark.c -# benchmark_LDADD = libgslrng.la ../err/libgslerr.la ../utils/libutils.la -# rng_dump_SOURCES = rng-dump.c -# rng_dump_LDADD = libgslrng.la ../err/libgslerr.la ../utils/libutils.la +benchmark_SOURCES = benchmark.c +benchmark_LDADD = libgslrng.la ../err/libgslerr.la ../utils/libutils.la +rng_dump_SOURCES = rng-dump.c +rng_dump_LDADD = libgslrng.la ../err/libgslerr.la ../utils/libutils.la diff --git rng/benchmark.c rng/benchmark.c index b6278494..71cf45c1 100644 --- rng/benchmark.c +++ rng/benchmark.c @@ -91,6 +91,8 @@ main (void) benchmark(gsl_rng_vax); benchmark(gsl_rng_waterman14); benchmark(gsl_rng_zuf); + benchmark(gsl_rng_jsf); + benchmark(gsl_rng_jsf64); return 0; } diff --git rng/gsl_rng.h rng/gsl_rng.h index 4ec55911..1e3dc61f 100644 --- rng/gsl_rng.h +++ rng/gsl_rng.h @@ -121,6 +121,8 @@ GSL_VAR const gsl_rng_type *gsl_rng_uni32; GSL_VAR const gsl_rng_type *gsl_rng_vax; GSL_VAR const gsl_rng_type *gsl_rng_waterman14; GSL_VAR const gsl_rng_type *gsl_rng_zuf; +GSL_VAR const gsl_rng_type *gsl_rng_jsf; +GSL_VAR const gsl_rng_type *gsl_rng_jsf64; const gsl_rng_type ** gsl_rng_types_setup(void); diff --git rng/jsf.c rng/jsf.c new file mode 100644 index 00000000..2ff98119 --- /dev/null +++ rng/jsf.c @@ -0,0 +1,144 @@ +/* + * rng/jsf.c + * + * Copyright (C) 2007 Bob Jenkins, placed into the Public domain. + * Copyright (C) 2020 Reini Urban, GSL packaging + * + * Jenkin's Small Fast pseudo-random number generator 2007 + * From https://burtleburtle.net/bob/rand/smallprng.html + * + * The 32bit variant passes all dieharder tests, the 64bit variant fails some. + */ + +#include +#include +#include +#include + +static unsigned long int jsf_get (void *vstate); +static double jsf_get_double (void *vstate); +static void jsf_set (void *vstate, unsigned long int s); +static unsigned long int jsf64_get (void *vstate); +static double jsf64_get_double (void *vstate); +static void jsf64_set (void *vstate, unsigned long int s); + +typedef struct { + uint32_t a; + uint32_t b; + uint32_t c; + uint32_t d; +} jsf_state_t; + +static inline uint32_t rotl32(const uint32_t x, int k) { + return (x << k) | (x >> (32 - k)); +} + +static unsigned long int jsf_get (void *vstate) +{ + + jsf_state_t *x = vstate; + uint32_t e; + e = x->a - rotl32(x->b, 27); + x->a = x->b ^ rotl32(x->c, 17); + x->b = x->c + x->d; + x->c = x->d + e; + x->d = e + x->a; + return (unsigned long int)x->d; +} + +static double jsf_get_double (void *vstate) +{ + return (double) jsf_get (vstate) / (double) UINT32_MAX; +} + +/* + (April 2009) Elias Yarrkov used a black-box solver to find some fixed + points of this generator (which produce cycles of length 1): + + { 0x00000000, 0x00000000, 0x00000000, 0x00000000 }, + { 0x77777777, 0x55555555, 0x11111111, 0x44444444 }, + { 0x5591F2E3, 0x69EBA6CD, 0x2A171E3D, 0x3FD48890 }, + { 0x47CB8D56, 0xAE9B35A7, 0x5C78F4A8, 0x522240FF } + + (January 2016) David Blackman found two more, and claims that these + six solutions are all that there are. I (Bob Jekins) have not verified. + + { 0x71AAC8F9, 0x66B4F5D3, 0x1E950B8F, 0x481FEA44 }, + { 0xAB23E5C6, 0xD3D74D9A, 0x542E3C7A, 0x7FA91120 } + */ +static void +jsf_set (void *vstate, unsigned long int s) +{ + jsf_state_t *x = (jsf_state_t *) vstate; + int i; + x->a = 0xf1ea5eed, x->b = x->c = x->d = s; + for (i=0; i<20; ++i) { + (void)jsf_get(x); + } + return; +} + +static const gsl_rng_type jsf_type = +{"jsf", /* name */ + UINT32_MAX, /* RAND_MAX */ + 0, /* RAND_MIN */ + sizeof (jsf_state_t), + &jsf_set, + &jsf_get, + &jsf_get_double}; + +/* ============ 64bit variant ============= */ + +typedef struct { + uint64_t a; + uint64_t b; + uint64_t c; + uint64_t d; +} jsf64_state_t; + +static inline uint64_t rotl64(const uint64_t x, int k) { + return (x << k) | (x >> (64 - k)); +} + +static unsigned long int jsf64_get (void *vstate) +{ + + jsf64_state_t *x = vstate; + uint64_t e; + e = x->a - rotl64(x->b, 7); + x->a = x->b ^ rotl64(x->c, 13); + x->b = x->c + rotl64(x->d, 37); + x->c = x->d + e; + x->d = e + x->a; + return (unsigned long int)x->d; +} + +static double jsf64_get_double (void *vstate) +{ + return (double) ((jsf64_get (vstate) >> 11) * 0x1.0p-53); +} + +static void +jsf64_set (void *vstate, unsigned long int s) +{ + /* Initialize automaton using specified seed. */ + jsf64_state_t *x = (jsf64_state_t *) vstate; + int i; + x->a = 0xf1ea5eed, x->b = x->c = x->d = s; + for (i=0; i<20; ++i) { + (void)jsf64_get(x); + } + return; +} + +static const gsl_rng_type jsf64_type = +{"jsf64", /* name */ + UINT64_MAX, /* RAND_MAX */ + 0, /* RAND_MIN */ + sizeof (jsf64_state_t), + &jsf64_set, + &jsf64_get, + &jsf64_get_double}; + +const gsl_rng_type *gsl_rng_jsf = &jsf_type; +const gsl_rng_type *gsl_rng_jsf64 = &jsf64_type; diff --git rng/test.c rng/test.c index f82ad1ae..e8a29ecd 100644 --- rng/test.c +++ rng/test.c @@ -181,6 +181,9 @@ main (void) rng_test (gsl_rng_ranf, 0, 10000, 2152890433UL); rng_test (gsl_rng_ranf, 2, 10000, 339327233); + rng_test (gsl_rng_jsf, 1, 10000, 210975159); + rng_test (gsl_rng_jsf64, 1, 10000, 10988883387291576526UL); + /* Test constant relationship between int and double functions */ for (r = rngs ; *r != 0; r++) @@ -221,6 +224,8 @@ main (void) rng_seed_test (gsl_rng_taus2); rng_seed_test (gsl_rng_gfsr4); + rng_seed_test (gsl_rng_jsf); + rng_seed_test (gsl_rng_jsf64); #endif exit (gsl_test_summary ()); @@ -247,7 +252,7 @@ rng_test (const gsl_rng_type * T, unsigned long int seed, unsigned int n, } status = (k != result); - gsl_test (status, "%s, %u steps (%u observed vs %u expected)", + gsl_test (status, "%s, %u steps (%lu observed vs %lu expected)", gsl_rng_name (r), n, k, result); gsl_rng_free (r); diff --git rng/types.c rng/types.c index 071edcf6..d550cf76 100644 --- rng/types.c +++ rng/types.c @@ -94,6 +94,8 @@ gsl_rng_types_setup (void) ADD(gsl_rng_vax); ADD(gsl_rng_waterman14); ADD(gsl_rng_zuf); + ADD(gsl_rng_jsf); + ADD(gsl_rng_jsf64); ADD(0); return gsl_rng_generator_types; -- 2.26.2