1 | //===-- runtime/random.cpp ------------------------------------------------===// |
2 | // |
3 | // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
4 | // See https://llvm.org/LICENSE.txt for license information. |
5 | // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
6 | // |
7 | //===----------------------------------------------------------------------===// |
8 | |
9 | // Implements the intrinsic subroutines RANDOM_INIT, RANDOM_NUMBER, and |
10 | // RANDOM_SEED. |
11 | |
12 | #include "flang/Runtime/random.h" |
13 | #include "lock.h" |
14 | #include "terminator.h" |
15 | #include "flang/Common/float128.h" |
16 | #include "flang/Common/leading-zero-bit-count.h" |
17 | #include "flang/Common/uint128.h" |
18 | #include "flang/Runtime/cpp-type.h" |
19 | #include "flang/Runtime/descriptor.h" |
20 | #include <algorithm> |
21 | #include <cmath> |
22 | #include <cstdint> |
23 | #include <limits> |
24 | #include <memory> |
25 | #include <random> |
26 | #include <time.h> |
27 | |
28 | namespace Fortran::runtime { |
29 | |
30 | // Newer "Minimum standard", recommended by Park, Miller, and Stockmeyer in |
31 | // 1993. Same as C++17 std::minstd_rand, but explicitly instantiated for |
32 | // permanence. |
33 | using Generator = |
34 | std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>; |
35 | |
36 | using GeneratedWord = typename Generator::result_type; |
37 | static constexpr std::uint64_t range{ |
38 | static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)}; |
39 | static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0}; |
40 | static constexpr int rangeBits{ |
41 | 64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo}; |
42 | |
43 | static Lock lock; |
44 | static Generator generator; |
45 | static std::optional<GeneratedWord> nextValue; |
46 | |
47 | // Call only with lock held |
48 | static GeneratedWord GetNextValue() { |
49 | GeneratedWord result; |
50 | if (nextValue.has_value()) { |
51 | result = *nextValue; |
52 | nextValue.reset(); |
53 | } else { |
54 | result = generator(); |
55 | } |
56 | return result; |
57 | } |
58 | |
59 | template <typename REAL, int PREC> |
60 | inline void Generate(const Descriptor &harvest) { |
61 | static constexpr std::size_t minBits{ |
62 | std::max<std::size_t>(a: PREC, b: 8 * sizeof(GeneratedWord))}; |
63 | using Int = common::HostUnsignedIntType<minBits>; |
64 | static constexpr std::size_t words{ |
65 | static_cast<std::size_t>(PREC + rangeBits - 1) / rangeBits}; |
66 | std::size_t elements{harvest.Elements()}; |
67 | SubscriptValue at[maxRank]; |
68 | harvest.GetLowerBounds(at); |
69 | { |
70 | CriticalSection critical{lock}; |
71 | for (std::size_t j{0}; j < elements; ++j) { |
72 | while (true) { |
73 | Int fraction{GetNextValue()}; |
74 | if constexpr (words > 1) { |
75 | for (std::size_t k{1}; k < words; ++k) { |
76 | static constexpr auto rangeMask{ |
77 | (GeneratedWord{1} << rangeBits) - 1}; |
78 | GeneratedWord word{(GetNextValue() - generator.min()) & rangeMask}; |
79 | fraction = (fraction << rangeBits) | word; |
80 | } |
81 | } |
82 | fraction >>= words * rangeBits - PREC; |
83 | REAL next{std::ldexp(static_cast<REAL>(fraction), -(PREC + 1))}; |
84 | if (next >= 0.0 && next < 1.0) { |
85 | *harvest.Element<REAL>(at) = next; |
86 | break; |
87 | } |
88 | } |
89 | harvest.IncrementSubscripts(at); |
90 | } |
91 | } |
92 | } |
93 | |
94 | extern "C" { |
95 | |
96 | void RTNAME(RandomInit)(bool repeatable, bool /*image_distinct*/) { |
97 | // TODO: multiple images and image_distinct: add image number |
98 | { |
99 | CriticalSection critical{lock}; |
100 | if (repeatable) { |
101 | generator.seed(s: 0); |
102 | } else { |
103 | #ifdef CLOCK_REALTIME |
104 | timespec ts; |
105 | clock_gettime(CLOCK_REALTIME, tp: &ts); |
106 | generator.seed(s: ts.tv_sec & ts.tv_nsec); |
107 | #else |
108 | generator.seed(time(nullptr)); |
109 | #endif |
110 | } |
111 | } |
112 | } |
113 | |
114 | void RTNAME(RandomNumber)( |
115 | const Descriptor &harvest, const char *source, int line) { |
116 | Terminator terminator{source, line}; |
117 | auto typeCode{harvest.type().GetCategoryAndKind()}; |
118 | RUNTIME_CHECK(terminator, typeCode && typeCode->first == TypeCategory::Real); |
119 | int kind{typeCode->second}; |
120 | switch (kind) { |
121 | // TODO: REAL (2 & 3) |
122 | case 4: |
123 | Generate<CppTypeFor<TypeCategory::Real, 4>, 24>(harvest); |
124 | return; |
125 | case 8: |
126 | Generate<CppTypeFor<TypeCategory::Real, 8>, 53>(harvest); |
127 | return; |
128 | case 10: |
129 | if constexpr (HasCppTypeFor<TypeCategory::Real, 10>) { |
130 | #if LDBL_MANT_DIG == 64 |
131 | Generate<CppTypeFor<TypeCategory::Real, 10>, 64>(harvest); |
132 | return; |
133 | #endif |
134 | } |
135 | break; |
136 | case 16: |
137 | if constexpr (HasCppTypeFor<TypeCategory::Real, 16>) { |
138 | #if LDBL_MANT_DIG == 113 |
139 | Generate<CppTypeFor<TypeCategory::Real, 16>, 113>(harvest); |
140 | return; |
141 | #endif |
142 | } |
143 | break; |
144 | } |
145 | terminator.Crash( |
146 | "not yet implemented: intrinsic: REAL(KIND=%d) in RANDOM_NUMBER" , kind); |
147 | } |
148 | |
149 | void RTNAME(RandomSeedSize)( |
150 | const Descriptor *size, const char *source, int line) { |
151 | if (!size || !size->raw().base_addr) { |
152 | RTNAME(RandomSeedDefaultPut)(); |
153 | return; |
154 | } |
155 | Terminator terminator{source, line}; |
156 | auto typeCode{size->type().GetCategoryAndKind()}; |
157 | RUNTIME_CHECK(terminator, |
158 | size->rank() == 0 && typeCode && |
159 | typeCode->first == TypeCategory::Integer); |
160 | int sizeArg{typeCode->second}; |
161 | switch (sizeArg) { |
162 | case 4: |
163 | *size->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1; |
164 | break; |
165 | case 8: |
166 | *size->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1; |
167 | break; |
168 | default: |
169 | terminator.Crash( |
170 | "not yet implemented: intrinsic: RANDOM_SEED(SIZE=): size %d\n" , |
171 | sizeArg); |
172 | } |
173 | } |
174 | |
175 | void RTNAME(RandomSeedPut)( |
176 | const Descriptor *put, const char *source, int line) { |
177 | if (!put || !put->raw().base_addr) { |
178 | RTNAME(RandomSeedDefaultPut)(); |
179 | return; |
180 | } |
181 | Terminator terminator{source, line}; |
182 | auto typeCode{put->type().GetCategoryAndKind()}; |
183 | RUNTIME_CHECK(terminator, |
184 | put->rank() == 1 && typeCode && |
185 | typeCode->first == TypeCategory::Integer && |
186 | put->GetDimension(0).Extent() >= 1); |
187 | int putArg{typeCode->second}; |
188 | GeneratedWord seed; |
189 | switch (putArg) { |
190 | case 4: |
191 | seed = *put->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>(); |
192 | break; |
193 | case 8: |
194 | seed = *put->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>(); |
195 | break; |
196 | default: |
197 | terminator.Crash( |
198 | "not yet implemented: intrinsic: RANDOM_SEED(PUT=): put %d\n" , putArg); |
199 | } |
200 | { |
201 | CriticalSection critical{lock}; |
202 | generator.seed(s: seed); |
203 | nextValue = seed; |
204 | } |
205 | } |
206 | |
207 | void RTNAME(RandomSeedDefaultPut)() { |
208 | // TODO: should this be time &/or image dependent? |
209 | { |
210 | CriticalSection critical{lock}; |
211 | generator.seed(s: 0); |
212 | } |
213 | } |
214 | |
215 | void RTNAME(RandomSeedGet)( |
216 | const Descriptor *get, const char *source, int line) { |
217 | if (!get || !get->raw().base_addr) { |
218 | RTNAME(RandomSeedDefaultPut)(); |
219 | return; |
220 | } |
221 | Terminator terminator{source, line}; |
222 | auto typeCode{get->type().GetCategoryAndKind()}; |
223 | RUNTIME_CHECK(terminator, |
224 | get->rank() == 1 && typeCode && |
225 | typeCode->first == TypeCategory::Integer && |
226 | get->GetDimension(0).Extent() >= 1); |
227 | int getArg{typeCode->second}; |
228 | GeneratedWord seed; |
229 | { |
230 | CriticalSection critical{lock}; |
231 | seed = GetNextValue(); |
232 | nextValue = seed; |
233 | } |
234 | switch (getArg) { |
235 | case 4: |
236 | *get->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed; |
237 | break; |
238 | case 8: |
239 | *get->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed; |
240 | break; |
241 | default: |
242 | terminator.Crash( |
243 | "not yet implemented: intrinsic: RANDOM_SEED(GET=): get %d\n" , getArg); |
244 | } |
245 | } |
246 | |
247 | void RTNAME(RandomSeed)(const Descriptor *size, const Descriptor *put, |
248 | const Descriptor *get, const char *source, int line) { |
249 | bool sizePresent = size && size->raw().base_addr; |
250 | bool putPresent = put && put->raw().base_addr; |
251 | bool getPresent = get && get->raw().base_addr; |
252 | if (sizePresent + putPresent + getPresent > 1) |
253 | Terminator{source, line}.Crash( |
254 | "RANDOM_SEED must have either 1 or no arguments" ); |
255 | if (sizePresent) |
256 | RTNAME(RandomSeedSize)(size, source, line); |
257 | else if (putPresent) |
258 | RTNAME(RandomSeedPut)(put, source, line); |
259 | else if (getPresent) |
260 | RTNAME(RandomSeedGet)(get, source, line); |
261 | else |
262 | RTNAME(RandomSeedDefaultPut)(); |
263 | } |
264 | |
265 | } // extern "C" |
266 | } // namespace Fortran::runtime |
267 | |