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 "random.h" 13 #include "cpp-type.h" 14 #include "descriptor.h" 15 #include "lock.h" 16 #include "flang/Common/leading-zero-bit-count.h" 17 #include "flang/Common/uint128.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("RANDOM_NUMBER(): unsupported REAL kind %d", kind); 112 } 113 } 114 115 void RTNAME(RandomSeedSize)( 116 const Descriptor &size, const char *source, int line) { 117 Terminator terminator{source, line}; 118 auto typeCode{size.type().GetCategoryAndKind()}; 119 RUNTIME_CHECK(terminator, 120 size.rank() == 0 && typeCode && typeCode->first == TypeCategory::Integer); 121 int kind{typeCode->second}; 122 switch (kind) { 123 case 4: 124 *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1; 125 break; 126 case 8: 127 *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1; 128 break; 129 default: 130 terminator.Crash("RANDOM_SEED(SIZE=): bad kind %d\n", kind); 131 } 132 } 133 134 void RTNAME(RandomSeedPut)( 135 const Descriptor &put, const char *source, int line) { 136 Terminator terminator{source, line}; 137 auto typeCode{put.type().GetCategoryAndKind()}; 138 RUNTIME_CHECK(terminator, 139 put.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer && 140 put.GetDimension(0).Extent() >= 1); 141 int kind{typeCode->second}; 142 GeneratedWord seed; 143 switch (kind) { 144 case 4: 145 seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>(); 146 break; 147 case 8: 148 seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>(); 149 break; 150 default: 151 terminator.Crash("RANDOM_SEED(PUT=): bad kind %d\n", kind); 152 } 153 { 154 CriticalSection critical{lock}; 155 generator.seed(seed); 156 } 157 } 158 159 void RTNAME(RandomSeedDefaultPut)() { 160 // TODO: should this be time &/or image dependent? 161 { 162 CriticalSection critical{lock}; 163 generator.seed(0); 164 } 165 } 166 167 void RTNAME(RandomSeedGet)( 168 const Descriptor &got, const char *source, int line) { 169 Terminator terminator{source, line}; 170 auto typeCode{got.type().GetCategoryAndKind()}; 171 RUNTIME_CHECK(terminator, 172 got.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer && 173 got.GetDimension(0).Extent() >= 1); 174 int kind{typeCode->second}; 175 GeneratedWord seed; 176 { 177 CriticalSection critical{lock}; 178 seed = generator(); 179 generator.seed(seed); 180 } 181 switch (kind) { 182 case 4: 183 *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed; 184 break; 185 case 8: 186 *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed; 187 break; 188 default: 189 terminator.Crash("RANDOM_SEED(GET=): bad kind %d\n", kind); 190 } 191 } 192 } // extern "C" 193 } // namespace Fortran::runtime 194