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