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
128
129
130
131
|
// Copyright Catch2 Authors
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE.txt or copy at
// https://www.boost.org/LICENSE_1_0.txt)
// SPDX-License-Identifier: BSL-1.0
#ifndef CATCH_UNIFORM_FLOATING_POINT_DISTRIBUTION_HPP_INCLUDED
#define CATCH_UNIFORM_FLOATING_POINT_DISTRIBUTION_HPP_INCLUDED
#include <catch2/internal/catch_random_floating_point_helpers.hpp>
#include <catch2/internal/catch_uniform_integer_distribution.hpp>
#include <cmath>
#include <type_traits>
namespace Catch {
namespace Detail {
#if defined( __GNUC__ ) || defined( __clang__ )
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
// The issue with overflow only happens with maximal ULP and HUGE
// distance, e.g. when generating numbers in [-inf, inf] for given
// type. So we only check for the largest possible ULP in the
// type, and return something that does not overflow to inf in 1 mult.
constexpr std::uint64_t calculate_max_steps_in_one_go(double gamma) {
if ( gamma == 1.99584030953472e+292 ) { return 9007199254740991; }
return static_cast<std::uint64_t>( -1 );
}
constexpr std::uint32_t calculate_max_steps_in_one_go(float gamma) {
if ( gamma == 2.028241e+31f ) { return 16777215; }
return static_cast<std::uint32_t>( -1 );
}
#if defined( __GNUC__ ) || defined( __clang__ )
# pragma GCC diagnostic pop
#endif
}
/**
* Implementation of uniform distribution on floating point numbers.
*
* Note that we support only `float` and `double` types, because these
* usually mean the same thing across different platform. `long double`
* varies wildly by platform and thus we cannot provide reproducible
* implementation. Also note that we don't implement all parts of
* distribution per standard: this distribution is not serializable, nor
* can the range be arbitrarily reset.
*
* The implementation also uses different approach than the one taken by
* `std::uniform_real_distribution`, where instead of generating a number
* between [0, 1) and then multiplying the range bounds with it, we first
* split the [a, b] range into a set of equidistributed floating point
* numbers, and then use uniform int distribution to pick which one to
* return.
*
* This has the advantage of guaranteeing uniformity (the multiplication
* method loses uniformity due to rounding when multiplying floats), except
* for small non-uniformity at one side of the interval, where we have
* to deal with the fact that not every interval is splittable into
* equidistributed floats.
*
* Based on "Drawing random floating-point numbers from an interval" by
* Frederic Goualard.
*/
template <typename FloatType>
class uniform_floating_point_distribution {
static_assert(std::is_floating_point<FloatType>::value, "...");
static_assert(!std::is_same<FloatType, long double>::value,
"We do not support long double due to inconsistent behaviour between platforms");
using WidthType = Detail::DistanceType<FloatType>;
FloatType m_a, m_b;
FloatType m_ulp_magnitude;
WidthType m_floats_in_range;
uniform_integer_distribution<WidthType> m_int_dist;
// In specific cases, we can overflow into `inf` when computing the
// `steps * g` offset. To avoid this, we don't offset by more than this
// in one multiply + addition.
WidthType m_max_steps_in_one_go;
// We don't want to do the magnitude check every call to `operator()`
bool m_a_has_leq_magnitude;
public:
using result_type = FloatType;
uniform_floating_point_distribution( FloatType a, FloatType b ):
m_a( a ),
m_b( b ),
m_ulp_magnitude( Detail::gamma( m_a, m_b ) ),
m_floats_in_range( Detail::count_equidistant_floats( m_a, m_b, m_ulp_magnitude ) ),
m_int_dist(0, m_floats_in_range),
m_max_steps_in_one_go( Detail::calculate_max_steps_in_one_go(m_ulp_magnitude)),
m_a_has_leq_magnitude(std::fabs(m_a) <= std::fabs(m_b))
{
assert( a <= b );
}
template <typename Generator>
result_type operator()( Generator& g ) {
WidthType steps = m_int_dist( g );
if ( m_a_has_leq_magnitude ) {
if ( steps == m_floats_in_range ) { return m_a; }
auto b = m_b;
while (steps > m_max_steps_in_one_go) {
b -= m_max_steps_in_one_go * m_ulp_magnitude;
steps -= m_max_steps_in_one_go;
}
return b - steps * m_ulp_magnitude;
} else {
if ( steps == m_floats_in_range ) { return m_b; }
auto a = m_a;
while (steps > m_max_steps_in_one_go) {
a += m_max_steps_in_one_go * m_ulp_magnitude;
steps -= m_max_steps_in_one_go;
}
return a + steps * m_ulp_magnitude;
}
}
result_type a() const { return m_a; }
result_type b() const { return m_b; }
};
} // end namespace Catch
#endif // CATCH_UNIFORM_FLOATING_POINT_DISTRIBUTION_HPP_INCLUDED
|
