diff --git a/flang/runtime/CMakeLists.txt b/flang/runtime/CMakeLists.txt --- a/flang/runtime/CMakeLists.txt +++ b/flang/runtime/CMakeLists.txt @@ -53,6 +53,7 @@ memory.cpp misc-intrinsic.cpp numeric.cpp + random.cpp reduction.cpp stat.cpp stop.cpp diff --git a/flang/runtime/random.h b/flang/runtime/random.h new file mode 100644 --- /dev/null +++ b/flang/runtime/random.h @@ -0,0 +1,30 @@ +//===-- 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 + +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 diff --git a/flang/runtime/random.cpp b/flang/runtime/random.cpp new file mode 100644 --- /dev/null +++ b/flang/runtime/random.cpp @@ -0,0 +1,193 @@ +//===-- 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 +#include +#include +#include +#include +#include +#include + +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; + +using GeneratedWord = typename Generator::result_type; +static constexpr std::uint64_t range{ + static_cast(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 +inline void Generate(const Descriptor &harvest) { + static constexpr std::size_t minBits{ + std::max(PREC, 8 * sizeof(GeneratedWord))}; + using Int = common::HostUnsignedIntType; + static constexpr std::size_t words{ + static_cast(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(at) = + std::ldexp(static_cast(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, 24>(harvest); + break; + case 8: + Generate, 53>(harvest); + break; +#if LONG_DOUBLE == 80 + case 10: + Generate, 64>(harvest); + break; +#elif LONG_DOUBLE == 128 + case 4: + Generate, 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>() = 1; + break; + case 8: + *size.OffsetElement>() = 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>(); + break; + case 8: + seed = *put.OffsetElement>(); + 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>() = seed; + break; + case 8: + *got.OffsetElement>() = seed; + break; + default: + terminator.Crash("RANDOM_SEED(GET=): bad kind %d\n", kind); + } +} +} // extern "C" +} // namespace Fortran::runtime diff --git a/flang/unittests/RuntimeGTest/CMakeLists.txt b/flang/unittests/RuntimeGTest/CMakeLists.txt --- a/flang/unittests/RuntimeGTest/CMakeLists.txt +++ b/flang/unittests/RuntimeGTest/CMakeLists.txt @@ -4,6 +4,7 @@ MiscIntrinsic.cpp Numeric.cpp NumericalFormatTest.cpp + Random.cpp Reduction.cpp RuntimeCrashTest.cpp ) diff --git a/flang/unittests/RuntimeGTest/Random.cpp b/flang/unittests/RuntimeGTest/Random.cpp new file mode 100644 --- /dev/null +++ b/flang/unittests/RuntimeGTest/Random.cpp @@ -0,0 +1,63 @@ +//===-- 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 + +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); +}