uniform-real-distribution.hpp 13.9 KB
Newer Older
1
// Copyright 2018-2019 Christoph Conrads
Christoph Conrads's avatar
Christoph Conrads committed
2
3
4
5
6
7
8
9
10
//
// 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_IMPL_UNIFORM_REAL_DISTRIBUTION_HPP
#define RADEMACHER_FPL_IMPL_UNIFORM_REAL_DISTRIBUTION_HPP

#include <cassert>
11
#include <limits>
12
#include <rademacher-fpl/bit.hpp>
13
#include <rademacher-fpl/math.hpp>
14
#include <rademacher-fpl/type-traits.hpp>
15
#include <rademacher-fpl/uniform-int-distribution.hpp>
Christoph Conrads's avatar
Christoph Conrads committed
16
17
#include <rademacher-fpl/uniform-real-distribution.hpp>
#include <random>
18
#include <type_traits>
Christoph Conrads's avatar
Christoph Conrads committed
19
20


21
22
23
24
25
26
27
#if __GNUC__
#define RADEMACHER_FPL_UNLIKELY(expr) __builtin_expect(!!(expr), 0)
#else
#define RADEMACHER_FPL_UNLIKELY(expr) expr
#endif


Christoph Conrads's avatar
Christoph Conrads committed
28
namespace rademacher_fpl {
29
namespace urd_impl {
Christoph Conrads's avatar
Christoph Conrads committed
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
template<typename T>
T fast_nextsmaller(T x)
{
	static_assert(std::is_floating_point<T>::value, "");

	using UInt = traits::as_unsigned_t<T>;
	using Integer = typename std::make_signed<UInt>::type;

	assert(x > 0 or rademacher_fpl::isfinite(x));

	auto i = Integer{0};

	std::memcpy(&i, &x, sizeof(x));

	i = i == 0 ? std::numeric_limits<Integer>::min() :
		i < 0  ? i+1 :
		i > 0  ? i-1 : -1
	;

	std::memcpy(&x, &i, sizeof(x));

	assert(not rademacher_fpl::isnan(x));

	return x;
}


template<typename T>
T nextsmaller(T x)
{
	return rademacher_fpl::nextafter(x, -std::numeric_limits<T>::infinity());
}

inline float nextsmaller(float x) { return fast_nextsmaller(x); }
inline double nextsmaller(double x) { return fast_nextsmaller(x); }



69
70
71
72
73
74
75

template<
	typename Unsigned,
	class Generator
>
Unsigned draw_exponent(Unsigned a, Unsigned b, Generator* p_generator)
{
76
	static_assert(std::is_integral<Unsigned>::value, "");
77
78
	static_assert(not std::is_signed<Unsigned>::value, "");

79
80
81
82
83
84
85
	namespace rfpl = rademacher_fpl;

	using result_t = typename Generator::result_type;
	using uint_t = rfpl::common_t<Unsigned, result_t>;

	static_assert(std::is_integral<result_t>::value, "");
	static_assert(not std::is_signed<result_t>::value, "");
86
87
88
89
90
91
92
93
94
95
96
97

	static_assert(not std::is_signed<uint_t>::value, "");

	// One may use small unsigned integer types to test the behavior when there
	// are more bins than digits. This assertion checks that we are not
	// /accidentally/ using an integer type wider than `Unsigned`.
	static_assert(
		std::numeric_limits<std::size_t>::digits
			< std::numeric_limits<result_t>::digits
		or std::is_same<Unsigned, uint_t>::value,
		""
	);
98

99
100
	assert(a < b);

101
	if(a + 1u == b)
102
103
104
		return a;


105
106
	constexpr auto zero = uint_t{0};
	constexpr auto one = uint_t{1};
107
	constexpr auto num_digits = uint_t(std::numeric_limits<uint_t>::digits);
108
	// the number of bins for *normal* numbers
109
	auto num_bins = b - std::max(Unsigned{1}, a) + 0u;
110
111
112
113
114
115
	auto denormal_p = a == 0;

	if(num_bins <= num_digits)
	{
		auto l = denormal_p ? zero : one;
		auto r = num_bins == num_digits
116
			? std::numeric_limits<uint_t>::max()
117
			: uint_t((one << num_bins) - 1u);
118
		auto virtual_bin = rfpl::impl::draw_uniform_int(p_generator, l, r);
119
		auto exponent = log2p1(virtual_bin) - (denormal_p ? 0u : 1u) + a;
120
121
122
123
124
125
126

		assert(exponent >= a);
		assert(exponent < b);

		return exponent;
	}

127
128
129
130
131
132
133
134
135
136
137
138
	// use only as many digits as returned by the generator (if it has a
	// power-of-two range) instead of drawing from all possible `uint_t` values.
	// this will speed up generators like `std::ranlux24` or generators
	// returning 32-bit values when `uint_t` is a 64-bit type.
	//
	// do not use constexpr here because of, e.g., boost generators
	const auto num_generator_digits = has_pow2_range<Generator>::value
		? rfpl::log2p1(uint_t(Generator::max() - Generator::min()))
		: num_digits
	;
	const auto d = num_generator_digits;

139
140
	while(true)
	{
141
		auto num_iterations = (num_bins-1u) / d;
142

143
		assert(num_bins > num_iterations * d);
144
145
146

		for(auto i = zero; i < num_iterations; ++i)
		{
147
148
149
150
151
152
153
			auto num_remaining_bins = num_bins - i * d;
			auto num_bins_per_group = num_remaining_bins - d;
			auto l = zero;
			auto r = has_pow2_range<Generator>::value
				? Generator::max()
				: std::numeric_limits<uint_t>::max()
			;
154
			auto group = rfpl::impl::draw_uniform_int(p_generator, l, r);
155
156
157
158

			if(group > 0)
			{
				auto bin =
159
					log2p1(group) - (denormal_p ? 0u : 1u) + num_bins_per_group;
160
161
				auto exponent = bin + a;

162
				assert(bin >= num_bins_per_group + (denormal_p ? 1u : 0u));
163
164
165
166
167
168
169
				assert(exponent >= a);
				assert(exponent < b);

				return exponent;
			}
		}

170
		auto num_remaining_bins = num_bins - num_iterations * d;
171
172

		assert(num_remaining_bins <= num_digits);
173
		assert(num_remaining_bins <= d);
174
175
176

		auto l = zero;
		auto r = num_remaining_bins == num_digits
177
			? std::numeric_limits<uint_t>::max()
178
			: uint_t((one << num_remaining_bins) - 1u);
179
		auto virtual_bin = rfpl::impl::draw_uniform_int(p_generator, l, r);
180
181
182

		if(virtual_bin > 0 or denormal_p)
		{
183
			auto bin = log2p1(virtual_bin) - (denormal_p ? 0u : 1u);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
			auto exponent = bin + a;

			assert(exponent >= a);
			assert(exponent < b);

			return exponent;
		}
	}
}


template<
	typename Real,
	class Generator,
	typename Unsigned
>
Unsigned draw_significand_like(Real, Generator* p_generator)
{
	assert(p_generator);

204
	using uint_t = Unsigned;
205
	using rademacher_fpl::impl::draw_uniform_int;
206

207
	constexpr auto one = uint_t{1};
208
	constexpr auto d = traits::num_significand_digits<Real>::value;
209
	constexpr auto l = uint_t{0};
210
	constexpr auto r = uint_t((one << d) - 1u);
211

212
	static_assert(d <= std::numeric_limits<uint_t>::digits, "");
213

214
	return draw_uniform_int(p_generator, l, r);
215
216
217
218
}



219
220
template<
	typename Real,
221
222
	class GeneratorE,
	class GeneratorS,
223
224
	typename Unsigned
>
225
226
Real make_uniform_random_value(
	Real a, Real b,
227
228
	GeneratorE* p_generator_exponent,
	GeneratorS* p_generator_significand
229
)
Christoph Conrads's avatar
Christoph Conrads committed
230
231
{
	static_assert(
232
		std::is_floating_point<Real>::value, ""
Christoph Conrads's avatar
Christoph Conrads committed
233
234
	);

235
	using uint_t = Unsigned;
236
	using rademacher_fpl::impl::draw_uniform_int;
237

238
239
240
	assert(p_generator_exponent);
	assert(p_generator_significand);

241
	constexpr auto num_digits =
242
		static_cast<uint_t>(std::numeric_limits<uint_t>::digits);
243
	constexpr auto num_significand_digits =
244
		uint_t{traits::num_significand_digits<Real>::value};
245
246
247

	static_assert(num_digits > num_significand_digits, "");

248
249
	if(RADEMACHER_FPL_UNLIKELY(
		not p_generator_exponent
250
251
		or not p_generator_significand
		or a >= b
252
253
		or isnan(a)
		or isnan(b)
254
		))
Christoph Conrads's avatar
Christoph Conrads committed
255
	{
256
		return std::numeric_limits<Real>::quiet_NaN();
Christoph Conrads's avatar
Christoph Conrads committed
257
258
	}

259
260
	constexpr auto zero = uint_t{0};
	constexpr auto one = uint_t{1};
261
262
	constexpr auto max_exponent = traits::max_exponent<Real>::value;
	constexpr auto max_significand = traits::max_significand<Real>::value;
Christoph Conrads's avatar
Christoph Conrads committed
263

264
	auto z = nextsmaller(b);
265

Christoph Conrads's avatar
Christoph Conrads committed
266
	// case 1
267
	if(a < 0 && z <= 0)
Christoph Conrads's avatar
Christoph Conrads committed
268
	{
269
270
271
272
		namespace rfpl = rademacher_fpl;

		const auto inf = std::numeric_limits<Real>::infinity();

273
274
		auto x = rfpl::abs(z);
		auto y = rfpl::nextafter(rfpl::abs(a), inf);
275
276
		auto c = make_uniform_random_value(
			x, y, p_generator_exponent, p_generator_significand
277
		);
Christoph Conrads's avatar
Christoph Conrads committed
278

279
		return -c;
Christoph Conrads's avatar
Christoph Conrads committed
280
281
282
	}

	// case 2
283
	if(a == -z)
Christoph Conrads's avatar
Christoph Conrads committed
284
	{
285
286
287
		auto m = make_uniform_random_value(
			Real{0}, b, p_generator_exponent, p_generator_significand
		);
288
		auto sign = draw_uniform_int(p_generator_significand, zero, one);
Christoph Conrads's avatar
Christoph Conrads committed
289
290
291
292
293

		return sign == 0 ? m : -m;
	}

	// case 3
294
	if(a >= 0 and exponent(a) == exponent(z))
Christoph Conrads's avatar
Christoph Conrads committed
295
296
	{
		auto l = significand(a);
297
		auto r = significand(z);
298
		auto significand = draw_uniform_int(p_generator_significand, l, r);
Christoph Conrads's avatar
Christoph Conrads committed
299
300
301
302
303

		return assemble_like(a, exponent(a), significand);
	}

	// case 4
304
	if(a >= 0 and significand(a) == 0 and significand(b) == 0)
Christoph Conrads's avatar
Christoph Conrads committed
305
	{
306
		assert(exponent(b) >= exponent(a) + 1);
Christoph Conrads's avatar
Christoph Conrads committed
307

308
		auto new_significand = draw_significand_like(a, p_generator_significand);
309
310
		auto exponent_a = uint_t{exponent(a)};
		auto exponent_b = uint_t{exponent(b)};
311
312
313
314
		auto new_exponent =
			draw_exponent(exponent_a, exponent_b, p_generator_exponent);

		return assemble_like(a, new_exponent, new_significand);
Christoph Conrads's avatar
Christoph Conrads committed
315
316
	}

317
	if(a >= 0)
318
319
320
	{
		assert(exponent(a) < exponent(b));

321
322
		constexpr auto n_s = one << num_significand_digits;

323
		auto a_denormal_p = exponent(a) == 0;
324
325
		auto e_a = uint_t{exponent(a)};
		auto e_b = uint_t{exponent(b)};
326
		auto p = uint_t(e_b - e_a);
327
328
329
330
331
332
333
		auto q = p - (a_denormal_p ? one : zero);

		auto n_l = n_s - significand(a);
		auto n_c = n_s;
		auto n_r = significand(b);

		if(q + num_significand_digits < num_digits)
334
		{
335
336
337
338
			auto w_l = 1;
			auto w_c = (one << q) - (a_denormal_p ? one : uint_t{2});
			auto w_r = one << q;

339
340
341
			auto t_l = n_l * w_l;
			auto t_c = n_c * w_c + t_l;
			auto t_r = n_r * w_r + t_c;
342

343
			auto l = zero;
344
345
			auto r = uint_t(t_r - 1u);
			const auto v = draw_uniform_int(p_generator_significand, l, r);
346

347
			constexpr auto e_max = uint_t{traits::max_exponent<Real>::value};
348
349

			auto compute_e_c =
350
				[a_denormal_p, e_a, t_l] (uint_t v) -> uint_t {
351
				auto shift = num_significand_digits + (a_denormal_p ? 0 : 1);
352
				auto w = uint_t(((v - t_l) >> shift) + 1);
353
354
355
				return e_a + log2p1(w);
			};
			auto s_c =
356
				             v < t_l ? v + significand(a) :
357
358
				t_l <= v and v < t_c ? uint_t((v - t_l) & (n_s - 1)) :
				t_c <= v and v < t_r ? uint_t((v - t_c) >> q) :
359
				true                 ? assert(false), one : one
360
			;
361
362
363
364
365
			auto e_c =
				             v < t_l ? e_a :
				t_l <= v and v < t_c ? compute_e_c(v) :
				t_c <= v and v < t_r ? e_b :
				true                 ? assert(false), e_max : e_max
366
367
			;

368
369
			auto c = assemble_like(a, e_c, s_c);

370
			assert(c >= a);
371
			assert(c < b);
372
373
374
375

			return c;
		}

376
		if(a_denormal_p)
377
		{
378
			auto l = zero;
379
			auto r = uint_t(n_s + n_r - 1u);
380
381
382

			while(true)
			{
383
				auto x = draw_uniform_int(p_generator_significand, l, r);
384
385
386
387
388
389
390
391
392
393
394
395
				auto s_c = x - (x < n_s ? zero : n_s);
				auto e_c = x < n_s
					? draw_exponent(e_a, e_b, p_generator_exponent)
					: e_b;

				if(e_c == e_a and s_c < significand(a))
					continue;

				auto c = assemble_like(a, e_c, s_c);

				return c;
			}
396
397
		}

398
		auto l = zero;
399
		auto r = uint_t(n_s + n_r - 1u);
400
401
402

		while(true)
		{
403
			auto vbin = draw_uniform_int(p_generator_significand, l, r);
404
405
406
407
408
409
			auto right_bin_p = vbin >= n_s;

			if(right_bin_p)
				return assemble_like(b, e_b, vbin-n_s);

			auto s = zero;
410
			auto t = uint_t(p + one);
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
			auto e_r = draw_exponent(s, t, p_generator_exponent);
			auto s_c = vbin;

			if(e_r == 0)
				continue;

			assert(e_r >= 1);
			assert(e_r <= p);

			auto e_c = e_r + e_a - 1;

			if(e_c == e_a and s_c < significand(a))
				continue;

			auto c = assemble_like(a, e_c, s_c);

			return c;
		}
429
430
	}

431
432
433
	assert(a < 0 and rademacher_fpl::nextafter(b, a) > 0);

	auto e_a = exponent(a);
434
435
	auto e_z = exponent(z);
	auto e_max = std::max(e_a, e_z);
436
437
438
439
440
441
442
	auto n_l = significand(a) + 1u;
	auto n_c = max_significand + 1u;
	auto n_r = significand(b);

	if(e_max + num_significand_digits + 2u < num_digits)
	{
		auto w_1 = e_a == 0 ? one : one << (e_a - 1u);
443
		auto w_0 = e_z == 0 ? one : one << (e_z - 1u);
444
445

		auto t_l = n_l * w_1;
446
		auto t_1 = t_l + (e_a > 0 ? n_c * w_1 : 0u);
447
		auto t_0 = t_1 + (e_z > 0 ? n_c * w_0 : 0u);
448
449
450
451
452
453
454
		auto t_r = t_0 + n_r * w_0;

		assert(t_l > 0);
		assert(t_1 >= t_l);
		assert(t_0 >= t_1);
		assert(t_r > t_0);
		assert(t_l == t_1 or e_a > 0);
455
		assert(t_1 == t_0 or e_z > 0);
456
457

		auto l = zero;
458
459
		auto r = uint_t(t_r - 1u);
		auto vbin = draw_uniform_int(p_generator_significand, l, r);
460

461
		auto signbit_c = vbin < t_1;
462
463
		auto e_c =
			vbin < t_l ? e_a :
464
			vbin < t_1 ? log2p1((vbin-t_l+1) & (w_1 - 1u)) :
465
			vbin < t_0 ? log2p1((vbin - t_1) & (w_0 - 1u)) :
466
			vbin < t_r ? e_z :
467
468
469
470
            true       ? assert(false), max_exponent : 0u
		;
		auto s_c =
			vbin < t_l ? (vbin - 0)   / w_1 :
471
			vbin < t_1 ? (vbin - t_l) / w_1 :
472
473
474
475
476
477
478
479
480
481
482
483
			vbin < t_0 ? (vbin - t_1) / w_0 :
			vbin < t_r ? (vbin - t_0) / w_0 :
			true       ? assert(false), max_significand : 0u
		;
		auto c = assemble_like(a, signbit_c, e_c, s_c);

		assert(c >= a);
		assert(c < b);

		return c;
	}

484
485
486
487
488
489
490
491
492
493
	assert(e_a > 0);
	assert(e_z > 0);

	auto n_a = significand(a) + 1u;
	auto n_1 = e_a == 0 ? 0u : max_significand + 1u;
	auto n_0 = e_z == 0 ? 0u : max_significand + 1u;
	auto n_z = significand(z) + 1u;

	auto restart_p = true;
	auto vbin_s = zero;
494
	auto signbit_c = false;
495
496
497
498
499

	while(restart_p)
	{
		vbin_s = std::numeric_limits<uint_t>::max();

500
		// decide on signbit
501
		auto e_diff = e_a > e_z ? e_a - e_z : e_z - e_a;
502
503
504
505
506
507
508
		auto vbin_signbit = std::numeric_limits<uint_t>::max();
		auto num_iterations = (e_diff+1u + num_digits-1u) / num_digits;

		for(auto i = one; i <= num_iterations; ++i)
		{
			auto l = zero;
			auto r = i == 1 and (e_diff+1u) % num_digits != 0
509
				? uint_t((1u << ((e_diff+1u) % num_digits)) - 1u)
510
511
512
				: std::numeric_limits<uint_t>::max()
			;

513
			vbin_signbit = draw_uniform_int(p_generator_exponent, l, r);
514
515
516
517
518
519
520
521
522
523

			if(i == 1 and vbin_signbit > r/2u + 1u)
			{
				i = zero;
				continue;
			}

			if(vbin_signbit != 0)
				break;
		}
524

525
		assert(vbin_signbit == 0 or vbin_signbit == 1 or e_a != e_z);
526
527

		signbit_c =
528
529
530
			(e_a == e_z and vbin_signbit==1)
			or (e_a > e_z and vbin_signbit != 0)
			or (e_z > e_a and vbin_signbit == 0)
531
532
533
		;

		auto l_s = zero;
534
		auto r_s = uint_t(std::max(n_a + n_1, n_0 + n_z) - 1u);
535

536
		vbin_s = draw_uniform_int(p_generator_significand, l_s, r_s);
537
		restart_p =
538
539
			(signbit_c and vbin_s >= n_a + n_1)
			or (not signbit_c and vbin_s >= n_0 + n_z)
540
541
542
		;
	}

543
544
545
546
	assert(vbin_s < n_a + n_1 or not signbit_c);
	assert(vbin_s < n_z + n_0 or signbit_c);
	assert(vbin_s <  n_a or e_a >  0 or not signbit_c);
	assert(vbin_s <  n_z or e_z >  0 or signbit_c);
547
548

	auto significand_c =
549
550
551
552
		    signbit_c and vbin_s <  n_a ? vbin_s :
		    signbit_c and vbin_s >= n_a ? vbin_s - n_a :
		not signbit_c and vbin_s <  n_z ? vbin_s :
		not signbit_c and vbin_s >= n_z ? vbin_s - n_z :
553
554
555
556
557
558
		true                                  ? max_significand + 1u : zero
	;
	auto draw_exp = [p_generator_exponent] (uint_t l, uint_t r) -> uint_t {
		return draw_exponent(l, r, p_generator_exponent);
	};
	auto exponent_c =
559
560
561
562
		    signbit_c and vbin_s <  n_a ? e_a :
		    signbit_c and vbin_s >= n_a ? draw_exp(zero, e_a) :
		not signbit_c and vbin_s <  n_z ? e_z :
		not signbit_c and vbin_s >= n_z ? draw_exp(zero, e_z) :
563
564
565
566
567
568
569
570
		true                              ? max_exponent + 1u : zero
	;
	auto c = assemble_like(a, signbit_c, exponent_c, significand_c);

	assert(a <= c);
	assert(c < b);

	return c;
Christoph Conrads's avatar
Christoph Conrads committed
571
572
573
574
575
576
}

}
}

#endif