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 <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>
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       while (true) {
58         Int fraction{generator()};
59         if constexpr (words > 1) {
60           for (std::size_t k{1}; k < words; ++k) {
61             static constexpr auto rangeMask{
62                 (GeneratedWord{1} << rangeBits) - 1};
63             GeneratedWord word{(generator() - generator.min()) & rangeMask};
64             fraction = (fraction << rangeBits) | word;
65           }
66         }
67         fraction >>= words * rangeBits - PREC;
68         REAL next{std::ldexp(static_cast<REAL>(fraction), -(PREC + 1))};
69         if (next >= 0.0 && next < 1.0) {
70           *harvest.Element<REAL>(at) = next;
71           break;
72         }
73       }
74       harvest.IncrementSubscripts(at);
75     }
76   }
77 }
78 
79 extern "C" {
80 
81 void RTNAME(RandomInit)(bool repeatable, bool /*image_distinct*/) {
82   // TODO: multiple images and image_distinct: add image number
83   {
84     CriticalSection critical{lock};
85     if (repeatable) {
86       generator.seed(0);
87     } else {
88       generator.seed(std::time(nullptr));
89     }
90   }
91 }
92 
93 void RTNAME(RandomNumber)(
94     const Descriptor &harvest, const char *source, int line) {
95   Terminator terminator{source, line};
96   auto typeCode{harvest.type().GetCategoryAndKind()};
97   RUNTIME_CHECK(terminator, typeCode && typeCode->first == TypeCategory::Real);
98   int kind{typeCode->second};
99   switch (kind) {
100   // TODO: REAL (2 & 3)
101   case 4:
102     Generate<CppTypeFor<TypeCategory::Real, 4>, 24>(harvest);
103     break;
104   case 8:
105     Generate<CppTypeFor<TypeCategory::Real, 8>, 53>(harvest);
106     break;
107 #if LONG_DOUBLE == 80
108   case 10:
109     Generate<CppTypeFor<TypeCategory::Real, 10>, 64>(harvest);
110     break;
111 #elif LONG_DOUBLE == 128
112   case 16:
113     Generate<CppTypeFor<TypeCategory::Real, 16>, 113>(harvest);
114     break;
115 #endif
116   default:
117     terminator.Crash(
118         "not yet implemented: RANDOM_NUMBER(): REAL kind %d", kind);
119   }
120 }
121 
122 void RTNAME(RandomSeedSize)(
123     const Descriptor &size, const char *source, int line) {
124   Terminator terminator{source, line};
125   auto typeCode{size.type().GetCategoryAndKind()};
126   RUNTIME_CHECK(terminator,
127       size.rank() == 0 && typeCode && typeCode->first == TypeCategory::Integer);
128   int kind{typeCode->second};
129   switch (kind) {
130   case 4:
131     *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1;
132     break;
133   case 8:
134     *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1;
135     break;
136   default:
137     terminator.Crash(
138         "not yet implemented: RANDOM_SEED(SIZE=): kind %d\n", kind);
139   }
140 }
141 
142 void RTNAME(RandomSeedPut)(
143     const Descriptor &put, const char *source, int line) {
144   Terminator terminator{source, line};
145   auto typeCode{put.type().GetCategoryAndKind()};
146   RUNTIME_CHECK(terminator,
147       put.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer &&
148           put.GetDimension(0).Extent() >= 1);
149   int kind{typeCode->second};
150   GeneratedWord seed;
151   switch (kind) {
152   case 4:
153     seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>();
154     break;
155   case 8:
156     seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>();
157     break;
158   default:
159     terminator.Crash("not yet implemented: RANDOM_SEED(PUT=): kind %d\n", kind);
160   }
161   {
162     CriticalSection critical{lock};
163     generator.seed(seed);
164   }
165 }
166 
167 void RTNAME(RandomSeedDefaultPut)() {
168   // TODO: should this be time &/or image dependent?
169   {
170     CriticalSection critical{lock};
171     generator.seed(0);
172   }
173 }
174 
175 void RTNAME(RandomSeedGet)(
176     const Descriptor &got, const char *source, int line) {
177   Terminator terminator{source, line};
178   auto typeCode{got.type().GetCategoryAndKind()};
179   RUNTIME_CHECK(terminator,
180       got.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer &&
181           got.GetDimension(0).Extent() >= 1);
182   int kind{typeCode->second};
183   GeneratedWord seed;
184   {
185     CriticalSection critical{lock};
186     seed = generator();
187     generator.seed(seed);
188   }
189   switch (kind) {
190   case 4:
191     *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed;
192     break;
193   case 8:
194     *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed;
195     break;
196   default:
197     terminator.Crash("not yet implemented: RANDOM_SEED(GET=): kind %d\n", kind);
198   }
199 }
200 } // extern "C"
201 } // namespace Fortran::runtime
202