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 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> 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 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 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 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 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 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 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