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 "random.h"
13 #include "cpp-type.h"
14 #include "descriptor.h"
15 #include "lock.h"
16 #include "flang/Common/leading-zero-bit-count.h"
17 #include "flang/Common/uint128.h"
18 #include <algorithm>
19 #include <cmath>
20 #include <cstdint>
21 #include <ctime>
22 #include <limits>
23 #include <memory>
24 #include <random>
25
26 namespace Fortran::runtime {
27
28 // Newer "Minimum standard", recommended by Park, Miller, and Stockmeyer in
29 // 1993. Same as C++17 std::minstd_rand, but explicitly instantiated for
30 // permanence.
31 using Generator =
32 std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>;
33
34 using GeneratedWord = typename Generator::result_type;
35 static constexpr std::uint64_t range{
36 static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)};
37 static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0};
38 static constexpr int rangeBits{
39 64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo};
40
41 static Lock lock;
42 static Generator generator;
43
44 template <typename REAL, int PREC>
Generate(const Descriptor & harvest)45 inline void Generate(const Descriptor &harvest) {
46 static constexpr std::size_t minBits{
47 std::max<std::size_t>(PREC, 8 * sizeof(GeneratedWord))};
48 using Int = common::HostUnsignedIntType<minBits>;
49 static constexpr std::size_t words{
50 static_cast<std::size_t>(PREC + rangeBits - 1) / rangeBits};
51 std::size_t elements{harvest.Elements()};
52 SubscriptValue at[maxRank];
53 harvest.GetLowerBounds(at);
54 {
55 CriticalSection critical{lock};
56 for (std::size_t j{0}; j < elements; ++j) {
57 Int fraction{generator()};
58 if constexpr (words > 1) {
59 for (std::size_t k{1}; k < words; ++k) {
60 static constexpr auto rangeMask{(GeneratedWord{1} << rangeBits) - 1};
61 GeneratedWord word{(generator() - generator.min()) & rangeMask};
62 fraction = (fraction << rangeBits) | word;
63 }
64 }
65 fraction >>= words * rangeBits - PREC;
66 *harvest.Element<REAL>(at) =
67 std::ldexp(static_cast<REAL>(fraction), -(PREC + 1));
68 harvest.IncrementSubscripts(at);
69 }
70 }
71 }
72
73 extern "C" {
74
RTNAME(RandomInit)75 void RTNAME(RandomInit)(bool repeatable, bool /*image_distinct*/) {
76 // TODO: multiple images and image_distinct: add image number
77 {
78 CriticalSection critical{lock};
79 if (repeatable) {
80 generator.seed(0);
81 } else {
82 generator.seed(std::time(nullptr));
83 }
84 }
85 }
86
RTNAME(RandomNumber)87 void RTNAME(RandomNumber)(
88 const Descriptor &harvest, const char *source, int line) {
89 Terminator terminator{source, line};
90 auto typeCode{harvest.type().GetCategoryAndKind()};
91 RUNTIME_CHECK(terminator, typeCode && typeCode->first == TypeCategory::Real);
92 int kind{typeCode->second};
93 switch (kind) {
94 // TODO: REAL (2 & 3)
95 case 4:
96 Generate<CppTypeFor<TypeCategory::Real, 4>, 24>(harvest);
97 break;
98 case 8:
99 Generate<CppTypeFor<TypeCategory::Real, 8>, 53>(harvest);
100 break;
101 #if LONG_DOUBLE == 80
102 case 10:
103 Generate<CppTypeFor<TypeCategory::Real, 10>, 64>(harvest);
104 break;
105 #elif LONG_DOUBLE == 128
106 case 16:
107 Generate<CppTypeFor<TypeCategory::Real, 16>, 113>(harvest);
108 break;
109 #endif
110 default:
111 terminator.Crash("RANDOM_NUMBER(): unsupported REAL kind %d", kind);
112 }
113 }
114
RTNAME(RandomSeedSize)115 void RTNAME(RandomSeedSize)(
116 const Descriptor &size, const char *source, int line) {
117 Terminator terminator{source, line};
118 auto typeCode{size.type().GetCategoryAndKind()};
119 RUNTIME_CHECK(terminator,
120 size.rank() == 0 && typeCode && typeCode->first == TypeCategory::Integer);
121 int kind{typeCode->second};
122 switch (kind) {
123 case 4:
124 *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1;
125 break;
126 case 8:
127 *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1;
128 break;
129 default:
130 terminator.Crash("RANDOM_SEED(SIZE=): bad kind %d\n", kind);
131 }
132 }
133
RTNAME(RandomSeedPut)134 void RTNAME(RandomSeedPut)(
135 const Descriptor &put, const char *source, int line) {
136 Terminator terminator{source, line};
137 auto typeCode{put.type().GetCategoryAndKind()};
138 RUNTIME_CHECK(terminator,
139 put.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer &&
140 put.GetDimension(0).Extent() >= 1);
141 int kind{typeCode->second};
142 GeneratedWord seed;
143 switch (kind) {
144 case 4:
145 seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>();
146 break;
147 case 8:
148 seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>();
149 break;
150 default:
151 terminator.Crash("RANDOM_SEED(PUT=): bad kind %d\n", kind);
152 }
153 {
154 CriticalSection critical{lock};
155 generator.seed(seed);
156 }
157 }
158
RTNAME(RandomSeedDefaultPut)159 void RTNAME(RandomSeedDefaultPut)() {
160 // TODO: should this be time &/or image dependent?
161 {
162 CriticalSection critical{lock};
163 generator.seed(0);
164 }
165 }
166
RTNAME(RandomSeedGet)167 void RTNAME(RandomSeedGet)(
168 const Descriptor &got, const char *source, int line) {
169 Terminator terminator{source, line};
170 auto typeCode{got.type().GetCategoryAndKind()};
171 RUNTIME_CHECK(terminator,
172 got.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer &&
173 got.GetDimension(0).Extent() >= 1);
174 int kind{typeCode->second};
175 GeneratedWord seed;
176 {
177 CriticalSection critical{lock};
178 seed = generator();
179 generator.seed(seed);
180 }
181 switch (kind) {
182 case 4:
183 *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed;
184 break;
185 case 8:
186 *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed;
187 break;
188 default:
189 terminator.Crash("RANDOM_SEED(GET=): bad kind %d\n", kind);
190 }
191 }
192 } // extern "C"
193 } // namespace Fortran::runtime
194