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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
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
171
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
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/* */
/* This file is part of the HiGHS linear optimization suite */
/* */
/* Available as open-source under the MIT License */
/* */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/**@file util/HighsRandom.h
* @brief Random number generators for HiGHS
*/
#ifndef UTIL_HIGHSRANDOM_H_
#define UTIL_HIGHSRANDOM_H_
#include <cassert>
#include "util/HighsHash.h"
/**
* @brief Class for HiGHS random number generators
*/
class HighsRandom {
private:
uint32_t drawUniform(uint32_t sup, int nbits) {
// draw uniformly in interval [0,sup) where nbits is the maximal number of
// bits the results can have we draw random numbers with nbits many bits
// until we get one that is in the desired range we first use all available
// output functions for the same state before we advance the state again. We
// expect nbits to be at most 32 for this 32 bit version.
assert(sup <= uint32_t{1} << nbits);
while (true) {
advance();
uint32_t lo = static_cast<uint32_t>(state);
uint32_t hi = state >> 32;
uint32_t val;
#define HIGHS_RAND_TRY_OUTPUT(n) \
val = static_cast<uint32_t>(HighsHashHelpers::pair_hash<n>(lo, hi) >> \
(64 - nbits)); \
if (val < sup) return val;
HIGHS_RAND_TRY_OUTPUT(0);
HIGHS_RAND_TRY_OUTPUT(1);
HIGHS_RAND_TRY_OUTPUT(2);
HIGHS_RAND_TRY_OUTPUT(3);
HIGHS_RAND_TRY_OUTPUT(4);
HIGHS_RAND_TRY_OUTPUT(5);
HIGHS_RAND_TRY_OUTPUT(6);
HIGHS_RAND_TRY_OUTPUT(7);
HIGHS_RAND_TRY_OUTPUT(9);
HIGHS_RAND_TRY_OUTPUT(10);
HIGHS_RAND_TRY_OUTPUT(11);
HIGHS_RAND_TRY_OUTPUT(12);
HIGHS_RAND_TRY_OUTPUT(13);
HIGHS_RAND_TRY_OUTPUT(14);
HIGHS_RAND_TRY_OUTPUT(15);
HIGHS_RAND_TRY_OUTPUT(16);
HIGHS_RAND_TRY_OUTPUT(17);
HIGHS_RAND_TRY_OUTPUT(18);
HIGHS_RAND_TRY_OUTPUT(19);
HIGHS_RAND_TRY_OUTPUT(20);
HIGHS_RAND_TRY_OUTPUT(21);
HIGHS_RAND_TRY_OUTPUT(22);
HIGHS_RAND_TRY_OUTPUT(23);
HIGHS_RAND_TRY_OUTPUT(24);
HIGHS_RAND_TRY_OUTPUT(25);
HIGHS_RAND_TRY_OUTPUT(26);
HIGHS_RAND_TRY_OUTPUT(27);
HIGHS_RAND_TRY_OUTPUT(28);
HIGHS_RAND_TRY_OUTPUT(29);
HIGHS_RAND_TRY_OUTPUT(30);
HIGHS_RAND_TRY_OUTPUT(31);
#undef HIGHS_RAND_TRY_OUTPUT
}
}
uint64_t drawUniform(uint64_t sup, int nbits) {
// 64 bit version for drawUniform. If result fits in 32 bits use 32 bit
// version above so that it is only called when we actually need more than
// 32 bits
if (nbits <= 32) return drawUniform(uint32_t(sup), nbits);
assert(sup <= uint64_t{1} << nbits);
while (true) {
advance();
uint32_t lo = static_cast<uint32_t>(state);
uint32_t hi = state >> 32;
uint64_t val;
#define HIGHS_RAND_TRY_OUTPUT(n) \
val = (HighsHashHelpers::pair_hash<2 * n>(lo, hi) >> (64 - nbits)) ^ \
(HighsHashHelpers::pair_hash<2 * n + 1>(lo, hi) >> 32); \
if (val < sup) return val;
HIGHS_RAND_TRY_OUTPUT(0);
HIGHS_RAND_TRY_OUTPUT(1);
HIGHS_RAND_TRY_OUTPUT(2);
HIGHS_RAND_TRY_OUTPUT(3);
HIGHS_RAND_TRY_OUTPUT(4);
HIGHS_RAND_TRY_OUTPUT(5);
HIGHS_RAND_TRY_OUTPUT(6);
HIGHS_RAND_TRY_OUTPUT(7);
HIGHS_RAND_TRY_OUTPUT(9);
HIGHS_RAND_TRY_OUTPUT(10);
HIGHS_RAND_TRY_OUTPUT(11);
HIGHS_RAND_TRY_OUTPUT(12);
HIGHS_RAND_TRY_OUTPUT(13);
HIGHS_RAND_TRY_OUTPUT(14);
HIGHS_RAND_TRY_OUTPUT(15);
#undef HIGHS_RAND_TRY_OUTPUT
}
}
public:
/**
* @brief Initialisations
*/
HighsRandom(HighsUInt seed = 0) { initialise(seed); }
/**
* @brief (Re-)initialise the random number generator
*/
void initialise(HighsUInt seed = 0) {
state = seed;
do {
state = HighsHashHelpers::pair_hash<0>(static_cast<uint32_t>(state),
state >> 32);
state ^= (HighsHashHelpers::pair_hash<1>(state >> 32, seed) >> 32);
} while (state == 0);
}
void advance() {
// advance state with simple xorshift the outputs are produced by applying
// strongly universal hash functions to the state so that the least
// significant bits are as well distributed as the most significant bits.
state ^= state >> 12;
state ^= state << 25;
state ^= state >> 27;
}
/**
* @brief Return a random integer between 0 and 2147483647
*/
HighsInt integer() {
advance();
#ifdef HIGHSINT64
// use 63 bits of the first hash result and use second hash for lower half
// of bits
return (HighsHashHelpers::pair_hash<0>(state, state >> 32) >> 1) ^
(HighsHashHelpers::pair_hash<1>(state, state >> 32) >> 32);
#else
// use 31 bits of the 64 bit result
return HighsHashHelpers::pair_hash<0>(static_cast<uint32_t>(state),
state >> 32) >>
33;
#endif
}
/**
* @brief Return a random integer between [0,sup)
*/
HighsInt integer(HighsInt sup) { // let overload resolution select the 32bit
// or the 64bit version
if (sup <= 1) return 0;
int nbits = HighsHashHelpers::log2i(HighsUInt(sup - 1)) + 1;
return static_cast<HighsInt>(drawUniform(HighsUInt(sup), nbits));
}
/**
* @brief Return a random integer between [min,sup)
*/
HighsInt integer(HighsInt min, HighsInt sup) {
return min + integer(sup - min);
}
/**
* @brief Return a random fraction - real in (0, 1)
*/
double fraction() {
advance();
// 52 bit output is in interval [0,2^52-1]
uint64_t output = (HighsHashHelpers::pair_hash<0>(
static_cast<uint32_t>(state), state >> 32) >>
(64 - 52)) ^
(HighsHashHelpers::pair_hash<1>(
static_cast<uint32_t>(state), state >> 32) >>
(64 - 26));
// compute (1+output) / (2^52+1) which is strictly between 0 and 1
return static_cast<double>(1 + output) * 2.2204460492503125e-16;
}
/**
* @brief Return a random fraction - real in [0, 1]
*/
double closedFraction() {
advance();
// 53 bit result is in interval [0,2^53-1]
uint64_t output = (HighsHashHelpers::pair_hash<0>(
static_cast<uint32_t>(state), state >> 32) >>
(64 - 53)) ^
(HighsHashHelpers::pair_hash<1>(
static_cast<uint32_t>(state), state >> 32) >>
32);
// compute output / (2^53-1) in double precision which is in the closed
// interval [0,1]
return static_cast<double>(output) * 1.1102230246251566e-16;
}
/**
* @brief Return a random real value in the interval [a,b]
*/
double real(double a, double b) { return a + (b - a) * closedFraction(); }
/**
* @brief Return a random bit
*/
bool bit() {
advance();
return HighsHashHelpers::pair_hash<0>(static_cast<uint32_t>(state),
state >> 32) >>
63;
}
/**
* @brief shuffle the given data array
*/
template <typename T>
void shuffle(T* data, HighsInt N) {
for (HighsInt i = N; i > 1; --i) {
HighsInt pos = integer(i);
std::swap(data[pos], data[i - 1]);
}
}
private:
uint64_t state;
};
#endif