Oops,
the first patch was missing. I thought it was only 3 commits, but it was 4.
On Fri, Oct 16, 2020 at 12:02 PM Reini Urban <reini.ur...@gmail.com> wrote:
> Hi
> I'm a GNU maintainer, and want to contribute some new rng's, because I
> stepped over them when updating the DIEHARDER testsuite.
> I'm also maintainer of smhasher.
>
> Do you need a seperate FSF copyright assignment? I'm already on file at
> savannah.
>
> --
> Reini Urban
>
--
Reini Urban
From 227dca23ff166604e4bbf4b540b3c8c935286371 Mon Sep 17 00:00:00 2001
From: Reini Urban <rur...@cpan.org>
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 <config.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <gsl/gsl_rng.h>
+
+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
