ranlux-over-seed-lfsr.cpp 2.95 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
// 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);
	}
}