flang/runtime/random-templates.h

   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 Generate(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 } // namespace Fortran::runtime::random
  87
  88 #endif // FORTRAN_RUNTIME_RANDOM_TEMPLATES_H_