Commit dc9e9853 authored by Christoph Conrads's avatar Christoph Conrads

Add code for RANLUX initialization problems

Add code making plots like in Figure 1 of

Matsumoto et al.: "Common Defects in Initialization of Pseudorandom
Number Generators". 2007.
parent 9910515c
#!/bin/sh
# Copyright 2019 Christoph Conrads
#
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
set -e
c++ -Wextra -Wall -std=c++11 -pedantic \
ranlux-over-seed-lcg.cpp -o ranlux-over-seed-lcg
c++ -Wextra -Wall -std=c++11 -pedantic \
-O3 ranlux-over-seed-counter.cpp -o ranlux-over-seed-counter
c++ -Wextra -Wall -std=c++11 -pedantic \
ranlux-over-seed-lfsr.cpp -o ranlux-over-seed-lfsr
# set the python matplotlib back-end
export MPLBACKEND="cairo"
./ranlux-over-seed-lcg | \
python3 plot-ranlux-over-seed.py \
'LCG initialization' \
ranlux-over-seed-lcg.png
./ranlux-over-seed-counter | \
python3 plot-ranlux-over-seed.py \
'Counter initialization, 10^6 values discarded' \
ranlux-over-seed-counter.png
./ranlux-over-seed-lfsr | \
python3 plot-ranlux-over-seed.py \
'LFSR initialization' \
ranlux-over-seed-lfsr.png
#!/usr/bin/python3
# Copyright 2019 Christoph Conrads
#
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
# This program reads data created by `ranlux-over-seed.cpp` from standard input
# and plots it.
import matplotlib.pyplot as pp
import numpy as np
import sys
def main():
xs = np.loadtxt(sys.stdin, dtype=int)
if xs.shape[1] != 3:
return 1
ns = xs[0,:]
xs = xs[1:,:]
pp.plot(xs[:,0], 'b.', linestyle="None", label='n={:d}'.format(ns[0]))
pp.plot(xs[:,1], 'g+', linestyle="None", label='n={:d}'.format(ns[1]))
pp.plot(xs[:,2], 'rx', linestyle="None", label='n={:d}'.format(ns[2]))
pp.ylim(top=1.025 * 2**24);
pp.ylim(bottom=-0.025 * 2**24);
pp.xlabel('Seed');
pp.ylabel('std::ranlux24_base output');
pp.legend()
title = '' if len(sys.argv) < 2 else sys.argv[1]
if title != '':
pp.title(title)
if len(sys.argv) > 2:
filename = sys.argv[2]
print('Saving plot to "{:s}"'.format(filename))
pp.savefig(filename)
else:
pp.show()
if __name__ == '__main__':
main()
// Copyright 2019 Christoph Conrads
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include <cassert>
#include <cstdio>
#define private public
#include <random>
// This code shows the n-th draw for a ranlux24 pseudo-random number generator
// with seeds from 1 to 100 after skipping a large number of initial values. If
// you plot drawn values over the seeds, patterns will emerge if the PRNG used
// to initialize the ranlux generator is similar in mathematical structure,
// e.g., a linear congruential generator.
int main()
{
constexpr auto n1 = 11u;
constexpr auto n2 = 17u;
constexpr auto n3 = 20u;
static_assert(n1 < n2, "");
static_assert(n2 < n3, "");
std::printf("%u %u %u\n", n1, n2, n3);
for(auto i = 0u; i < 100u; ++i)
{
auto gen = std::ranlux24_base();
std::fill_n(gen._M_x, 24u, 0u);
for(auto j = 0u; j < 12u; ++j)
gen._M_x[j] = i + 1u;
gen.discard(1000u*1000u);
gen.discard(n1-1u);
auto x1 = gen();
gen.discard(n2-n1-1u);
auto x2 = gen();
gen.discard(n3-n2-1u);
auto x3 = gen();
std::printf("%lu %lu %lu\n", x1, x2, x3);
}
}
// Copyright 2019 Christoph Conrads
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include <cassert>
#include <cstdio>
#include <random>
// This code shows the n-th draw for a ranlux24 pseudo-random number generator
// with seeds from 1 to 100. If you plot drawn values over the seeds, patterns
// will emerge if the PRNG used to initialize the ranlux generator is similar in
// mathematical structure, e.g., a linear congruential generator.
int main()
{
constexpr auto n1 = 11u;
constexpr auto n2 = 17u;
constexpr auto n3 = 20u;
static_assert(n1 < n2, "");
static_assert(n2 < n3, "");
std::printf("%u %u %u\n", n1, n2, n3);
for(auto i = 0u; i < 100u; ++i)
{
auto seed = i + 1u;
auto gen = std::ranlux24_base(seed);
gen.discard(n1-1u);
auto x1 = gen();
gen.discard(n2-n1-1u);
auto x2 = gen();
gen.discard(n3-n2-1u);
auto x3 = gen();
std::printf("%lu %lu %lu\n", x1, x2, x3);
}
}
// Copyright 2019 Christoph Conrads
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include <cassert>
#include <cstdio>
#include <cstdint>
// we have to modify the ranlux24 state
#define private public
#include <random>
// This code shows the n-th draw for a ranlux24 pseudo-random number generator
// with seeds from 1 to 100 after skipping a large number of initial values. If
// you plot drawn values over the seeds, *no* patterns should emerge because the
// ranlux generator is initialized with a PRNG with different mathematical
// structure (a linear feedback shift generator).
//
// References:
// * D. Blackman, S. Vigna: "Scrambled Linear Pseudorandom Number Generators."
// 2018. arXiv:1805.01407v1.
// * M. Matsumoto et al.: "Defects in Initialization of Pseudorandom Number
// Generators." In: ACM Transactions on Modeling and Computer Simulation,
// Vol. 17, No. 4, Article 15. 2007. DOI: 10.1145/1276927.1276928
constexpr std::uint32_t rotl(std::uint32_t x, unsigned k)
{
return assert(k < 32), (x << k) | (x >> (32 - k));
}
/**
* This class implements a xoshiro128+ pseudo-random number generator.
*
* D. Blackman, S. Vigna: Scrambled Linear Pseudorandom Number Generators. 2018.
* arXiv:1805.01407v1
*
* See Table 4 of the paper for TestU01 BigCrush results (column "S" shows
* "systemic" failures and column "R" contains "repeated" failures").
*/
struct xoshiro128plus
{
using result_type = std::uint32_t;
explicit xoshiro128plus(std::uint32_t seed) :
xoshiro128plus(std::uint64_t{seed})
{}
explicit xoshiro128plus(std::uint64_t seed) :
// initialize state with some ones here to avoid LFSR zeroland
state_{
0, ~std::uint32_t{0}, std::uint32_t(seed>>32), std::uint32_t(seed)
}
{
// escape from zeroland, cf. Section 9 in Blackman/Vigna (2018)
// this generator is fast so generously discard a few values
discard(128);
}
std::uint32_t operator() ()
{
auto retval = state_[0] + state_[3];
const uint32_t t = state_[1] << 9;
state_[2] ^= state_[0];
state_[3] ^= state_[1];
state_[1] ^= state_[2];
state_[0] ^= state_[3];
state_[2] ^= t;
state_[3] = rotl(state_[3], 11);
return retval;
}
void discard(std::size_t n)
{
for(auto i = std::size_t{0}; i < n; ++i)
(*this)();
}
std::uint32_t state_[4];
};
int main()
{
constexpr auto n1 = 11u;
constexpr auto n2 = 17u;
constexpr auto n3 = 20u;
static_assert(n1 < n2, "");
static_assert(n2 < n3, "");
std::printf("%u %u %u\n", n1, n2, n3);
for(auto i = 0u; i < 100u; ++i)
{
auto seed_gen = xoshiro128plus(i);
auto gen = std::ranlux24_base();
for(auto j = 0u; j < 24u; ++j)
{
gen._M_x[j] = seed_gen() >> 8;
}
gen.discard(n1-1u);
auto x1 = gen();
gen.discard(n2-n1-1u);
auto x2 = gen();
gen.discard(n3-n2-1u);
auto x3 = gen();
std::printf("%lu %lu %lu\n", x1, x2, x3);
}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment