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
28namespace 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.
33using Generator =
34 std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>;
35
36using GeneratedWord = typename Generator::result_type;
37static constexpr std::uint64_t range{
38 static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)};
39static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0};
40static constexpr int rangeBits{
41 64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo};
42
43static Lock lock;
44static Generator generator;
45static std::optional<GeneratedWord> nextValue;
46
47// Call only with lock held
48static 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
59template <typename REAL, int PREC>
60inline 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
94extern "C" {
95
96void 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
114void 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
149void 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
175void 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
207void RTNAME(RandomSeedDefaultPut)() {
208 // TODO: should this be time &/or image dependent?
209 {
210 CriticalSection critical{lock};
211 generator.seed(s: 0);
212 }
213}
214
215void 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
247void 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

source code of flang/runtime/random.cpp