random-number-engine.hpp 9.33 KB
Newer Older
Christoph Conrads's avatar
Christoph Conrads committed
1 2 3 4 5 6 7 8 9
// 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/.

#ifndef RADEMACHER_FPL_RANDOM_NUMBER_ENGINE_HPP
#define RADEMACHER_FPL_RANDOM_NUMBER_ENGINE_HPP

10
#include <cassert>
Christoph Conrads's avatar
Christoph Conrads committed
11 12 13 14
#include <cstddef>
#include <cstdint>
#include <limits>
#include <random>
15
#include <cstring>
Christoph Conrads's avatar
Christoph Conrads committed
16 17
#include <type_traits>

Christoph Conrads's avatar
Christoph Conrads committed
18 19 20 21 22 23
#if defined(__SIZEOF_INT128__) && __SIZEOF_INT128__ == 16
#define RANLUX_TOOLS_HAS_INT128 1
#else
#define RANLUX_TOOLS_HAS_INT128 0
#endif

Christoph Conrads's avatar
Christoph Conrads committed
24

Christoph Conrads's avatar
Christoph Conrads committed
25
namespace ranlux_tools {
Christoph Conrads's avatar
Christoph Conrads committed
26 27 28 29 30 31
namespace impl_random_number_engine
{
	constexpr std::uint32_t rotl(std::uint32_t x, unsigned k)
	{
		return (x << k) | (x >> (32 - k));
	}
Christoph Conrads's avatar
Christoph Conrads committed
32

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

	// minimize the amount of code needed to switch between AWC/SWB
	// implementations by providing a dummy integer to be used whenever int128
	// is unavailable.
	// this type allows us to use a simple if statement to switch the
	// implementation and the optimizer should take care of the rest.
	struct dummy_integer_t
	{
		dummy_integer_t(std::uintmax_t) {}

		operator std::uintmax_t() const { return 0; }
	};

	dummy_integer_t operator+ (dummy_integer_t, dummy_integer_t) { return 0; }
	dummy_integer_t operator+ (dummy_integer_t, std::uintmax_t) { return 0; }
	dummy_integer_t operator- (dummy_integer_t, dummy_integer_t) { return 0; }
	dummy_integer_t operator- (dummy_integer_t, std::uintmax_t) { return 0; }
	dummy_integer_t operator>> (dummy_integer_t, std::size_t) { return 0;}


Christoph Conrads's avatar
Christoph Conrads committed
53 54 55 56 57 58 59
#if RANLUX_TOOLS_HAS_INT128
#if __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wpedantic"
#endif
	template<std::size_t w>
	using big_integer_t = typename std::conditional<
60
		w == 8u,  std::uint16_t, typename std::conditional<
Christoph Conrads's avatar
Christoph Conrads committed
61 62
		w == 16u, std::uint32_t, typename std::conditional<
		w == 32u, std::uint64_t, typename std::conditional<
63
		w == 64u, unsigned __int128, dummy_integer_t
64
	>::type>::type>::type>::type;
Christoph Conrads's avatar
Christoph Conrads committed
65 66 67 68 69 70
#if __GNUC__
#pragma GCC diagnostic pop
#endif
#else
	template<std::size_t w>
	using big_integer_t = typename std::conditional<
71
		w == 8u,  std::uint16_t, typename std::conditional<
Christoph Conrads's avatar
Christoph Conrads committed
72
		w == 16u, std::uint32_t, typename std::conditional<
73
		w == 32u, std::uint64_t, dummy_integer_t
74
	>::type>::type>::type;
Christoph Conrads's avatar
Christoph Conrads committed
75 76
#endif

Christoph Conrads's avatar
Christoph Conrads committed
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 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
}

/**
 * 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
 * "systematic" failures and column "R" contains "repeated" failures"). The
 * public domain code for xoshiro128plus by Blackman and Vigna can be downloaded
 * from the xoshiro/xoroshiro website:
 *
 *   http://xoshiro.di.unimi.it
 */
struct xoshiro128plus
{
	using result_type = std::uint32_t;

	static constexpr result_type max() {
		return std::numeric_limits<result_type>::max();
	}
	static constexpr result_type min() {
		return std::numeric_limits<result_type>::min();
	}


	explicit xoshiro128plus(std::uint64_t seed=0u) :
		// 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() ()
	{
		using impl_random_number_engine::rotl;

		auto retval = state_[0] + state_[3];
		auto 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(unsigned long long n)
	{
		for(auto i = 0ull; i < n; ++i)
			(*this)();
	}


	std::uint32_t state_[4];
};


145

Christoph Conrads's avatar
Christoph Conrads committed
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
/**
 * This class implements an add-with-carry (AWC) pseudo-random number generator.
 *
 * This implementation supports only the word size 32 with 32-bit unsigned
 * integer types. This generator is initialized with the aid of a xoshiro128+
 * PRNG to avoid correlated outputs for AWC PRNGs with similar seeds.
 *
 * References:
 * * G. Marsaglia, A. Zaman: "A New Class of Random Number Generators." In: The
 *   Annals of Applied Probability, Vol. 1, No. 3, pp. 462--469. 1991.
 * * 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
 */
template<typename T, std::size_t w, std::size_t s, std::size_t r>
struct add_with_carry_engine
{
	static_assert(std::is_integral<T>::value, "");
	static_assert(not std::is_signed<T>::value, "");
	static_assert(w <= std::numeric_limits<T>::digits, "");
	static_assert(s > 0u, "");
	static_assert(r > s, "");

	// prototype implementation
	static_assert(w == std::numeric_limits<T>::digits, "");
171
	static_assert(w == 8u or w == 16u or w == 32u or w == 64u, "");
Christoph Conrads's avatar
Christoph Conrads committed
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201

	using result_type = T;

	static constexpr auto long_lag = r;
	static constexpr auto short_lag = s;
	static constexpr auto word_size = w;
	// `std::subtract_with_carry_engine::default_seed` is 19780503. I have no
	// idea where this value is coming from so I use its larger prime factor
	// instead.
	static constexpr auto default_seed = std::uint32_t{387853};


	static constexpr T max() { return std::numeric_limits<T>::max(); }
	static constexpr T min() { return std::numeric_limits<T>::min(); }


	explicit add_with_carry_engine(std::uint64_t seed=default_seed)
	{
		auto gen = xoshiro128plus(seed);

		for(auto& x : xs_)
			x = gen();

		// ensure entering a periodic sequence
		discard(r);
	}


	T operator() ()
	{
Christoph Conrads's avatar
Christoph Conrads committed
202
		using BigInt = impl_random_number_engine::big_integer_t<w>;
203
		using DummyInt = impl_random_number_engine::dummy_integer_t;
204 205 206 207 208 209

		auto i = index_;
		auto j = index_ >= s ? index_ - s : index_ + r - s;

		assert(carry_ == 0 or carry_ == 1);

210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
		if(std::is_same<BigInt, DummyInt>::value)
		{
			auto x = xs_[i];
			auto y = T(x + xs_[j]);
			auto z = T(y + carry_);

			xs_[i] = z;
			carry_ = y < x or z < y ? 1u : 0u;
		}
		else
		{
			auto x = BigInt{xs_[i]} + xs_[j] + carry_;

			xs_[i] = x;
			carry_ = x >> w;
		}
226 227 228

		index_ = index_ + 1u == r ? 0u : index_ + 1u;

229
		return xs_[i];
230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
	}



	void discard(unsigned long long n)
	{
		for(auto i = 0ull; i < n; ++i)
			(*this)();
	}


	std::size_t index_ = 0;
	T carry_ = 0;
	T xs_[r] = { 0 };
};


247 248

template<typename T, std::size_t w, std::size_t p, std::size_t q>
249 250
struct subtract_with_borrow_engine
{
251 252 253 254 255 256 257 258

	static constexpr auto word_size = w;
	static constexpr auto long_lag = p > q ? p : q;
	static constexpr auto short_lag = p > q ? q : p;

	static constexpr auto r = long_lag;
	static constexpr auto s = short_lag;

259 260 261
	static_assert(std::is_integral<T>::value, "");
	static_assert(not std::is_signed<T>::value, "");
	static_assert(w <= std::numeric_limits<T>::digits, "");
262 263
	static_assert(short_lag > 0u, "");
	static_assert(long_lag > short_lag, "");
264 265 266

	// prototype implementation
	static_assert(w == std::numeric_limits<T>::digits, "");
267
	static_assert(w == 8u or w == 16u or w == 32u or w == 64u, "");
268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294

	using result_type = T;

	// `std::subtract_with_carry_engine::default_seed` is 19780503. I have no
	// idea where this value is coming from so I use its larger prime factor
	// instead.
	static constexpr auto default_seed = std::uint32_t{387853};


	static constexpr T max() { return std::numeric_limits<T>::max(); }
	static constexpr T min() { return std::numeric_limits<T>::min(); }


	explicit subtract_with_borrow_engine(std::uint64_t seed=default_seed)
	{
		auto gen = xoshiro128plus(seed);

		for(auto& x : xs_)
			x = gen();

		// ensure entering a periodic sequence
		discard(r);
	}


	T operator() ()
	{
Christoph Conrads's avatar
Christoph Conrads committed
295
		using BigInt = impl_random_number_engine::big_integer_t<w>;
296
		using DummyInt = impl_random_number_engine::dummy_integer_t;
Christoph Conrads's avatar
Christoph Conrads committed
297 298 299 300 301 302

		auto i = index_;
		auto j = index_ >= s ? index_ - s : index_ + r - s;

		assert(carry_ == 0 or carry_ == 1);

303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
		if(std::is_same<BigInt, DummyInt>::value)
		{
			auto x = p > q ? xs_[i] : xs_[j];
			auto y = p > q ? T(x - xs_[j]) : T(x - xs_[i]);
			auto z = T(y - carry_);

			xs_[i] = z;
			carry_ = y > x or z > y ? 1u : 0u;
		}
		else
		{
			auto x = p > q
				? BigInt{xs_[i]} - xs_[j] - carry_
				: BigInt{xs_[j]} - xs_[i] - carry_
			;
			auto c = -T(x >> w);

			xs_[i] = x;
			carry_ = c;
		}
Christoph Conrads's avatar
Christoph Conrads committed
323 324 325

		index_ = index_ + 1u == r ? 0u : index_ + 1u;

326
		return xs_[i];
Christoph Conrads's avatar
Christoph Conrads committed
327 328 329 330 331 332 333 334 335 336 337
	}


	void discard(unsigned long long n)
	{
		for(auto i = 0ull; i < n; ++i)
			(*this)();
	}


	std::size_t index_ = 0;
338
	T carry_ = 0;
339
	T xs_[long_lag] = { 0 };
Christoph Conrads's avatar
Christoph Conrads committed
340 341 342
};


343

Christoph Conrads's avatar
Christoph Conrads committed
344 345 346
using ranlux8_base = subtract_with_borrow_engine<std::uint8_t, 8u, 5u, 8u>;
using ranlux8 =  std::discard_block_engine<ranlux8_base, 67u, 8u>;
using fast_ranlux8  = std::discard_block_engine<ranlux8_base, 17u, 8u>;
347

Christoph Conrads's avatar
Christoph Conrads committed
348 349 350
using ranlux16_base = add_with_carry_engine<std::uint16_t, 16u, 2u, 9u>;
using ranlux16 = std::discard_block_engine<ranlux16_base, 97u, 9u>;
using fast_ranlux16 = std::discard_block_engine<ranlux16_base, 23u, 9u>;
351

Christoph Conrads's avatar
Christoph Conrads committed
352 353 354
using ranlux32_base = add_with_carry_engine<std::uint32_t, 32u, 3u, 16u>;
using ranlux32 = std::discard_block_engine<ranlux32_base, 277u, 16u>;
using fast_ranlux32 = std::discard_block_engine<ranlux32_base, 71u, 16u>;
355

Christoph Conrads's avatar
Christoph Conrads committed
356 357 358
using ranlux64_base = subtract_with_borrow_engine<std::uint64_t, 64u, 62u, 3u>;
using ranlux64 = std::discard_block_engine<ranlux64_base, 1303u, 62u>;
using fast_ranlux64 = std::discard_block_engine<ranlux64_base, 331u, 62u>;
Christoph Conrads's avatar
Christoph Conrads committed
359 360 361 362

}

#endif