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
|
// Copyright (c) the JPEG XL Project Authors. All rights reserved.
//
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
#ifndef LIB_JXL_BASE_RANDOM_
#define LIB_JXL_BASE_RANDOM_
// Random number generator + distributions.
// We don't use <random> because the implementation (and thus results) differs
// between libstdc++ and libc++.
#include <stdint.h>
#include <string.h>
#include <algorithm>
#include <cmath>
#include "lib/jxl/base/status.h"
namespace jxl {
struct Rng {
explicit Rng(size_t seed)
: s{static_cast<uint64_t>(0x94D049BB133111EBull),
static_cast<uint64_t>(0xBF58476D1CE4E5B9ull) + seed} {}
// Xorshift128+ adapted from xorshift128+-inl.h
uint64_t operator()() {
uint64_t s1 = s[0];
const uint64_t s0 = s[1];
const uint64_t bits = s1 + s0; // b, c
s[0] = s0;
s1 ^= s1 << 23;
s1 ^= s0 ^ (s1 >> 18) ^ (s0 >> 5);
s[1] = s1;
return bits;
}
// Uniformly distributed int64_t in [begin, end), under the assumption that
// `end-begin` is significantly smaller than 1<<64, otherwise there is some
// bias.
int64_t UniformI(int64_t begin, int64_t end) {
JXL_DASSERT(end > begin);
return static_cast<int64_t>((*this)() %
static_cast<uint64_t>(end - begin)) +
begin;
}
// Same as UniformI, but for uint64_t.
uint64_t UniformU(uint64_t begin, uint64_t end) {
JXL_DASSERT(end > begin);
return (*this)() % (end - begin) + begin;
}
// Uniformly distributed float in [begin, end) range. Note: only 23 bits of
// randomness.
float UniformF(float begin, float end) {
float f;
// Bits of a random [1, 2) float.
uint32_t u = ((*this)() >> (64 - 23)) | 0x3F800000;
static_assert(sizeof(f) == sizeof(u),
"Float and U32 must have the same size");
memcpy(&f, &u, sizeof(f));
// Note: (end-begin) * f + (2*begin-end) may fail to return a number >=
// begin.
return (end - begin) * (f - 1.0f) + begin;
}
// Bernoulli trial
bool Bernoulli(float p) { return UniformF(0, 1) < p; }
// State for geometric distributions.
// The stored value is inv_log_1mp
using GeometricDistribution = float;
static GeometricDistribution MakeGeometric(float p) {
return 1.0 / std::log(1 - p);
}
uint32_t Geometric(const GeometricDistribution& dist) {
float f = UniformF(0, 1);
float inv_log_1mp = dist;
float log = std::log(1 - f) * inv_log_1mp;
return static_cast<uint32_t>(log);
}
template <typename T>
void Shuffle(T* t, size_t n) {
for (size_t i = 0; i + 1 < n; i++) {
size_t a = UniformU(i, n);
std::swap(t[a], t[i]);
}
}
private:
uint64_t s[2];
};
} // namespace jxl
#endif // LIB_JXL_BASE_RANDOM_
|
