Eigen  5.0.1
 
Loading...
Searching...
No Matches
RandomImpl.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2024 Charles Schlosser <cs.schlosser@gmail.com>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9// SPDX-License-Identifier: MPL-2.0
10
11#ifndef EIGEN_RANDOM_IMPL_H
12#define EIGEN_RANDOM_IMPL_H
13
14// IWYU pragma: private
15#include "./InternalHeaderCheck.h"
16
17namespace Eigen {
18
19namespace internal {
20
21/****************************************************************************
22 * Implementation of random *
23 ****************************************************************************/
24
25// The SplitMix64 finalizer (Steele, Lea and Flood, "Fast splittable pseudorandom number generators", OOPSLA
26// 2014): a bijection on 64-bit words with full avalanche, Stafford's "Mix13" variant of the MurmurHash3 finalizer.
27// Applied to a Weyl sequence (state += 0x9E3779B97F4A7C15) it is the SplitMix64 generator; applied to a seed plus
28// an index it seeds independent per-index streams.
29EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE numext::uint64_t splitmix64_mix(numext::uint64_t z) {
30 z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
31 z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
32 return z ^ (z >> 31);
33}
34
35template <typename Scalar, bool IsComplex, bool IsInteger>
36struct random_default_impl {};
37
38template <typename Scalar>
39struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
40
41template <typename Scalar>
42inline Scalar random(const Scalar& x, const Scalar& y) {
43 return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y);
44}
45
46template <typename Scalar>
47inline Scalar random() {
48 return EIGEN_MATHFUNC_IMPL(random, Scalar)::run();
49}
50
51// TODO: replace or provide alternatives to this, e.g. std::random_device
52struct eigen_random_device {
53 using ReturnType = int;
54 static constexpr int Entropy = floor_log2((unsigned int)(RAND_MAX) + 1);
55 static constexpr ReturnType Highest = RAND_MAX;
56 static EIGEN_DEVICE_FUNC inline ReturnType run() { return std::rand(); }
57};
58
59// Fill a built-in unsigned integer with numRandomBits beginning with the least significant bit
60template <typename Scalar>
61struct random_bits_impl {
62 EIGEN_STATIC_ASSERT(std::is_unsigned<Scalar>::value, SCALAR MUST BE A BUILT - IN UNSIGNED INTEGER)
63 using RandomDevice = eigen_random_device;
64 using RandomReturnType = typename RandomDevice::ReturnType;
65 static constexpr int kTotalBits = sizeof(Scalar) * CHAR_BIT;
66 static constexpr int kEntropy = plain_enum_min(kTotalBits, RandomDevice::Entropy);
67 // return a Scalar filled with numRandomBits beginning from the least significant bit
68 static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) {
69 eigen_assert((numRandomBits >= 0) && (numRandomBits <= kTotalBits));
70 Scalar randomBits = 0;
71 for (int filledBits = 0; filledBits < numRandomBits; filledBits += kEntropy) {
72 Scalar r = static_cast<Scalar>(RandomDevice::run());
73 int remainingBits = numRandomBits - filledBits;
74 if (remainingBits < kEntropy) {
75 // clear the excess bits to avoid UB and rounding bias
76 r >>= kEntropy - remainingBits;
77 }
78 randomBits |= r << filledBits;
79 }
80 return randomBits;
81 }
82};
83
84template <typename BitsType>
85EIGEN_DEVICE_FUNC inline BitsType getRandomBits(int numRandomBits) {
86 return random_bits_impl<BitsType>::run(numRandomBits);
87}
88
89// random implementation for a built-in floating point type
90template <typename Scalar, bool BuiltIn = std::is_floating_point<Scalar>::value>
91struct random_float_impl {
92 using BitsType = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type;
93 static constexpr EIGEN_DEVICE_FUNC inline int mantissaBits() {
94 const int digits = NumTraits<Scalar>::digits();
95 return digits - 1;
96 }
97 static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) {
98 eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissaBits());
99 BitsType randomBits = getRandomBits<BitsType>(numRandomBits);
100 // if fewer than MantissaBits is requested, shift them to the left
101 randomBits <<= (mantissaBits() - numRandomBits);
102 // randomBits is in the half-open interval [2,4)
103 randomBits |= numext::bit_cast<BitsType>(Scalar(2));
104 // result is in the half-open interval [-1,1)
105 Scalar result = numext::bit_cast<Scalar>(randomBits) - Scalar(3);
106 return result;
107 }
108};
109// random implementation for a custom floating point type
110// uses double as the implementation with a mantissa with a size equal to either the target scalar's mantissa or that of
111// double, whichever is smaller
112template <typename Scalar>
113struct random_float_impl<Scalar, false> {
114 static EIGEN_DEVICE_FUNC inline int mantissaBits() {
115 const int digits = NumTraits<Scalar>::digits();
116 constexpr int kDoubleDigits = NumTraits<double>::digits();
117 return numext::mini(digits, kDoubleDigits) - 1;
118 }
119 static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) {
120 eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissaBits());
121 Scalar result = static_cast<Scalar>(random_float_impl<double>::run(numRandomBits));
122 return result;
123 }
124};
125
126#if !EIGEN_COMP_NVCC
127// random implementation for long double
128// this specialization is not compatible with double-double scalars
129template <bool Specialize = (sizeof(long double) == 2 * sizeof(uint64_t)) &&
130 ((std::numeric_limits<long double>::digits != (2 * std::numeric_limits<double>::digits)))>
131struct random_longdouble_impl {
132 static constexpr int Size = sizeof(long double);
133 static constexpr EIGEN_DEVICE_FUNC int mantissaBits() { return NumTraits<long double>::digits() - 1; }
134 static EIGEN_DEVICE_FUNC inline long double run(int numRandomBits) {
135 eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissaBits());
136 EIGEN_USING_STD(memcpy);
137 int numLowBits = numext::mini(numRandomBits, 64);
138 int numHighBits = numext::maxi(numRandomBits - 64, 0);
139 uint64_t randomBits[2];
140 long double result = 2.0L;
141 memcpy(&randomBits, &result, Size);
142#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
143 randomBits[0] |= getRandomBits<uint64_t>(numLowBits);
144 randomBits[1] |= getRandomBits<uint64_t>(numHighBits);
145#elif __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
146 randomBits[0] |= getRandomBits<uint64_t>(numHighBits);
147 randomBits[1] |= getRandomBits<uint64_t>(numLowBits);
148#else
149#error Unexpected or undefined __BYTE_ORDER__
150#endif
151 memcpy(&result, &randomBits, Size);
152 result -= 3.0L;
153 return result;
154 }
155};
156template <>
157struct random_longdouble_impl<false> {
158 static constexpr EIGEN_DEVICE_FUNC int mantissaBits() { return NumTraits<double>::digits() - 1; }
159 static EIGEN_DEVICE_FUNC inline long double run(int numRandomBits) {
160 return static_cast<long double>(random_float_impl<double>::run(numRandomBits));
161 }
162};
163template <>
164struct random_float_impl<long double> : random_longdouble_impl<> {};
165#endif
166
167template <typename Scalar>
168struct random_default_impl<Scalar, false, false> {
169 using Impl = random_float_impl<Scalar>;
170 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y, int numRandomBits) {
171 Scalar half_x = Scalar(0.5) * x;
172 Scalar half_y = Scalar(0.5) * y;
173 Scalar result = (half_x + half_y) + (half_y - half_x) * run(numRandomBits);
174 // result is in the half-open interval [x, y) -- provided that x < y
175 return result;
176 }
177 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) {
178 return run(x, y, Impl::mantissaBits());
179 }
180 static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { return Impl::run(numRandomBits); }
181 static EIGEN_DEVICE_FUNC inline Scalar run() { return run(Impl::mantissaBits()); }
182};
183
184template <typename Scalar, bool IsSigned = NumTraits<Scalar>::IsSigned, bool BuiltIn = std::is_integral<Scalar>::value>
185struct random_int_impl;
186
187// random implementation for a built-in unsigned integer type
188template <typename Scalar>
189struct random_int_impl<Scalar, false, true> {
190 static constexpr int kTotalBits = sizeof(Scalar) * CHAR_BIT;
191 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) {
192 if (y <= x) return x;
193 Scalar range = y - x;
194 // handle edge case where [x,y] spans the entire range of Scalar
195 if (range == NumTraits<Scalar>::highest()) return run();
196 Scalar count = range + 1;
197 // calculate the number of random bits needed to fill range
198 int numRandomBits = log2_ceil(count);
199 Scalar randomBits;
200 do {
201 randomBits = getRandomBits<Scalar>(numRandomBits);
202 // if the random draw is outside [0, range], try again (rejection sampling)
203 // in the worst-case scenario, the probability of rejection is: 1/2 - 1/2^numRandomBits < 50%
204 } while (randomBits >= count);
205 Scalar result = x + randomBits;
206 return result;
207 }
208 static EIGEN_DEVICE_FUNC inline Scalar run() { return getRandomBits<Scalar>(kTotalBits); }
209};
210
211// random implementation for a built-in signed integer type
212template <typename Scalar>
213struct random_int_impl<Scalar, true, true> {
214 static constexpr int kTotalBits = sizeof(Scalar) * CHAR_BIT;
215 // avoid implicit integral promotion to `int`
216 using BitsType = std::conditional_t<(sizeof(Scalar) < sizeof(int)), unsigned int, std::make_unsigned_t<Scalar> >;
217 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) {
218 if (y <= x) return x;
219 // Avoid overflow by representing `range` as an unsigned type
220 BitsType range = static_cast<BitsType>(y) - static_cast<BitsType>(x);
221 BitsType randomBits = random_int_impl<BitsType>::run(0, range);
222 // Avoid overflow in the case where `x` is negative and there is a large range so
223 // `randomBits` would also be negative if cast to `Scalar` first.
224 Scalar result = static_cast<Scalar>(static_cast<BitsType>(x) + randomBits);
225 return result;
226 }
227 static EIGEN_DEVICE_FUNC inline Scalar run() { return static_cast<Scalar>(getRandomBits<BitsType>(kTotalBits)); }
228};
229
230// todo: custom integers
231template <typename Scalar, bool IsSigned>
232struct random_int_impl<Scalar, IsSigned, false> {
233 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar&, const Scalar&) { return run(); }
234 static EIGEN_DEVICE_FUNC inline Scalar run() {
235 eigen_assert(std::false_type::value && "RANDOM FOR CUSTOM INTEGERS NOT YET SUPPORTED");
236 return Scalar(0);
237 }
238};
239
240template <typename Scalar>
241struct random_default_impl<Scalar, false, true> : random_int_impl<Scalar> {};
242
243template <>
244struct random_impl<bool> {
245 static EIGEN_DEVICE_FUNC inline bool run(const bool& x, const bool& y) {
246 if (y <= x) return x;
247 return run();
248 }
249 static EIGEN_DEVICE_FUNC inline bool run() { return getRandomBits<unsigned>(1) != 0; }
250};
251
252template <typename Scalar>
253struct random_default_impl<Scalar, true, false> {
254 using RealScalar = typename NumTraits<Scalar>::Real;
255 using Impl = random_impl<RealScalar>;
256 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y, int numRandomBits) {
257 return Scalar(Impl::run(x.real(), y.real(), numRandomBits), Impl::run(x.imag(), y.imag(), numRandomBits));
258 }
259 static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) {
260 return Scalar(Impl::run(x.real(), y.real()), Impl::run(x.imag(), y.imag()));
261 }
262 static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) {
263 return Scalar(Impl::run(numRandomBits), Impl::run(numRandomBits));
264 }
265 static EIGEN_DEVICE_FUNC inline Scalar run() { return Scalar(Impl::run(), Impl::run()); }
266};
267
268} // namespace internal
269} // namespace Eigen
270
271#endif // EIGEN_RANDOM_IMPL_H