flang/runtime/random.cpp

   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 "terminator.h"
  15 #include "flang/Common/float128.h"
  16 #include "flang/Common/leading-zero-bit-count.h"
  17 #include "flang/Common/uint128.h"
  18 #include "flang/Runtime/cpp-type.h"
  19 #include "flang/Runtime/descriptor.h"
  20 #include <algorithm>
  21 #include <cmath>
  22 #include <cstdint>
  23 #include <limits>
  24 #include <memory>
  25 #include <random>
  26 #include <time.h>
  27
  28 namespace Fortran::runtime {
  29
  30 // Newer "Minimum standard", recommended by Park, Miller, and Stockmeyer in
  31 // 1993. Same as C++17 std::minstd_rand, but explicitly instantiated for
  32 // permanence.
  33 using Generator =
  34     std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>;
  35
  36 using GeneratedWord = typename Generator::result_type;
  37 static constexpr std::uint64_t range{
  38     static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)};
  39 static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0};
  40 static constexpr int rangeBits{
  41     64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo};
  42
  43 static Lock lock;
  44 static Generator generator;
  45 static std::optional<GeneratedWord> nextValue;
  46
  47 // Call only with lock held
  48 static GeneratedWord GetNextValue() {
  49   GeneratedWord result;
  50   if (nextValue.has_value()) {
  51     result = *nextValue;
  52     nextValue.reset();
  53   } else {
  54     result = generator();
  55   }
  56   return result;
  57 }
  58
  59 template <typename REAL, int PREC>
  60 inline void Generate(const Descriptor &harvest) {
  61   static constexpr std::size_t minBits{
  62       std::max<std::size_t>(PREC, 8 * sizeof(GeneratedWord))};
  63   using Int = common::HostUnsignedIntType<minBits>;
  64   static constexpr std::size_t words{
  65       static_cast<std::size_t>(PREC + rangeBits - 1) / rangeBits};
  66   std::size_t elements{harvest.Elements()};
  67   SubscriptValue at[maxRank];
  68   harvest.GetLowerBounds(at);
  69   {
  70     CriticalSection critical{lock};
  71     for (std::size_t j{0}; j < elements; ++j) {
  72       while (true) {
  73         Int fraction{GetNextValue()};
  74         if constexpr (words > 1) {
  75           for (std::size_t k{1}; k < words; ++k) {
  76             static constexpr auto rangeMask{
  77                 (GeneratedWord{1} << rangeBits) - 1};
  78             GeneratedWord word{(GetNextValue() - generator.min()) & rangeMask};
  79             fraction = (fraction << rangeBits) | word;
  80           }
  81         }
  82         fraction >>= words * rangeBits - PREC;
  83         REAL next{std::ldexp(static_cast<REAL>(fraction), -(PREC + 1))};
  84         if (next >= 0.0 && next < 1.0) {
  85           *harvest.Element<REAL>(at) = next;
  86           break;
  87         }
  88       }
  89       harvest.IncrementSubscripts(at);
  90     }
  91   }
  92 }
  93
  94 extern "C" {
  95
  96 void RTNAME(RandomInit)(bool repeatable, bool /*image_distinct*/) {
  97   // TODO: multiple images and image_distinct: add image number
  98   {
  99     CriticalSection critical{lock};
 100     if (repeatable) {
 101       generator.seed(0);
 102     } else {
 103 #ifdef CLOCK_REALTIME
 104       timespec ts;
 105       clock_gettime(CLOCK_REALTIME, &ts);
 106       generator.seed(ts.tv_sec & ts.tv_nsec);
 107 #else
 108       generator.seed(time(nullptr));
 109 #endif
 110     }
 111   }
 112 }
 113
 114 void RTNAME(RandomNumber)(
 115     const Descriptor &harvest, const char *source, int line) {
 116   Terminator terminator{source, line};
 117   auto typeCode{harvest.type().GetCategoryAndKind()};
 118   RUNTIME_CHECK(terminator, typeCode && typeCode->first == TypeCategory::Real);
 119   int kind{typeCode->second};
 120   switch (kind) {
 121   // TODO: REAL (2 & 3)
 122   case 4:
 123     Generate<CppTypeFor<TypeCategory::Real, 4>, 24>(harvest);
 124     return;
 125   case 8:
 126     Generate<CppTypeFor<TypeCategory::Real, 8>, 53>(harvest);
 127     return;
 128   case 10:
 129     if constexpr (HasCppTypeFor<TypeCategory::Real, 10>) {
 130 #if LDBL_MANT_DIG == 64
 131       Generate<CppTypeFor<TypeCategory::Real, 10>, 64>(harvest);
 132       return;
 133 #endif
 134     }
 135     break;
 136   case 16:
 137     if constexpr (HasCppTypeFor<TypeCategory::Real, 16>) {
 138 #if LDBL_MANT_DIG == 113
 139       Generate<CppTypeFor<TypeCategory::Real, 16>, 113>(harvest);
 140       return;
 141 #endif
 142     }
 143     break;
 144   }
 145   terminator.Crash("not yet implemented: RANDOM_NUMBER(): REAL kind %d", kind);
 146 }
 147
 148 void RTNAME(RandomSeedSize)(
 149     const Descriptor *size, const char *source, int line) {
 150   if (!size || !size->raw().base_addr) {
 151     RTNAME(RandomSeedDefaultPut)();
 152     return;
 153   }
 154   Terminator terminator{source, line};
 155   auto typeCode{size->type().GetCategoryAndKind()};
 156   RUNTIME_CHECK(terminator,
 157       size->rank() == 0 && typeCode &&
 158           typeCode->first == TypeCategory::Integer);
 159   int kind{typeCode->second};
 160   switch (kind) {
 161   case 4:
 162     *size->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1;
 163     break;
 164   case 8:
 165     *size->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1;
 166     break;
 167   default:
 168     terminator.Crash(
 169         "not yet implemented: RANDOM_SEED(SIZE=): kind %d\n", kind);
 170   }
 171 }
 172
 173 void RTNAME(RandomSeedPut)(
 174     const Descriptor *put, const char *source, int line) {
 175   if (!put || !put->raw().base_addr) {
 176     RTNAME(RandomSeedDefaultPut)();
 177     return;
 178   }
 179   Terminator terminator{source, line};
 180   auto typeCode{put->type().GetCategoryAndKind()};
 181   RUNTIME_CHECK(terminator,
 182       put->rank() == 1 && typeCode &&
 183           typeCode->first == TypeCategory::Integer &&
 184           put->GetDimension(0).Extent() >= 1);
 185   int kind{typeCode->second};
 186   GeneratedWord seed;
 187   switch (kind) {
 188   case 4:
 189     seed = *put->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>();
 190     break;
 191   case 8:
 192     seed = *put->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>();
 193     break;
 194   default:
 195     terminator.Crash("not yet implemented: RANDOM_SEED(PUT=): kind %d\n", kind);
 196   }
 197   {
 198     CriticalSection critical{lock};
 199     generator.seed(seed);
 200     nextValue = seed;
 201   }
 202 }
 203
 204 void RTNAME(RandomSeedDefaultPut)() {
 205   // TODO: should this be time &/or image dependent?
 206   {
 207     CriticalSection critical{lock};
 208     generator.seed(0);
 209   }
 210 }
 211
 212 void RTNAME(RandomSeedGet)(
 213     const Descriptor *get, const char *source, int line) {
 214   if (!get || !get->raw().base_addr) {
 215     RTNAME(RandomSeedDefaultPut)();
 216     return;
 217   }
 218   Terminator terminator{source, line};
 219   auto typeCode{get->type().GetCategoryAndKind()};
 220   RUNTIME_CHECK(terminator,
 221       get->rank() == 1 && typeCode &&
 222           typeCode->first == TypeCategory::Integer &&
 223           get->GetDimension(0).Extent() >= 1);
 224   int kind{typeCode->second};
 225   GeneratedWord seed;
 226   {
 227     CriticalSection critical{lock};
 228     seed = GetNextValue();
 229     nextValue = seed;
 230   }
 231   switch (kind) {
 232   case 4:
 233     *get->OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed;
 234     break;
 235   case 8:
 236     *get->OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed;
 237     break;
 238   default:
 239     terminator.Crash("not yet implemented: RANDOM_SEED(GET=): kind %d\n", kind);
 240   }
 241 }
 242
 243 void RTNAME(RandomSeed)(const Descriptor *size, const Descriptor *put,
 244     const Descriptor *get, const char *source, int line) {
 245   bool sizePresent = size && size->raw().base_addr;
 246   bool putPresent = put && put->raw().base_addr;
 247   bool getPresent = get && get->raw().base_addr;
 248   if (sizePresent + putPresent + getPresent > 1)
 249     Terminator{source, line}.Crash(
 250         "RANDOM_SEED must have either 1 or no arguments");
 251   if (sizePresent)
 252     RTNAME(RandomSeedSize)(size, source, line);
 253   else if (putPresent)
 254     RTNAME(RandomSeedPut)(put, source, line);
 255   else if (getPresent)
 256     RTNAME(RandomSeedGet)(get, source, line);
 257   else
 258     RTNAME(RandomSeedDefaultPut)();
 259 }
 260
 261 } // extern "C"
 262 } // namespace Fortran::runtime