1 //===-- runtime/random-templates.h ------------------------------*- C++ -*-===// 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 #ifndef FORTRAN_RUNTIME_RANDOM_TEMPLATES_H_ 10 #define FORTRAN_RUNTIME_RANDOM_TEMPLATES_H_ 11 12 #include "lock.h" 13 #include "numeric-templates.h" 14 #include "flang/Common/optional.h" 15 #include "flang/Runtime/descriptor.h" 16 #include <algorithm> 17 #include <random> 18 19 namespace Fortran::runtime::random { 20 21 // Newer "Minimum standard", recommended by Park, Miller, and Stockmeyer in 22 // 1993. Same as C++17 std::minstd_rand, but explicitly instantiated for 23 // permanence. 24 using Generator = 25 std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>; 26 27 using GeneratedWord = typename Generator::result_type; 28 static constexpr std::uint64_t range{ 29 static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)}; 30 static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0}; 31 static constexpr int rangeBits{ 32 64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo}; 33 34 extern Lock lock; 35 extern Generator generator; 36 extern Fortran::common::optional<GeneratedWord> nextValue; 37 38 // Call only with lock held 39 static GeneratedWord GetNextValue() { 40 GeneratedWord result; 41 if (nextValue.has_value()) { 42 result = *nextValue; 43 nextValue.reset(); 44 } else { 45 result = generator(); 46 } 47 return result; 48 } 49 50 template <typename REAL, int PREC> 51 inline void GenerateReal(const Descriptor &harvest) { 52 static constexpr std::size_t minBits{ 53 std::max<std::size_t>(PREC, 8 * sizeof(GeneratedWord))}; 54 using Int = common::HostUnsignedIntType<minBits>; 55 static constexpr std::size_t words{ 56 static_cast<std::size_t>(PREC + rangeBits - 1) / rangeBits}; 57 std::size_t elements{harvest.Elements()}; 58 SubscriptValue at[maxRank]; 59 harvest.GetLowerBounds(at); 60 { 61 CriticalSection critical{lock}; 62 for (std::size_t j{0}; j < elements; ++j) { 63 while (true) { 64 Int fraction{GetNextValue()}; 65 if constexpr (words > 1) { 66 for (std::size_t k{1}; k < words; ++k) { 67 static constexpr auto rangeMask{ 68 (GeneratedWord{1} << rangeBits) - 1}; 69 GeneratedWord word{(GetNextValue() - generator.min()) & rangeMask}; 70 fraction = (fraction << rangeBits) | word; 71 } 72 } 73 fraction >>= words * rangeBits - PREC; 74 REAL next{ 75 LDEXPTy<REAL>::compute(static_cast<REAL>(fraction), -(PREC + 1))}; 76 if (next >= 0.0 && next < 1.0) { 77 *harvest.Element<REAL>(at) = next; 78 break; 79 } 80 } 81 harvest.IncrementSubscripts(at); 82 } 83 } 84 } 85 86 template <typename UINT> 87 inline void GenerateUnsigned(const Descriptor &harvest) { 88 static constexpr std::size_t words{ 89 (8 * sizeof(UINT) + rangeBits - 1) / rangeBits}; 90 std::size_t elements{harvest.Elements()}; 91 SubscriptValue at[maxRank]; 92 harvest.GetLowerBounds(at); 93 { 94 CriticalSection critical{lock}; 95 for (std::size_t j{0}; j < elements; ++j) { 96 UINT next{static_cast<UINT>(GetNextValue())}; 97 if constexpr (words > 1) { 98 for (std::size_t k{1}; k < words; ++k) { 99 next <<= rangeBits; 100 next |= GetNextValue(); 101 } 102 } 103 *harvest.Element<UINT>(at) = next; 104 harvest.IncrementSubscripts(at); 105 } 106 } 107 } 108 109 } // namespace Fortran::runtime::random 110 111 #endif // FORTRAN_RUNTIME_RANDOM_TEMPLATES_H_ 112