[flang] RANDOM_NUMBER, RANDOM_SEED, RANDOM_INIT in runtime

[klausler]

Summary

Add APIs, initial non-coarray implementations, and unit
tests for the intrinsic subroutines for pseudo-random
number generation.

[jeanPerier]

LGTM

flang/runtime/CMakeLists.txt

add_flang_library(FortranRuntime
   iostat.cpp
   io-api.cpp
   io-error.cpp
   io-stmt.cpp
   main.cpp
   memory.cpp
   misc-intrinsic.cpp
   numeric.cpp
+  random.cpp
   reduction.cpp
   stat.cpp
   stop.cpp
   terminator.cpp
   tools.cpp
   transformational.cpp
   type-code.cpp
   unit.cpp
   unit-map.cpp
 
   LINK_LIBS
   FortranDecimal
 )

flang/runtime/random.h

+//===-- runtime/random.h --------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+// Intrinsic subroutines RANDOM_INIT, RANDOM_NUMBER, and RANDOM_SEED.
+
+#include "entry-names.h"
+#include <cstdint>
+
+namespace Fortran::runtime {
+class Descriptor;
+extern "C" {
+
+void RTNAME(RandomInit)(bool repeatable, bool image_distinct);
+
+void RTNAME(RandomNumber)(
+    const Descriptor &harvest, const char *source, int line);
+
+// Subroutine RANDOM_SEED can be called with at most one of its optional
+// arguments; they each (plus the default case) map to these entry points.
+void RTNAME(RandomSeedSize)(const Descriptor &, const char *source, int line);
+void RTNAME(RandomSeedPut)(const Descriptor &, const char *source, int line);
+void RTNAME(RandomSeedGet)(const Descriptor &, const char *source, int line);
+void RTNAME(RandomSeedDefaultPut)();
+} // extern "C"
+} // namespace Fortran::runtime

flang/runtime/random.cpp

+//===-- runtime/random.cpp ------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+// Implements the intrinsic subroutines RANDOM_INIT, RANDOM_NUMBER, and
+// RANDOM_SEED.
+
+#include "random.h"
+#include "cpp-type.h"
+#include "descriptor.h"
+#include "lock.h"
+#include "flang/Common/leading-zero-bit-count.h"
+#include "flang/Common/uint128.h"
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <ctime>
+#include <limits>
+#include <memory>
+#include <random>
+
+namespace Fortran::runtime {
+
+// Newer "Minimum standard", recommended by Park, Miller, and Stockmeyer in
+// 1993. Same as C++17 std::minstd_rand, but explicitly instantiated for
+// permanence.
+using Generator =
+    std::linear_congruential_engine<std::uint_fast32_t, 48271, 0, 2147483647>;
+
+using GeneratedWord = typename Generator::result_type;
+static constexpr std::uint64_t range{
+    static_cast<std::uint64_t>(Generator::max() - Generator::min() + 1)};
+static constexpr bool rangeIsPowerOfTwo{(range & (range - 1)) == 0};
+static constexpr int rangeBits{
+    64 - common::LeadingZeroBitCount(range) - !rangeIsPowerOfTwo};
+
+static Lock lock;
+static Generator generator;
+
+template <typename REAL, int PREC>
+inline void Generate(const Descriptor &harvest) {
+  static constexpr std::size_t minBits{
+      std::max<std::size_t>(PREC, 8 * sizeof(GeneratedWord))};
+  using Int = common::HostUnsignedIntType<minBits>;
+  static constexpr std::size_t words{
+      static_cast<std::size_t>(PREC + rangeBits - 1) / rangeBits};
+  std::size_t elements{harvest.Elements()};
+  SubscriptValue at[maxRank];
+  harvest.GetLowerBounds(at);
+  {
+    CriticalSection critical{lock};
+    for (std::size_t j{0}; j < elements; ++j) {
+      Int fraction{generator()};
+      if constexpr (words > 1) {
+        for (std::size_t k{1}; k < words; ++k) {
+          static constexpr Int rangeMask{(Int{1} << rangeBits) - 1};
+          GeneratedWord word{(generator() - generator.min()) & rangeMask};
+          fraction = (fraction << rangeBits) | word;
+        }
+      }
+      fraction >>= words * rangeBits - PREC;
+      *harvest.Element<REAL>(at) =
+          std::ldexp(static_cast<REAL>(fraction), -(PREC + 1));
+      harvest.IncrementSubscripts(at);
+    }
+  }
+}
+
+extern "C" {
+
+void RTNAME(RandomInit)(bool repeatable, bool /*image_distinct*/) {
+  // TODO: multiple images and image_distinct: add image number
+  {
+    CriticalSection critical{lock};
+    if (repeatable) {
+      generator.seed(0);
+    } else {
+      generator.seed(std::time(nullptr));
+    }
+  }
+}
+
+void RTNAME(RandomNumber)(
+    const Descriptor &harvest, const char *source, int line) {
+  Terminator terminator{source, line};
+  auto typeCode{harvest.type().GetCategoryAndKind()};
+  RUNTIME_CHECK(terminator, typeCode && typeCode->first == TypeCategory::Real);
+  int kind{typeCode->second};
+  switch (kind) {
+  // TODO: REAL (2 & 3)
+  case 4:
+    Generate<CppTypeFor<TypeCategory::Real, 4>, 24>(harvest);
+    break;
+  case 8:
+    Generate<CppTypeFor<TypeCategory::Real, 8>, 53>(harvest);
+    break;
+#if LONG_DOUBLE == 80
+  case 10:
+    Generate<CppTypeFor<TypeCategory::Real, 10>, 64>(harvest);
+    break;
+#elif LONG_DOUBLE == 128
+  case 4:
+    Generate<CppTypeFor<TypeCategory::Real, 16>, 113>(harvest);
+    break;
+#endif
+  default:
+    terminator.Crash("RANDOM_NUMBER(): unsupported REAL kind %d", kind);
+  }
+}
+
+void RTNAME(RandomSeedSize)(
+    const Descriptor &size, const char *source, int line) {
+  Terminator terminator{source, line};
+  auto typeCode{size.type().GetCategoryAndKind()};
+  RUNTIME_CHECK(terminator,
+      size.rank() == 0 && typeCode && typeCode->first == TypeCategory::Integer);
+  int kind{typeCode->second};
+  switch (kind) {
+  case 4:
+    *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = 1;
+    break;
+  case 8:
+    *size.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = 1;
+    break;
+  default:
+    terminator.Crash("RANDOM_SEED(SIZE=): bad kind %d\n", kind);
+  }
+}
+
+void RTNAME(RandomSeedPut)(
+    const Descriptor &put, const char *source, int line) {
+  Terminator terminator{source, line};
+  auto typeCode{put.type().GetCategoryAndKind()};
+  RUNTIME_CHECK(terminator,
+      put.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer &&
+          put.GetDimension(0).Extent() >= 1);
+  int kind{typeCode->second};
+  GeneratedWord seed;
+  switch (kind) {
+  case 4:
+    seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>();
+    break;
+  case 8:
+    seed = *put.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>();
+    break;
+  default:
+    terminator.Crash("RANDOM_SEED(PUT=): bad kind %d\n", kind);
+  }
+  {
+    CriticalSection critical{lock};
+    generator.seed(seed);
+  }
+}
+
+void RTNAME(RandomSeedDefaultPut)() {
+  // TODO: should this be time &/or image dependent?
+  {
+    CriticalSection critical{lock};
+    generator.seed(0);
+  }
+}
+
+void RTNAME(RandomSeedGet)(
+    const Descriptor &got, const char *source, int line) {
+  Terminator terminator{source, line};
+  auto typeCode{got.type().GetCategoryAndKind()};
+  RUNTIME_CHECK(terminator,
+      got.rank() == 1 && typeCode && typeCode->first == TypeCategory::Integer &&
+          got.GetDimension(0).Extent() >= 1);
+  int kind{typeCode->second};
+  GeneratedWord seed;
+  {
+    CriticalSection critical{lock};
+    seed = generator();
+    generator.seed(seed);
+  }
+  switch (kind) {
+  case 4:
+    *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 4>>() = seed;
+    break;
+  case 8:
+    *got.OffsetElement<CppTypeFor<TypeCategory::Integer, 8>>() = seed;
+    break;
+  default:
+    terminator.Crash("RANDOM_SEED(GET=): bad kind %d\n", kind);
+  }
+}
+} // extern "C"
+} // namespace Fortran::runtime

flang/unittests/RuntimeGTest/CMakeLists.txt

 add_flang_unittest(FlangRuntimeTests
   CharacterTest.cpp
   CrashHandlerFixture.cpp
   MiscIntrinsic.cpp
   Numeric.cpp
   NumericalFormatTest.cpp
+  Random.cpp
   Reduction.cpp
   RuntimeCrashTest.cpp
 )
 
 target_link_libraries(FlangRuntimeTests
   PRIVATE
   FortranRuntime
 )

flang/unittests/RuntimeGTest/Random.cpp

+//===-- flang/unittests/RuntimeGTest/Random.cpp -----------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "../../runtime/random.h"
+#include "gtest/gtest.h"
+#include "../../runtime/descriptor.h"
+#include "../../runtime/type-code.h"
+#include <cmath>
+
+using namespace Fortran::runtime;
+
+TEST(RandomNumber, Real4) {
+  StaticDescriptor<1> statDesc;
+  Descriptor &harvest{statDesc.descriptor()};
+  static constexpr int n{10000};
+  float xs[n]{};
+  SubscriptValue extent[1]{n};
+  harvest.Establish(TypeCategory::Real, 4, xs, 1, extent);
+  RTNAME(RandomNumber)(harvest, __FILE__, __LINE__);
+  double sum{0};
+  for (int j{0}; j < n; ++j) {
+    sum += xs[j];
+  }
+  double mean{sum / n};
+  std::fprintf(stderr, "mean of %d random numbers: %g\n", n, mean);
+  EXPECT_GE(mean, 0.95 * 0.5); // mean of uniform dist [0..1] is of course 0.5
+  EXPECT_LE(mean, 1.05 * 0.5);
+  double sumsq{0};
+  for (int j{0}; j < n; ++j) {
+    double diff{xs[j] - mean};
+    sumsq += diff * diff;
+  }
+  double sdev{std::sqrt(sumsq / n)};
+  std::fprintf(stderr, "stddev of %d random numbers: %g\n", n, sdev);
+  double expect{1.0 / std::sqrt(12.0)}; // stddev of uniform dist [0..1]
+  EXPECT_GE(sdev, 0.95 * expect);
+  EXPECT_LT(sdev, 1.05 * expect);
+}
+
+TEST(RandomNumber, RandomSeed) {
+  StaticDescriptor<1> statDesc[2];
+  Descriptor &desc{statDesc[0].descriptor()};
+  std::int32_t n;
+  desc.Establish(TypeCategory::Integer, 4, &n, 0, nullptr);
+  RTNAME(RandomSeedSize)(desc, __FILE__, __LINE__);
+  EXPECT_EQ(n, 1);
+  SubscriptValue extent[1]{1};
+  desc.Establish(TypeCategory::Integer, 4, &n, 1, extent);
+  RTNAME(RandomSeedGet)(desc, __FILE__, __LINE__);
+  Descriptor &harvest{statDesc[1].descriptor()};
+  float x;
+  harvest.Establish(TypeCategory::Real, 4, &x, 1, extent);
+  RTNAME(RandomNumber)(harvest, __FILE__, __LINE__);
+  float got{x};
+  RTNAME(RandomSeedPut)(desc, __FILE__, __LINE__); // n from RandomSeedGet()
+  RTNAME(RandomNumber)(harvest, __FILE__, __LINE__);
+  EXPECT_EQ(x, got);
+}