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