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
19namespace 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.
24using Generator =
25 std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>;
26
27using GeneratedWord = typename Generator::result_type;
28static constexpr std::uint64_t range{
29 static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)};
30static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0};
31static constexpr int rangeBits{
32 64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo};
33
34extern Lock lock;
35extern Generator generator;
36extern Fortran::common::optional<GeneratedWord> nextValue;
37
38// Call only with lock held
39static 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
50template <typename REAL, int PREC>
51inline void Generate(const Descriptor &harvest) {
52 static constexpr std::size_t minBits{
53 std::max<std::size_t>(a: PREC, b: 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} // namespace Fortran::runtime::random
87
88#endif // FORTRAN_RUNTIME_RANDOM_TEMPLATES_H_
89

source code of flang/runtime/random-templates.h