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
|
// 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
#pragma once
#include <Progmasoft/Catch3/Internal/RandomFloatingPointHelpers.hpp>
#include <Progmasoft/Catch3/Internal/UniformIntegerDistribution.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
|