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