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>(PREC, 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(0); 102 } else { 103 #ifdef CLOCK_REALTIME 104 timespec ts; 105 clock_gettime(CLOCK_REALTIME, &ts); 106 generator.seed(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("not yet implemented: RANDOM_NUMBER(): REAL kind %d", kind); 146 } 147 148 void RTNAME(RandomSeedSize)( 149 const Descriptor *size, const char *source, int line) { 150 if (!size || !size->raw().base_addr) { 151 RTNAME(RandomSeedDefaultPut)(); 152 return; 153 } 154 Terminator terminator{source, line}; 155 auto typeCode{size->type().GetCategoryAndKind()}; 156 RUNTIME_CHECK(terminator, 157 size->rank() == 0 && typeCode && 158 typeCode->first == TypeCategory::Integer); 159 int kind{typeCode->second}; 160 switch (kind) { 161 case 4: 162 *size->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1; 163 break; 164 case 8: 165 *size->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1; 166 break; 167 default: 168 terminator.Crash( 169 "not yet implemented: RANDOM_SEED(SIZE=): kind %d\n", kind); 170 } 171 } 172 173 void RTNAME(RandomSeedPut)( 174 const Descriptor *put, const char *source, int line) { 175 if (!put || !put->raw().base_addr) { 176 RTNAME(RandomSeedDefaultPut)(); 177 return; 178 } 179 Terminator terminator{source, line}; 180 auto typeCode{put->type().GetCategoryAndKind()}; 181 RUNTIME_CHECK(terminator, 182 put->rank() == 1 && typeCode && 183 typeCode->first == TypeCategory::Integer && 184 put->GetDimension(0).Extent() >= 1); 185 int kind{typeCode->second}; 186 GeneratedWord seed; 187 switch (kind) { 188 case 4: 189 seed = *put->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>(); 190 break; 191 case 8: 192 seed = *put->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>(); 193 break; 194 default: 195 terminator.Crash("not yet implemented: RANDOM_SEED(PUT=): kind %d\n", kind); 196 } 197 { 198 CriticalSection critical{lock}; 199 generator.seed(seed); 200 nextValue = seed; 201 } 202 } 203 204 void RTNAME(RandomSeedDefaultPut)() { 205 // TODO: should this be time &/or image dependent? 206 { 207 CriticalSection critical{lock}; 208 generator.seed(0); 209 } 210 } 211 212 void RTNAME(RandomSeedGet)( 213 const Descriptor *get, const char *source, int line) { 214 if (!get || !get->raw().base_addr) { 215 RTNAME(RandomSeedDefaultPut)(); 216 return; 217 } 218 Terminator terminator{source, line}; 219 auto typeCode{get->type().GetCategoryAndKind()}; 220 RUNTIME_CHECK(terminator, 221 get->rank() == 1 && typeCode && 222 typeCode->first == TypeCategory::Integer && 223 get->GetDimension(0).Extent() >= 1); 224 int kind{typeCode->second}; 225 GeneratedWord seed; 226 { 227 CriticalSection critical{lock}; 228 seed = GetNextValue(); 229 nextValue = seed; 230 } 231 switch (kind) { 232 case 4: 233 *get->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed; 234 break; 235 case 8: 236 *get->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed; 237 break; 238 default: 239 terminator.Crash("not yet implemented: RANDOM_SEED(GET=): kind %d\n", kind); 240 } 241 } 242 243 void RTNAME(RandomSeed)(const Descriptor *size, const Descriptor *put, 244 const Descriptor *get, const char *source, int line) { 245 bool sizePresent = size && size->raw().base_addr; 246 bool putPresent = put && put->raw().base_addr; 247 bool getPresent = get && get->raw().base_addr; 248 if (sizePresent + putPresent + getPresent > 1) 249 Terminator{source, line}.Crash( 250 "RANDOM_SEED must have either 1 or no arguments"); 251 if (sizePresent) 252 RTNAME(RandomSeedSize)(size, source, line); 253 else if (putPresent) 254 RTNAME(RandomSeedPut)(put, source, line); 255 else if (getPresent) 256 RTNAME(RandomSeedGet)(get, source, line); 257 else 258 RTNAME(RandomSeedDefaultPut)(); 259 } 260 261 } // extern "C" 262 } // namespace Fortran::runtime 263