[flang] Rework host runtime folding and enable REAL(2) folding with it.

[jeanPerier]

Summary

• Rework the host runtime table so that it is constexpr to avoid having to construct it and to store/propagate it.
• Make the interface simpler (remove many templates and a file)
• Enable 16bits float folding using 32bits float host runtime
• Move StaticMultimapView into its own header to use it for host folding

[klausler]

There's a lot of work here and it's quite well done; thank you.

I suggest that you add meinersbur as a reviewer to check that this code doesn't run afoul of some MSVC bug.

flang/include/flang/Common/static-multimap-view.h

43

The argument could be const std::array<V,N> &array.

48

This could be a constexpr check in the constructor.

51

We use GetRange naming in Common.

Could GetRange be private?

86

C++20 has constexpr versions of std::equal_range, std::lower_bound, std::binary_search, &c.

If you could use them now, would this new class be necessary?

If we will eventually use them, maybe what's needed here would be just a place-holding constexpr function or two, possibly conditional on the C++ language standard version.

flang/include/flang/Evaluate/common.h

11–12

Does removing this member from the folding context make them cheap to construct again?

flang/lib/Evaluate/fold-implementation.h

88

return std::nullopt; is a little more clear, I think.

flang/lib/Evaluate/intrinsics-library.cpp

142

Please capitalize the names of these functions; thanks.

flang/lib/Lower/IntrinsicCall.cpp

507–508

I know that this isn't part of this change, but maybe this argument should be a const std::array<> &.

flang/unittests/Evaluate/folding.cpp

55–56

I wonder if this constant should be checked; I worry about C++ compilation on a target machine that likes to flush subnormals to zero. Maybe a hexadecimal floating point constant should be used, or even a reinterpreted 0x00400000.

(The value should be 5.877471754111438E-39 if you want a double constant that converts to and back from the desired value; the value you have above is slightly too large and converts to a double with a lower bit set that's lost in the conversion to float in most rounding modes, but why take chances.)

[PeteSteinfeld]

Thanks for working on this, Jean. I have to confess that I don't understand most of the code. But it all built and tested correctly for me. It would be nice to have some constant folding tests that fold REAL(2) operands.

Where do things stand on 80 bit and 128 bit REALs?

[jeanPerier]

• Simplify StaticMultimapView (use std::equal_range and removed unused member functions)
• Add REAL(2) and REAL(3) folding test
• Enable X87 float folding with IEEE 128 bits float runtime if available on the host.
• Use RealKindForPrecision instead of TypeOf
• Add a compile time verfication that StaticMultimapView is sorted
• Capitalize some functions. std::nullopt instead of {}.
• Use raw subnormal in flushing test instead of decimal input.

[jeanPerier]

Remove static_assert from StaticMultimapView

The last update added a compile time assert that the array
is sorted before constructing a StaticMultimapView. It worked well with clang
but was actually most likely illegal C++ and failed with g++ since there is
no way to guarantee that the content of the array passed in template argument
is constexpr (only its address is guaranteed to be so). Clang++ seems to be
clever enough to propagate the constexpr aspect of the array address it gets
in template argument, but this is not standard C++.
Instead Verify() has to be called next to every instantiation of the map
where it is visible to the compiler that the map was built at compile time...

Also avoid the early return since constexpr functions should only have one return.

[jeanPerier]

Thanks for the reviews, I have updated the patch where I could.

flang/include/flang/Common/static-multimap-view.h

43

I do not know a way to do that without having the tables declared as std::array which I find annoying because you cannot nicely delegate the array length to the compiler (e.g. T x[]{...} vs std::array<T, /* need to compute/ guess that */> x{...}).

48

Thanks, I tried hard to make that happen in the ctor/a static Create helper, but there is no way I could find to have it constexpr check that way since the constexpr aspect does not actually guarantee that this will be called with a constexpr array at compile time, so static_assert are rejected.
Instead I added checks after each instantiation. That's still better than nothing and actually caught a bessel_y0/y1 not sorted correctly in pgmath header.

51

GetRange was removed.

86

What is mainly missing for me here is a constexpr std::is_storted (also C++20), as well as something like std::span (C++20) that (from what I understand) can "erase" the lenght from the base container type and still allow constexpr iteration on it.
GetRange used to be called in constexpr in lowering to get the runtime but that is not the case anymore, so I already switched to std::equal_range and remove some other member functions that removed ~20 lines of boilerplate.

flang/unittests/Evaluate/folding.cpp

55–56

Thanks for the tip, I used the raw evaluate::Real constructor to avoid any issue with C++ compilers here.

[jeanPerier]

In D88981#2321900, @PeteSteinfeld wrote:

Where do things stand on 80 bit and 128 bit REALs?

I have updated the patch to allow using IEEE 128bits to fold 80 bits (X87) floats when IEEE 128 bits is available but not X87.
However 80/128 bits folding that require host runtime is currently not expected to work on all host at all (hence I did not add tests yet).
The folding code is designed to take advantage of when 80/128bits can be safely folded using long double without assuming too much about the platform. It is also designed to allow platform where 80/128bits cannot be safely folded to still compile/run flang normally and only refuse to compile programs that require such host folding. More precisely:

• 80 bits folding works if the host long double is the X87 or IEEE 128bits float (for this last case we would convert 80bits to 128bits before folding).
• 128 bits folding works if the host long double is the IEEE 128 bits.
• When the host long double type does not allow folding 80bits/128bits, Fold() will return the original expressions, and a warning message (non fatal) saying " intrinsic acos(real(kind=16)) cannot be folded on host". If the expressions must be a constant, compilation should then failed complaining it did not get a constant after Fold() calls. Otherwise, compilation will continue normally (leaving the expression unfolded since we do not have to, but warning the user that we could not fold it due to host limitation).

Hosts that needs to fold 80bits/128bits will need to modify host.h to map some non standard types (e.g. __float128) in host.h provided the c++ standard math library works with them (or they will need to provide an ad-hoc runtime with it in intrinsics-library.cpp).

I tested that 80 bits folding works OK on X86-64 linux. I expect 80bits and 128bits folding to work on Aarch64, but did not test it. On Power8/9, I think that depends on how the OS is compiled there. By default long double maps to IBM double double format, but the OS/libc/gcc... can also be compiled to map long double to IEEE 128bits in which case folding 80bits and 128bits should work. As you see, all this is very dependent on the platform.

Note that this only applies to folding intrinsics for which there is no abstract implementation in include/flang/Evaluate/real.h and for which we have to resort to host runtime, otherwise all REALs can be folded regardless of the platform.

[jeanPerier]

@Meinersbur, I have added you as a reviewer since this patch is touching a lot of template/constexpr features so it is likely it will hit MSVC errors. If you have the opportunity to test it ahead that would be great, otherwise beware that it may lead to MSVC regressions.

[Meinersbur]

The patch compiles successfully with msvc (with a patch to trunk that I still need to upload a patch for).

[jeanPerier]

The patch compiles successfully with msvc (with a patch to trunk that I still need to upload a patch for).

Thanks for testing this @Meinersbur !

flang/include/flang/Evaluate/common.h

11–12

Yes, FoldingContext are 100 times cheaper to construct according to my measurements. This improves fcvs f18 -fparse-only time by on 12% on average.

FoldingContext ctor 100x speedup
With f18 compiled with gcc 8.3 in release mode on an Intel Xeon Gold 6148, I measured 0.05ms per FoldingContext construction before vs 0.00005ms with this patch (average of 10000 ctor calls in one run. I reproduced runs 10times and got stable results). Measurement were done by instrumenting the code (https://github.com/jeanPerier/llvm-project/commit/f511284b54805aa314c1316f9143d0d0cbaa522d).

Given FoldingContext are constructed for every function call check when an explicit interface that can translate in x4 speed-up one time f18 -fparse-only on carefully designed tests like:

 real, parameter :: x = 0.5

 ! Each following line semantic analysis end-up in 3 FoldingContext ctor call 
 real, parameter :: y1 = acos(x)
 real, parameter :: y2 = acos(x)
 ! ... repeated 9997 times 
 real, parameter :: y10000 = acos(x)
end

I measured 2s before vs 0.5s with this patch (time f18 -fparse-only real time).

Host folding 1.2x slowdown
However, there is a 20% time penalty with this patch per fold with host runtime (most likely due to the added encapsulation/decapsulation of Scalar to/from Expr<SomeExpr> in the folder). I measured the time spent in Evaluate/fold-real.cpp FoldIntrinsicFunction on the test file above. We spent 1.3usec per fold before vs 1.6usec with this patch (average of the 10000 folds, repeated 10 times). Given a for this is at the usec level, it is negligible on scalar fold since we create 3 FoldingContext per expressions. For array expressions, that can lead to overall slowdown in the compilation (that will never be bigger than 20%). For instance I could measure a 1% overall slow-down in a program folding acos( a_10000_element_array) (93ms before vs 94 now).

Conclusion: 12% overall parsing+semantics speed-up on real code
Regarding fcvs time f18 -fparse-only fm*.f real time went from 4.3s to 3.8s (ten run average). So this has a visible impact on real code.
Since scalar folding is much more widespread than huge array folding, the patch seems a win to me.

flang/include/flang/Common/static-multimap-view.h

+//===-- include/flang/Common/static-multimap-view.h -------------*- 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef FORTRAN_COMMON_STATIC_MULTIMAP_VIEW_H_
+#define FORTRAN_COMMON_STATIC_MULTIMAP_VIEW_H_
+#include <algorithm>
+#include <utility>
+
+/// StaticMultimapView is a constexpr friendly multimap implementation over
+/// sorted constexpr arrays. As the View name suggests, it does not duplicate
+/// the sorted array but only brings range and search concepts over it. It
+/// mainly erases the array size from the type and ensures the array is sorted
+/// at compile time. When C++20 brings std::span and constexpr std::is_sorted,
+/// this can most likely be replaced by those.
+
+namespace Fortran::common {
+
+template <typename V> class StaticMultimapView {
+public:
+  using Key = typename V::Key;
+  using const_iterator = const V *;
+
+  constexpr const_iterator begin() const { return begin_; }
+  constexpr const_iterator end() const { return end_; }
+  // Be sure to static_assert(map.Verify(), "must be sorted"); for
+  // every instance constexpr created. Sadly this cannot be done in
+  // the ctor since there is no way to know whether the ctor is actually
+  // called at compile time or not.
+  template <std::size_t N>
+  constexpr StaticMultimapView(const V (&array)[N])
+      : begin_{&array[0]}, end_{&array[0] + N} {}
+
+  // std::equal_range will be constexpr in C++20 only, so far there is actually
+  // no need for equal_range to be constexpr anyway.
+  std::pair<const_iterator, const_iterator> equal_range(const Key &key) const {
+    return std::equal_range(begin_, end_, key);
+  }
+
+  // Check that the array is sorted. This used to assert at compile time that
+  // the array is indeed sorted. When C++20 is required for flang,
+  // std::is_sorted can be used here since it will be constexpr.
+  constexpr bool Verify() const {
+    const V *lastSeen{begin_};
+    bool isSorted{true};
+    for (const auto *x{begin_}; x != end_; ++x) {
+      isSorted &= lastSeen->key <= x->key;
+      lastSeen = x;
+    }
+    return isSorted;
+  }
+
+private:
+  const_iterator begin_{nullptr};
+  const_iterator end_{nullptr};
+};
+} // namespace Fortran::common
+#endif // FORTRAN_COMMON_STATIC_MULTIMAP_VIEW_H_

flang/include/flang/Evaluate/common.h

 //===-- include/flang/Evaluate/common.h -------------------------*- 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
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef FORTRAN_EVALUATE_COMMON_H_
 #define FORTRAN_EVALUATE_COMMON_H_
 
-#include "intrinsics-library.h"
 #include "flang/Common/Fortran.h"
 #include "flang/Common/default-kinds.h"
 #include "flang/Common/enum-set.h"
 #include "flang/Common/idioms.h"
 #include "flang/Common/indirection.h"
 #include "flang/Common/restorer.h"
 #include "flang/Parser/char-block.h"
 #include "flang/Parser/message.h"
 #include <cinttypes>
public:
   const parser::ContextualMessages &messages() const { return messages_; }
   const common::IntrinsicTypeDefaultKinds &defaults() const {
     return defaults_;
   }
   Rounding rounding() const { return rounding_; }
   bool flushSubnormalsToZero() const { return flushSubnormalsToZero_; }
   bool bigEndian() const { return bigEndian_; }
   const semantics::DerivedTypeSpec *pdtInstance() const { return pdtInstance_; }
-  const HostIntrinsicProceduresLibrary &hostIntrinsicsLibrary() const {
-    return hostIntrinsicsLibrary_;
-  }
   const evaluate::IntrinsicProcTable &intrinsics() const { return intrinsics_; }
 
   ConstantSubscript &StartImpliedDo(parser::CharBlock, ConstantSubscript = 1);
   std::optional<ConstantSubscript> GetImpliedDo(parser::CharBlock) const;
   void EndImpliedDo(parser::CharBlock);
 
   std::map<parser::CharBlock, ConstantSubscript> &impliedDos() {
     return impliedDos_;
   }
 
   common::Restorer<const semantics::DerivedTypeSpec *> WithPDTInstance(
       const semantics::DerivedTypeSpec &spec) {
     return common::ScopedSet(pdtInstance_, &spec);
   }
 
 private:
   parser::ContextualMessages messages_;
   const common::IntrinsicTypeDefaultKinds &defaults_;
   const IntrinsicProcTable &intrinsics_;
   Rounding rounding_{defaultRounding};
   bool flushSubnormalsToZero_{false};
   bool bigEndian_{false};
   const semantics::DerivedTypeSpec *pdtInstance_{nullptr};
   std::map<parser::CharBlock, ConstantSubscript> impliedDos_;
-  HostIntrinsicProceduresLibrary hostIntrinsicsLibrary_;
 };
 
 void RealFlagWarnings(FoldingContext &, const RealFlags &, const char *op);
 } // namespace Fortran::evaluate
 #endif // FORTRAN_EVALUATE_COMMON_H_

flang/include/flang/Evaluate/intrinsics-library.h

 //===-- include/flang/Evaluate/intrinsics-library.h -------------*- 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
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef FORTRAN_EVALUATE_INTRINSICS_LIBRARY_H_
 #define FORTRAN_EVALUATE_INTRINSICS_LIBRARY_H_
 
 // Defines structures to be used in F18 for folding intrinsic function with host
-// runtime libraries. To avoid unnecessary header circular dependencies, the
-// actual implementation of the templatized member function are defined in
-// intrinsics-library-templates.h The header at hand is meant to be included by
-// files that need to define intrinsic runtime data structure but that do not
-// use them directly. To actually use the runtime data structures,
-// intrinsics-library-templates.h must be included.
+// runtime libraries.
 
 #include <functional>
-#include <map>
 #include <optional>
 #include <string>
 #include <vector>
 
 namespace Fortran::evaluate {
 class FoldingContext;
-
-using TypeCode = unsigned char;
-
-template <typename TR, typename... TA> using FuncPointer = TR (*)(TA...);
-// This specific type signature prevents GCC complaining about function casts.
-using GenericFunctionPointer = void (*)(void);
-
-enum class PassBy { Ref, Val };
-template <typename TA, PassBy Pass = PassBy::Ref> struct ArgumentInfo {
-  using Type = TA;
-  static constexpr PassBy pass{Pass};
-};
-
-template <typename TR, typename... ArgInfo> struct Signature {
-  // Note valid template argument are of form
-  //<TR, ArgumentInfo<TA, PassBy>...> where TA and TR belong to RuntimeTypes.
-  // RuntimeTypes is a type union defined in intrinsics-library-templates.h to
-  // avoid circular dependencies. Argument of type void cannot be passed by
-  // value. So far TR cannot be a pointer.
-  const std::string name;
-};
-
-struct IntrinsicProcedureRuntimeDescription {
-  const std::string name;
-  const TypeCode returnType;
-  const std::vector<TypeCode> argumentsType;
-  const std::vector<PassBy> argumentsPassedBy;
-  const bool isElemental;
-  const GenericFunctionPointer callable;
-  // Construct from description using host independent types (RuntimeTypes)
-  template <typename TR, typename... ArgInfo>
-  IntrinsicProcedureRuntimeDescription(
-      const Signature<TR, ArgInfo...> &signature, bool isElemental = false);
-};
-
-// HostRuntimeIntrinsicProcedure allows host runtime function to be called for
-// constant folding.
-struct HostRuntimeIntrinsicProcedure : IntrinsicProcedureRuntimeDescription {
-  // Construct from runtime pointer with host types (float, double....)
-  template <typename HostTR, typename... HostTA>
-  HostRuntimeIntrinsicProcedure(const std::string &name,
-      FuncPointer<HostTR, HostTA...> func, bool isElemental = false);
-  HostRuntimeIntrinsicProcedure(
-      const IntrinsicProcedureRuntimeDescription &rteProc,
-      GenericFunctionPointer handle)
-      : IntrinsicProcedureRuntimeDescription{rteProc}, handle{handle} {}
-  GenericFunctionPointer handle;
-};
-
-// Defines a wrapper type that indirects calls to host runtime functions.
-// Valid ConstantContainer are Scalar (only for elementals) and Constant.
-template <template <typename> typename ConstantContainer, typename TR,
-    typename... TA>
-using HostProcedureWrapper = std::function<ConstantContainer<TR>(
-    FoldingContext &, ConstantContainer<TA>...)>;
-
-// HostIntrinsicProceduresLibrary is a data structure that holds
-// HostRuntimeIntrinsicProcedure elements. It is meant for constant folding.
-// When queried for an intrinsic procedure, it can return a callable object that
-// implements this intrinsic if a host runtime function pointer for this
-// intrinsic was added to its data structure.
-class HostIntrinsicProceduresLibrary {
-public:
-  HostIntrinsicProceduresLibrary();
-  void AddProcedure(HostRuntimeIntrinsicProcedure &&sym) {
-    const std::string name{sym.name};
-    procedures_.insert(std::make_pair(name, std::move(sym)));
-  }
-  bool HasEquivalentProcedure(
-      const IntrinsicProcedureRuntimeDescription &sym) const;
-  template <template <typename> typename ConstantContainer, typename TR,
-      typename... TA>
-  std::optional<HostProcedureWrapper<ConstantContainer, TR, TA...>>
-  GetHostProcedureWrapper(const std::string &name) const;
-
-private:
-  std::multimap<std::string, const HostRuntimeIntrinsicProcedure> procedures_;
-};
-
+class DynamicType;
+struct SomeType;
+template <typename> class Expr;
+
+// Define a callable type that is used to fold scalar intrinsic function using
+// host runtime. These callables are responsible for the conversions between
+// host types and Fortran abstract types (Scalar<T>). They also deal with
+// floating point environment (To set it up to match the Fortran compiling
+// options and to clean it up after the call). Floating point errors are
+// reported to the FoldingContext. For 16bits float types, 32bits float host
+// runtime plus conversions may be used to build the host wrappers if no 16bits
+// runtime is available. IEEE 128bits float may also be used for x87 float.
+// Potential conversion overflows are reported by the HostRuntimeWrapper in the
+// FoldingContext.
+using HostRuntimeWrapper = std::function<Expr<SomeType>(
+    FoldingContext &, std::vector<Expr<SomeType>> &&)>;
+
+// Returns the folder using host runtime given the intrinsic function name,
+// result and argument types. Nullopt if no host runtime is available for such
+// intrinsic function.
+std::optional<HostRuntimeWrapper> GetHostRuntimeWrapper(const std::string &name,
+    DynamicType resultType, const std::vector<DynamicType> &argTypes);
 } // namespace Fortran::evaluate
 #endif // FORTRAN_EVALUATE_INTRINSICS_LIBRARY_H_

flang/lib/Evaluate/fold-complex.cpp

Expr<Type<TypeCategory::Complex, KIND>> FoldIntrinsicFunction(
   ActualArguments &args{funcRef.arguments()};
   auto *intrinsic{std::get_if<SpecificIntrinsic>(&funcRef.proc().u)};
   CHECK(intrinsic);
   std::string name{intrinsic->name};
   if (name == "acos" || name == "acosh" || name == "asin" || name == "asinh" ||
       name == "atan" || name == "atanh" || name == "cos" || name == "cosh" ||
       name == "exp" || name == "log" || name == "sin" || name == "sinh" ||
       name == "sqrt" || name == "tan" || name == "tanh") {
-    if (auto callable{context.hostIntrinsicsLibrary()
-                          .GetHostProcedureWrapper<Scalar, T, T>(name)}) {
+    if (auto callable{GetHostRuntimeWrapper<T, T>(name)}) {
       return FoldElementalIntrinsic<T, T>(
           context, std::move(funcRef), *callable);
     } else {
       context.messages().Say(
           "%s(complex(kind=%d)) cannot be folded on host"_en_US, name, KIND);
     }
   } else if (name == "conjg") {
     return FoldElementalIntrinsic<T, T>(

flang/lib/Evaluate/fold-implementation.h

 //===-- lib/Evaluate/fold-implementation.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
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef FORTRAN_EVALUATE_FOLD_IMPLEMENTATION_H_
 #define FORTRAN_EVALUATE_FOLD_IMPLEMENTATION_H_
 
 #include "character.h"
 #include "host.h"
 #include "int-power.h"
-#include "intrinsics-library-templates.h"
 #include "flang/Common/indirection.h"
 #include "flang/Common/template.h"
 #include "flang/Common/unwrap.h"
 #include "flang/Evaluate/characteristics.h"
 #include "flang/Evaluate/common.h"
 #include "flang/Evaluate/constant.h"
 #include "flang/Evaluate/expression.h"
 #include "flang/Evaluate/fold.h"
 #include "flang/Evaluate/formatting.h"
+#include "flang/Evaluate/intrinsics-library.h"
 #include "flang/Evaluate/intrinsics.h"
 #include "flang/Evaluate/shape.h"
 #include "flang/Evaluate/tools.h"
 #include "flang/Evaluate/traverse.h"
 #include "flang/Evaluate/type.h"
 #include "flang/Parser/message.h"
 #include "flang/Semantics/scope.h"
 #include "flang/Semantics/symbol.h"

 
 private:
   FoldingContext &context_;
 };
 
 std::optional<Constant<SubscriptInteger>> GetConstantSubscript(
     FoldingContext &, Subscript &, const NamedEntity &, int dim);
 
+// Helper to use host runtime on scalars for folding.
+template <typename TR, typename... TA>
+std::optional<std::function<Scalar<TR>(FoldingContext &, Scalar<TA>...)>>
+GetHostRuntimeWrapper(const std::string &name) {
+  std::vector<DynamicType> argTypes{TA{}.GetType()...};
+  if (auto hostWrapper{GetHostRuntimeWrapper(name, TR{}.GetType(), argTypes)}) {
+    return [hostWrapper](
+               FoldingContext &context, Scalar<TA>... args) -> Scalar<TR> {
+      std::vector<Expr<SomeType>> genericArgs{
+          AsGenericExpr(Constant<TA>{args})...};
+      return GetScalarConstantValue<TR>(
+          (*hostWrapper)(context, std::move(genericArgs)))
+          .value();
+    };
+  }
+  return std::nullopt;
+}
+
 // FoldOperation() rewrites expression tree nodes.
 // If there is any possibility that the rewritten node will
 // not have the same representation type, the result of
 // FoldOperation() will be packaged in an Expr<> of the same
 // specific type.
 
 // no-op base case
 template <typename A>
if constexpr (T::category == TypeCategory::Integer) {
       } else if (power.overflow) {
         context.messages().Say("INTEGER(%d) power overflowed"_en_US, T::kind);
       } else if (power.zeroToZero) {
         context.messages().Say(
             "INTEGER(%d) 0**0 is not defined"_en_US, T::kind);
       }
       return Expr<T>{Constant<T>{power.power}};
     } else {
-      if (auto callable{context.hostIntrinsicsLibrary()
-                            .GetHostProcedureWrapper<Scalar, T, T, T>("pow")}) {
+      if (auto callable{GetHostRuntimeWrapper<T, T, T>("pow")}) {
         return Expr<T>{
             Constant<T>{(*callable)(context, folded->first, folded->second)}};
       } else {
         context.messages().Say(
             "Power for %s cannot be folded on host"_en_US, T{}.AsFortran());
       }
     }
   }

flang/lib/Evaluate/fold-real.cpp

if (name == "acos" || name == "acosh" || name == "asin" || name == "asinh" ||
       (name == "atan" && args.size() == 1) || name == "atanh" ||
       name == "bessel_j0" || name == "bessel_j1" || name == "bessel_y0" ||
       name == "bessel_y1" || name == "cos" || name == "cosh" || name == "erf" ||
       name == "erfc" || name == "erfc_scaled" || name == "exp" ||
       name == "gamma" || name == "log" || name == "log10" ||
       name == "log_gamma" || name == "sin" || name == "sinh" ||
       name == "sqrt" || name == "tan" || name == "tanh") {
     CHECK(args.size() == 1);
-    if (auto callable{context.hostIntrinsicsLibrary()
-                          .GetHostProcedureWrapper<Scalar, T, T>(name)}) {
+    if (auto callable{GetHostRuntimeWrapper<T, T>(name)}) {
       return FoldElementalIntrinsic<T, T>(
           context, std::move(funcRef), *callable);
     } else {
       context.messages().Say(
           "%s(real(kind=%d)) cannot be folded on host"_en_US, name, KIND);
     }
   } else if (name == "amax0" || name == "amin0" || name == "amin1" ||
       name == "amax1" || name == "dmin1" || name == "dmax1") {
     return RewriteSpecificMINorMAX(context, std::move(funcRef));
   } else if (name == "atan" || name == "atan2" || name == "hypot" ||
       name == "mod") {
     std::string localName{name == "atan" ? "atan2" : name};
     CHECK(args.size() == 2);
-    if (auto callable{
-            context.hostIntrinsicsLibrary()
-                .GetHostProcedureWrapper<Scalar, T, T, T>(localName)}) {
+    if (auto callable{GetHostRuntimeWrapper<T, T, T>(localName)}) {
       return FoldElementalIntrinsic<T, T, T>(
           context, std::move(funcRef), *callable);
     } else {
       context.messages().Say(
           "%s(real(kind=%d), real(kind%d)) cannot be folded on host"_en_US,
           name, KIND, KIND);
     }
   } else if (name == "bessel_jn" || name == "bessel_yn") {
     if (args.size() == 2) { // elemental
       // runtime functions use int arg
       using Int4 = Type<TypeCategory::Integer, 4>;
-      if (auto callable{
-              context.hostIntrinsicsLibrary()
-                  .GetHostProcedureWrapper<Scalar, T, Int4, T>(name)}) {
+      if (auto callable{GetHostRuntimeWrapper<T, Int4, T>(name)}) {
         return FoldElementalIntrinsic<T, Int4, T>(
             context, std::move(funcRef), *callable);
       } else {
         context.messages().Say(
             "%s(integer(kind=4), real(kind=%d)) cannot be folded on host"_en_US,
             name, KIND);
       }
     }
   } else if (name == "abs") {
     // Argument can be complex or real
     if (auto *x{UnwrapExpr<Expr<SomeReal>>(args[0])}) {
       return FoldElementalIntrinsic<T, T>(
           context, std::move(funcRef), &Scalar<T>::ABS);
     } else if (auto *z{UnwrapExpr<Expr<SomeComplex>>(args[0])}) {
-      if (auto callable{
-              context.hostIntrinsicsLibrary()
-                  .GetHostProcedureWrapper<Scalar, T, ComplexT>("abs")}) {
+      if (auto callable{GetHostRuntimeWrapper<T, ComplexT>("abs")}) {
         return FoldElementalIntrinsic<T, ComplexT>(
             context, std::move(funcRef), *callable);
       } else {
         context.messages().Say(
             "abs(complex(kind=%d)) cannot be folded on host"_en_US, KIND);
       }
     } else {
       common::die(" unexpected argument type inside abs");

flang/lib/Evaluate/host.h

if constexpr (FTN_T::category == TypeCategory::Complex &&
     // complex
     return HostType<FTN_T>{CastFortranToHost<typename FTN_T::Part>(x.REAL()),
         CastFortranToHost<typename FTN_T::Part>(x.AIMAG())};
   } else {
     return *reinterpret_cast<const HostType<FTN_T> *>(&x);
   }
 }
 
-template <typename T> struct BiggerOrSameHostTypeHelper {
-  using Type =
-      std::conditional_t<HostTypeExists<T>(), HostType<T>, UnsupportedType>;
-  using FortranType = T;
-};
-
-template <typename FTN_T>
-using BiggerOrSameHostType = typename BiggerOrSameHostTypeHelper<FTN_T>::Type;
-template <typename FTN_T>
-using BiggerOrSameFortranTypeSupportedOnHost =
-    typename BiggerOrSameHostTypeHelper<FTN_T>::FortranType;
-
-template <typename... T> constexpr inline bool BiggerOrSameHostTypeExists() {
-  return (... && (!std::is_same_v<BiggerOrSameHostType<T>, UnsupportedType>));
-}
-
 // Defining the actual mapping
 template <> struct HostTypeHelper<Type<TypeCategory::Integer, 1>> {
   using Type = std::int8_t;
 };
 
 template <> struct HostTypeHelper<Type<TypeCategory::Integer, 2>> {
   using Type = std::int16_t;
 };

 #else
   using Type = UnsupportedType;
 #endif
 };
 
 // TODO no mapping to host types are defined currently for 16bits float
 // It should be defined when gcc/clang have a better support for it.
 
-template <> struct HostTypeHelper<Type<TypeCategory::Real, 4>> {
-  // IEEE 754 64bits
+template <>
+struct HostTypeHelper<
+    Type<TypeCategory::Real, common::RealKindForPrecision(24)>> {
+  // IEEE 754 32bits
   using Type = std::conditional_t<sizeof(float) == 4 &&
           std::numeric_limits<float>::is_iec559,
       float, UnsupportedType>;
 };
 
-template <> struct HostTypeHelper<Type<TypeCategory::Real, 8>> {
+template <>
+struct HostTypeHelper<
+    Type<TypeCategory::Real, common::RealKindForPrecision(53)>> {
   // IEEE 754 64bits
   using Type = std::conditional_t<sizeof(double) == 8 &&
           std::numeric_limits<double>::is_iec559,
       double, UnsupportedType>;
 };
 
-template <> struct HostTypeHelper<Type<TypeCategory::Real, 10>> {
+template <>
+struct HostTypeHelper<
+    Type<TypeCategory::Real, common::RealKindForPrecision(64)>> {
   // X87 80bits
   using Type = std::conditional_t<sizeof(long double) >= 10 &&
           std::numeric_limits<long double>::digits == 64 &&
           std::numeric_limits<long double>::max_exponent == 16384,
       long double, UnsupportedType>;
 };
 
-template <> struct HostTypeHelper<Type<TypeCategory::Real, 16>> {
+template <>
+struct HostTypeHelper<
+    Type<TypeCategory::Real, common::RealKindForPrecision(113)>> {
   // IEEE 754 128bits
   using Type = std::conditional_t<sizeof(long double) == 16 &&
           std::numeric_limits<long double>::digits == 113 &&
           std::numeric_limits<long double>::max_exponent == 16384,
       long double, UnsupportedType>;
 };
 
 template <int KIND> struct HostTypeHelper<Type<TypeCategory::Complex, KIND>> {

 
 template <typename HOST_T>
 using FortranType = typename FortranTypeHelper<HOST_T>::Type;
 
 template <typename... HT> constexpr inline bool FortranTypeExists() {
   return (... && (!std::is_same_v<FortranType<HT>, UnknownType>));
 }
 
-// Utility to find "bigger" types that exist on host. By bigger, it is meant
-// that the bigger type can represent all the values of the smaller types
-// without information loss.
-template <TypeCategory cat, int KIND> struct NextBiggerReal {
-  using Type = void;
-};
-template <TypeCategory cat> struct NextBiggerReal<cat, 2> {
-  using Type = Fortran::evaluate::Type<cat, 4>;
-};
-template <TypeCategory cat> struct NextBiggerReal<cat, 3> {
-  using Type = Fortran::evaluate::Type<cat, 4>;
-};
-template <TypeCategory cat> struct NextBiggerReal<cat, 4> {
-  using Type = Fortran::evaluate::Type<cat, 8>;
-};
-
-template <TypeCategory cat> struct NextBiggerReal<cat, 8> {
-  using Type = Fortran::evaluate::Type<cat, 10>;
-};
-
-template <TypeCategory cat> struct NextBiggerReal<cat, 10> {
-  using Type = Fortran::evaluate::Type<cat, 16>;
-};
-
-template <int KIND>
-struct BiggerOrSameHostTypeHelper<Type<TypeCategory::Real, KIND>> {
-  using T = Fortran::evaluate::Type<TypeCategory::Real, KIND>;
-  using NextT = typename NextBiggerReal<TypeCategory::Real, KIND>::Type;
-  using Type = std::conditional_t<HostTypeExists<T>(), HostType<T>,
-      typename BiggerOrSameHostTypeHelper<NextT>::Type>;
-  using FortranType = std::conditional_t<HostTypeExists<T>(), T,
-      typename BiggerOrSameHostTypeHelper<NextT>::FortranType>;
-};
-
-template <int KIND>
-struct BiggerOrSameHostTypeHelper<Type<TypeCategory::Complex, KIND>> {
-  using T = Fortran::evaluate::Type<TypeCategory::Complex, KIND>;
-  using NextT = typename NextBiggerReal<TypeCategory::Complex, KIND>::Type;
-  using Type = std::conditional_t<HostTypeExists<T>(), HostType<T>,
-      typename BiggerOrSameHostTypeHelper<NextT>::Type>;
-  using FortranType = std::conditional_t<HostTypeExists<T>(), T,
-      typename BiggerOrSameHostTypeHelper<NextT>::FortranType>;
-};
 } // namespace host
 } // namespace Fortran::evaluate
 
 #endif // FORTRAN_EVALUATE_HOST_H_

flang/lib/Evaluate/intrinsics-library-templates.h

-//===-- lib/Evaluate/intrinsics-library-templates.h -------------*- 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
-//
-//===----------------------------------------------------------------------===//
-
-#ifndef FORTRAN_EVALUATE_INTRINSICS_LIBRARY_TEMPLATES_H_
-#define FORTRAN_EVALUATE_INTRINSICS_LIBRARY_TEMPLATES_H_
-
-// This header defines the actual implementation of the templatized member
-// function of the structures defined in intrinsics-library.h. It should only be
-// included if these member functions are used, else intrinsics-library.h is
-// sufficient. This is to avoid circular dependencies. The below implementation
-// cannot be defined in .cpp file because it would be too cumbersome to decide
-// which version should be instantiated in a generic way.
-
-#include "host.h"
-#include "flang/Common/template.h"
-#include "flang/Evaluate/intrinsics-library.h"
-#include "flang/Evaluate/type.h"
-
-#include <tuple>
-#include <type_traits>
-
-namespace Fortran::evaluate {
-
-// Define meaningful types for the runtime
-using RuntimeTypes = evaluate::AllIntrinsicTypes;
-
-template <typename T, typename... TT> struct IndexInTupleHelper {};
-template <typename T, typename... TT>
-struct IndexInTupleHelper<T, std::tuple<TT...>> {
-  static constexpr TypeCode value{common::TypeIndex<T, TT...>};
-};
-
-static_assert(
-    std::tuple_size_v<RuntimeTypes> < std::numeric_limits<TypeCode>::max(),
-    "TypeCode is too small");
-template <typename T>
-inline constexpr TypeCode typeCodeOf{
-    IndexInTupleHelper<T, RuntimeTypes>::value};
-
-template <TypeCode n>
-using RuntimeTypeOf = typename std::tuple_element_t<n, RuntimeTypes>;
-
-template <typename TA, PassBy Pass>
-using HostArgType = std::conditional_t<Pass == PassBy::Ref,
-    std::add_lvalue_reference_t<std::add_const_t<host::HostType<TA>>>,
-    host::HostType<TA>>;
-
-template <typename TR, typename... ArgInfo>
-using HostFuncPointer = FuncPointer<host::HostType<TR>,
-    HostArgType<typename ArgInfo::Type, ArgInfo::pass>...>;
-
-// Software Subnormal Flushing helper.
-template <typename T> struct Flusher {
-  // Only flush floating-points. Forward other scalars untouched.
-  static constexpr inline const Scalar<T> &FlushSubnormals(const Scalar<T> &x) {
-    return x;
-  }
-};
-template <int Kind> struct Flusher<Type<TypeCategory::Real, Kind>> {
-  using T = Type<TypeCategory::Real, Kind>;
-  static constexpr inline Scalar<T> FlushSubnormals(const Scalar<T> &x) {
-    return x.FlushSubnormalToZero();
-  }
-};
-template <int Kind> struct Flusher<Type<TypeCategory::Complex, Kind>> {
-  using T = Type<TypeCategory::Complex, Kind>;
-  static constexpr inline Scalar<T> FlushSubnormals(const Scalar<T> &x) {
-    return x.FlushSubnormalToZero();
-  }
-};
-
-// Callable factory
-template <typename TR, typename... ArgInfo> struct CallableHostWrapper {
-  static Scalar<TR> scalarCallable(FoldingContext &context,
-      HostFuncPointer<TR, ArgInfo...> func,
-      const Scalar<typename ArgInfo::Type> &...x) {
-    if constexpr (host::HostTypeExists<TR, typename ArgInfo::Type...>()) {
-      host::HostFloatingPointEnvironment hostFPE;
-      hostFPE.SetUpHostFloatingPointEnvironment(context);
-      host::HostType<TR> hostResult{};
-      Scalar<TR> result{};
-      if (context.flushSubnormalsToZero() &&
-          !hostFPE.hasSubnormalFlushingHardwareControl()) {
-        hostResult = func(host::CastFortranToHost<typename ArgInfo::Type>(
-            Flusher<typename ArgInfo::Type>::FlushSubnormals(x))...);
-        result = Flusher<TR>::FlushSubnormals(
-            host::CastHostToFortran<TR>(hostResult));
-      } else {
-        hostResult =
-            func(host::CastFortranToHost<typename ArgInfo::Type>(x)...);
-        result = host::CastHostToFortran<TR>(hostResult);
-      }
-      if (!hostFPE.hardwareFlagsAreReliable()) {
-        CheckFloatingPointIssues(hostFPE, result);
-      }
-      hostFPE.CheckAndRestoreFloatingPointEnvironment(context);
-      return result;
-    } else {
-      common::die("Internal error: Host does not supports this function type."
-                  "This should not have been called for folding");
-    }
-  }
-  static constexpr inline auto MakeScalarCallable() { return &scalarCallable; }
-
-  static void CheckFloatingPointIssues(
-      host::HostFloatingPointEnvironment &hostFPE, const Scalar<TR> &x) {
-    if constexpr (TR::category == TypeCategory::Complex ||
-        TR::category == TypeCategory::Real) {
-      if (x.IsNotANumber()) {
-        hostFPE.SetFlag(RealFlag::InvalidArgument);
-      } else if (x.IsInfinite()) {
-        hostFPE.SetFlag(RealFlag::Overflow);
-      }
-    }
-  }
-};
-
-template <typename TR, typename... TA>
-inline GenericFunctionPointer ToGenericFunctionPointer(
-    FuncPointer<TR, TA...> f) {
-  return reinterpret_cast<GenericFunctionPointer>(f);
-}
-
-template <typename TR, typename... TA>
-inline FuncPointer<TR, TA...> FromGenericFunctionPointer(
-    GenericFunctionPointer g) {
-  return reinterpret_cast<FuncPointer<TR, TA...>>(g);
-}
-
-template <typename TR, typename... ArgInfo>
-IntrinsicProcedureRuntimeDescription::IntrinsicProcedureRuntimeDescription(
-    const Signature<TR, ArgInfo...> &signature, bool isElemental)
-    : name{signature.name}, returnType{typeCodeOf<TR>},
-      argumentsType{typeCodeOf<typename ArgInfo::Type>...},
-      argumentsPassedBy{ArgInfo::pass...}, isElemental{isElemental},
-      callable{ToGenericFunctionPointer(
-          CallableHostWrapper<TR, ArgInfo...>::MakeScalarCallable())} {}
-
-template <typename HostTA> static constexpr inline PassBy PassByMethod() {
-  if constexpr (std::is_pointer_v<std::decay_t<HostTA>> ||
-      std::is_lvalue_reference_v<HostTA>) {
-    return PassBy::Ref;
-  }
-  return PassBy::Val;
-}
-
-template <typename HostTA>
-using ArgInfoFromHostType =
-    ArgumentInfo<host::FortranType<std::remove_pointer_t<std::decay_t<HostTA>>>,
-        PassByMethod<HostTA>()>;
-
-template <typename HostTR, typename... HostTA>
-using SignatureFromHostFuncPointer =
-    Signature<host::FortranType<HostTR>, ArgInfoFromHostType<HostTA>...>;
-
-template <typename HostTR, typename... HostTA>
-HostRuntimeIntrinsicProcedure::HostRuntimeIntrinsicProcedure(
-    const std::string &name, FuncPointer<HostTR, HostTA...> func,
-    bool isElemental)
-    : IntrinsicProcedureRuntimeDescription(
-          SignatureFromHostFuncPointer<HostTR, HostTA...>{name}, isElemental),
-      handle{ToGenericFunctionPointer(func)} {}
-
-template <template <typename> typename ConstantContainer, typename TR,
-    typename... TA>
-std::optional<HostProcedureWrapper<ConstantContainer, TR, TA...>>
-HostIntrinsicProceduresLibrary::GetHostProcedureWrapper(
-    const std::string &name) const {
-  if constexpr (host::HostTypeExists<TR, TA...>()) {
-    auto rteProcRange{procedures_.equal_range(name)};
-    const TypeCode resTypeCode{typeCodeOf<TR>};
-    const std::vector<TypeCode> argTypes{typeCodeOf<TA>...};
-    const size_t nargs{argTypes.size()};
-    for (auto iter{rteProcRange.first}; iter != rteProcRange.second; ++iter) {
-      if (nargs == iter->second.argumentsType.size() &&
-          resTypeCode == iter->second.returnType &&
-          (!std::is_same_v<ConstantContainer<TR>, Scalar<TR>> ||
-              iter->second.isElemental)) {
-        bool match{true};
-        int pos{0};
-        for (auto const &type : argTypes) {
-          if (type != iter->second.argumentsType[pos++]) {
-            match = false;
-            break;
-          }
-        }
-        if (match) {
-          return {HostProcedureWrapper<ConstantContainer, TR, TA...>{
-              [=](FoldingContext &context,
-                  const ConstantContainer<TA> &...args) {
-                auto callable{FromGenericFunctionPointer<ConstantContainer<TR>,
-                    FoldingContext &, GenericFunctionPointer,
-                    const ConstantContainer<TA> &...>(iter->second.callable)};
-                return callable(context, iter->second.handle, args...);
-              }}};
-        }
-      }
-    }
-  }
-  return std::nullopt;
-}
-
-} // namespace Fortran::evaluate
-#endif // FORTRAN_EVALUATE_INTRINSICS_LIBRARY_TEMPLATES_H_

flang/lib/Evaluate/intrinsics-library.cpp

 //===-- lib/Evaluate/intrinsics-library.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
 //
 //===----------------------------------------------------------------------===//
 
 // This file defines host runtime functions that can be used for folding
 // intrinsic functions.
-// The default HostIntrinsicProceduresLibrary is built with <cmath> and
+// The default host runtime folders are built with <cmath> and
 // <complex> functions that are guaranteed to exist from the C++ standard.
 
-#include "intrinsics-library-templates.h"
+#include "flang/Evaluate/intrinsics-library.h"
+#include "fold-implementation.h"
+#include "host.h"
+#include "flang/Common/static-multimap-view.h"
+#include "flang/Evaluate/expression.h"
 #include <cmath>
 #include <complex>
+#include <functional>
+#include <type_traits>
 
 namespace Fortran::evaluate {
 
-// Note: argument passing is ignored in equivalence
-bool HostIntrinsicProceduresLibrary::HasEquivalentProcedure(
-    const IntrinsicProcedureRuntimeDescription &sym) const {
-  const auto rteProcRange{procedures_.equal_range(sym.name)};
-  const size_t nargs{sym.argumentsType.size()};
-  for (auto iter{rteProcRange.first}; iter != rteProcRange.second; ++iter) {
-    if (nargs == iter->second.argumentsType.size() &&
-        sym.returnType == iter->second.returnType &&
-        (sym.isElemental || iter->second.isElemental)) {
-      bool match{true};
-      int pos{0};
-      for (const auto &type : sym.argumentsType) {
-        if (type != iter->second.argumentsType[pos++]) {
-          match = false;
-          break;
-        }
-      }
-      if (match) {
-        return true;
-      }
-    }
+// Define a vector like class that can hold an arbitrary number of
+// Dynamic type and be built at compile time. This is like a
+// std::vector<DynamicType>, but constexpr only.
+template <typename... FortranType> struct TypeVectorStorage {
+  static constexpr DynamicType values[]{FortranType{}.GetType()...};
+  static constexpr const DynamicType *start{&values[0]};
+  static constexpr const DynamicType *end{start + sizeof...(FortranType)};
+};
+template <> struct TypeVectorStorage<> {
+  static constexpr const DynamicType *start{nullptr}, *end{nullptr};
+};
+struct TypeVector {
+  template <typename... FortranType> static constexpr TypeVector Create() {
+    using storage = TypeVectorStorage<FortranType...>;
+    return TypeVector{storage::start, storage::end, sizeof...(FortranType)};
+  }
+  constexpr size_t size() const { return size_; };
+  using const_iterator = const DynamicType *;
+  constexpr const_iterator begin() const { return startPtr; }
+  constexpr const_iterator end() const { return endPtr; }
+  const DynamicType &operator[](size_t i) const { return *(startPtr + i); }
+
+  const DynamicType *startPtr{nullptr};
+  const DynamicType *endPtr{nullptr};
+  const size_t size_;
+};
+inline bool operator==(
+    const TypeVector &lhs, const std::vector<DynamicType> &rhs) {
+  if (lhs.size() != rhs.size()) {
+    return false;
+  }
+  for (size_t i{0}; i < lhs.size(); ++i) {
+    if (lhs[i] != rhs[i]) {
+      return false;
+
+// HostRuntimeFunction holds a pointer to a Folder function that can fold
+// a Fortran scalar intrinsic using host runtime functions (e.g libm).
+// The folder take care of all conversions between Fortran types and the related
+// host types as well as setting and cleaning-up the floating point environment.
+// HostRuntimeFunction are intended to be built at compile time (members are all
+// constexpr constructible) so that they can be stored in a compile time static
+// map.
+struct HostRuntimeFunction {
+  using Folder = Expr<SomeType> (*)(
+      FoldingContext &, std::vector<Expr<SomeType>> &&);
+  using Key = std::string_view;
+  // Needed for implicit compare with keys.
+  constexpr operator Key() const { return key; }
+  // Name of the related Fortran intrinsic.
+  Key key;
+  // DynamicType of the Expr<SomeType> returns by folder.
+  DynamicType resultType;
+  // DynamicTypes expected for the Expr<SomeType> arguments of the folder.
+  // The folder will crash if provided arguments of different types.
+  TypeVector argumentTypes;
+  // Folder to be called to fold the intrinsic with host runtime. The provided
+  // Expr<SomeType> arguments must wrap scalar constants of the type described
+  // in argumentTypes, otherwise folder will crash. Any floating point issue
+  // raised while executing the host runtime will be reported in FoldingContext
+  // messages.
+  Folder folder;
+};
+
+// Translate a host function type signature (template arguments) into a
+// constexpr data representation based on Fortran DynamicType that can be
+// stored.
+template <typename TR, typename... TA> using FuncPointer = TR (*)(TA...);
+template <typename T> struct FuncTypeAnalyzer {};
+template <typename HostTR, typename... HostTA>
+struct FuncTypeAnalyzer<FuncPointer<HostTR, HostTA...>> {
+  static constexpr DynamicType result{host::FortranType<HostTR>{}.GetType()};
+  static constexpr TypeVector arguments{
+      TypeVector::Create<host::FortranType<HostTA>...>()};
+};
+
+// Define helpers to deal with host floating environment.
+template <typename TR>
+static void CheckFloatingPointIssues(
+    host::HostFloatingPointEnvironment &hostFPE, const Scalar<TR> &x) {
+  if constexpr (TR::category == TypeCategory::Complex ||
+      TR::category == TypeCategory::Real) {
+    if (x.IsNotANumber()) {
+      hostFPE.SetFlag(RealFlag::InvalidArgument);
+    } else if (x.IsInfinite()) {
+      hostFPE.SetFlag(RealFlag::Overflow);
+    }
+  }
+}
+// Software Subnormal Flushing helper.
+// Only flush floating-points. Forward other scalars untouched.
+// Software flushing is only performed if hardware flushing is not available
+// because it may not result in the same behavior as hardware flushing.
+// Some runtime implementations are "working around" subnormal flushing to
+// return results that they deem better than returning the result they would
+// with a null argument. An example is logf that should return -inf if arguments
+// are flushed to zero, but some implementations return -1.03972076416015625e2_4
+// for all subnormal values instead. It is impossible to reproduce this with the
+// simple software flushing below.
+template <typename T>
+static constexpr inline const Scalar<T> FlushSubnormals(Scalar<T> &&x) {
+  if constexpr (T::category == TypeCategory::Real ||
+      T::category == TypeCategory::Complex) {
+    return x.FlushSubnormalToZero();
+  }
+  return x;
+}
+
+// This is the kernel called by all HostRuntimeFunction folders, it convert the
+// Fortran Expr<SomeType> to the host runtime function argument types, calls
+// the runtime function, and wrap back the result into an Expr<SomeType>.
+// It deals with host floating point environment set-up and clean-up.
+template <typename FuncType, typename TR, typename... TA, size_t... I>
+static Expr<SomeType> ApplyHostFunctionHelper(FuncType func,
+    FoldingContext &context, std::vector<Expr<SomeType>> &&args,
+    std::index_sequence<I...>) {
+  host::HostFloatingPointEnvironment hostFPE;
+  hostFPE.SetUpHostFloatingPointEnvironment(context);
+  host::HostType<TR> hostResult{};
+  Scalar<TR> result{};
+  std::tuple<Scalar<TA>...> scalarArgs{
+      GetScalarConstantValue<TA>(args[I]).value()...};
+  if (context.flushSubnormalsToZero() &&
+      !hostFPE.hasSubnormalFlushingHardwareControl()) {
+    hostResult = func(host::CastFortranToHost<TA>(
+        FlushSubnormals<TA>(std::move(std::get<I>(scalarArgs))))...);
+    result = FlushSubnormals<TR>(host::CastHostToFortran<TR>(hostResult));
+  } else {
+    hostResult = func(host::CastFortranToHost<TA>(std::get<I>(scalarArgs))...);
+    result = host::CastHostToFortran<TR>(hostResult);
+  if (!hostFPE.hardwareFlagsAreReliable()) {
+    CheckFloatingPointIssues<TR>(hostFPE, result);
   }
-  return false;
+  hostFPE.CheckAndRestoreFloatingPointEnvironment(context);
+  return AsGenericExpr(Constant<TR>(std::move(result)));
+}
+template <typename HostTR, typename... HostTA>
+Expr<SomeType> ApplyHostFunction(FuncPointer<HostTR, HostTA...> func,
+    FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
+  return ApplyHostFunctionHelper<decltype(func), host::FortranType<HostTR>,
+      host::FortranType<HostTA>...>(
+      func, context, std::move(args), std::index_sequence_for<HostTA...>{});
 }
 
-// Map numerical intrinsic to  <cmath>/<complex> functions
+// FolderFactory builds a HostRuntimeFunction for the host runtime function
+// passed as a template argument.
+// Its static member function "fold" is the resulting folder. It captures the
+// host runtime function pointer and pass it to the host runtime function folder
+// kernel.
+template <typename HostFuncType, HostFuncType func> class FolderFactory {
+public:
+  static constexpr HostRuntimeFunction Create(const std::string_view &name) {
+    return HostRuntimeFunction{name, FuncTypeAnalyzer<HostFuncType>::result,
+        FuncTypeAnalyzer<HostFuncType>::arguments, &Fold};
+  }
+
+private:
+  static Expr<SomeType> Fold(
+      FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
+    return ApplyHostFunction(func, context, std::move(args));
+  }
+};
 
-// Define which host runtime functions will be used for folding
+// Define host runtime libraries that can be used for folding and
+// fill their description if they are available.
+enum class LibraryVersion { Libm, PgmathFast, PgmathRelaxed, PgmathPrecise };
+template <typename HostT, LibraryVersion> struct HostRuntimeLibrary {
+  // When specialized, this class holds a static constexpr table containing
+  // all the HostRuntimeLibrary for functions of library LibraryVersion
+  // that returns a value of type HostT.
+};
+
+using HostRuntimeMap = common::StaticMultimapView<HostRuntimeFunction>;
 
+// Map numerical intrinsic to  <cmath>/<complex> functions
 template <typename HostT>
-static void AddLibmRealHostProcedures(
-    HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
+struct HostRuntimeLibrary<HostT, LibraryVersion::Libm> {
   using F = FuncPointer<HostT, HostT>;
   using F2 = FuncPointer<HostT, HostT, HostT>;
-  HostRuntimeIntrinsicProcedure libmSymbols[]{
-      {"acos", F{std::acos}, true},
-      {"acosh", F{std::acosh}, true},
-      {"asin", F{std::asin}, true},
-      {"asinh", F{std::asinh}, true},
-      {"atan", F{std::atan}, true},
-      {"atan2", F2{std::atan2}, true},
-      {"atanh", F{std::atanh}, true},
-      {"cos", F{std::cos}, true},
-      {"cosh", F{std::cosh}, true},
-      {"erf", F{std::erf}, true},
-      {"erfc", F{std::erfc}, true},
-      {"exp", F{std::exp}, true},
-      {"gamma", F{std::tgamma}, true},
-      {"hypot", F2{std::hypot}, true},
-      {"log", F{std::log}, true},
-      {"log10", F{std::log10}, true},
-      {"log_gamma", F{std::lgamma}, true},
-      {"mod", F2{std::fmod}, true},
-      {"pow", F2{std::pow}, true},
-      {"sin", F{std::sin}, true},
-      {"sinh", F{std::sinh}, true},
-      {"sqrt", F{std::sqrt}, true},
-      {"tan", F{std::tan}, true},
-      {"tanh", F{std::tanh}, true},
+  using ComplexToRealF = FuncPointer<HostT, const std::complex<HostT> &>;
+  static constexpr HostRuntimeFunction table[]{
+      FolderFactory<ComplexToRealF, ComplexToRealF{std::abs}>::Create("abs"),
+      FolderFactory<F, F{std::acos}>::Create("acos"),
+      FolderFactory<F, F{std::acosh}>::Create("acosh"),
+      FolderFactory<F, F{std::asin}>::Create("asin"),
+      FolderFactory<F, F{std::asinh}>::Create("asinh"),
+      FolderFactory<F, F{std::atan}>::Create("atan"),
+      FolderFactory<F2, F2{std::atan2}>::Create("atan2"),
+      FolderFactory<F, F{std::atanh}>::Create("atanh"),
+      FolderFactory<F, F{std::cos}>::Create("cos"),
+      FolderFactory<F, F{std::cosh}>::Create("cosh"),
+      FolderFactory<F, F{std::erf}>::Create("erf"),
+      FolderFactory<F, F{std::erfc}>::Create("erfc"),
+      FolderFactory<F, F{std::exp}>::Create("exp"),
+      FolderFactory<F, F{std::tgamma}>::Create("gamma"),
+      FolderFactory<F2, F2{std::hypot}>::Create("hypot"),
+      FolderFactory<F, F{std::log}>::Create("log"),
+      FolderFactory<F, F{std::log10}>::Create("log10"),
+      FolderFactory<F, F{std::lgamma}>::Create("log_gamma"),
+      FolderFactory<F2, F2{std::fmod}>::Create("mod"),
+      FolderFactory<F2, F2{std::pow}>::Create("pow"),
+      FolderFactory<F, F{std::sin}>::Create("sin"),
+      FolderFactory<F, F{std::sinh}>::Create("sinh"),
+      FolderFactory<F, F{std::sqrt}>::Create("sqrt"),
+      FolderFactory<F, F{std::tan}>::Create("tan"),
+      FolderFactory<F, F{std::tanh}>::Create("tanh"),
   };
   // Note: cmath does not have modulo and erfc_scaled equivalent
 
   // Note regarding  lack of bessel function support:
   // C++17 defined standard Bessel math functions std::cyl_bessel_j
   // and std::cyl_neumann that can be used for Fortran j and y
   // bessel functions. However, they are not yet implemented in
   // clang libc++ (ok in GNU libstdc++). C maths functions j0...
   // are not C standard but a GNU extension so they are not used
   // to avoid introducing incompatibilities.
   // Use libpgmath to get bessel function folding support.
   // TODO:  Add Bessel functions when possible.
-
-  for (auto sym : libmSymbols) {
-    if (!hostIntrinsicLibrary.HasEquivalentProcedure(sym)) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  }
-}
-
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
 template <typename HostT>
-static void AddLibmComplexHostProcedures(
-    HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
+struct HostRuntimeLibrary<std::complex<HostT>, LibraryVersion::Libm> {
   using F = FuncPointer<std::complex<HostT>, const std::complex<HostT> &>;
   using F2 = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
       const std::complex<HostT> &>;
-  using F2a = FuncPointer<std::complex<HostT>, const HostT &,
+  using F2A = FuncPointer<std::complex<HostT>, const HostT &,
       const std::complex<HostT> &>;
-  using F2b = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
+  using F2B = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
       const HostT &>;
-  HostRuntimeIntrinsicProcedure libmSymbols[]{
-      {"abs", FuncPointer<HostT, const std::complex<HostT> &>{std::abs}, true},
-      {"acos", F{std::acos}, true},
-      {"acosh", F{std::acosh}, true},
-      {"asin", F{std::asin}, true},
-      {"asinh", F{std::asinh}, true},
-      {"atan", F{std::atan}, true},
-      {"atanh", F{std::atanh}, true},
-      {"cos", F{std::cos}, true},
-      {"cosh", F{std::cosh}, true},
-      {"exp", F{std::exp}, true},
-      {"log", F{std::log}, true},
-      {"pow", F2{std::pow}, true},
-      {"pow", F2a{std::pow}, true},
-      {"pow", F2b{std::pow}, true},
-      {"sin", F{std::sin}, true},
-      {"sinh", F{std::sinh}, true},
-      {"sqrt", F{std::sqrt}, true},
-      {"tan", F{std::tan}, true},
-      {"tanh", F{std::tanh}, true},
-  };
-
-  for (auto sym : libmSymbols) {
-    if (!hostIntrinsicLibrary.HasEquivalentProcedure(sym)) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  }
-}
-
-[[maybe_unused]] static void InitHostIntrinsicLibraryWithLibm(
-    HostIntrinsicProceduresLibrary &lib) {
-  if constexpr (host::FortranTypeExists<float>()) {
-    AddLibmRealHostProcedures<float>(lib);
-  }
-  if constexpr (host::FortranTypeExists<double>()) {
-    AddLibmRealHostProcedures<double>(lib);
-  }
-  if constexpr (host::FortranTypeExists<long double>()) {
-    AddLibmRealHostProcedures<long double>(lib);
-  }
-
-  if constexpr (host::FortranTypeExists<std::complex<float>>()) {
-    AddLibmComplexHostProcedures<float>(lib);
-  }
-  if constexpr (host::FortranTypeExists<std::complex<double>>()) {
-    AddLibmComplexHostProcedures<double>(lib);
-  }
-  if constexpr (host::FortranTypeExists<std::complex<long double>>()) {
-    AddLibmComplexHostProcedures<long double>(lib);
-  }
-}
+  static constexpr HostRuntimeFunction table[]{
+      FolderFactory<F, F{std::acos}>::Create("acos"),
+      FolderFactory<F, F{std::acosh}>::Create("acosh"),
+      FolderFactory<F, F{std::asin}>::Create("asin"),
+      FolderFactory<F, F{std::asinh}>::Create("asinh"),
+      FolderFactory<F, F{std::atan}>::Create("atan"),
+      FolderFactory<F, F{std::atanh}>::Create("atanh"),
+      FolderFactory<F, F{std::cos}>::Create("cos"),
+      FolderFactory<F, F{std::cosh}>::Create("cosh"),
+      FolderFactory<F, F{std::exp}>::Create("exp"),
+      FolderFactory<F, F{std::log}>::Create("log"),
+      FolderFactory<F2, F2{std::pow}>::Create("pow"),
+      FolderFactory<F2A, F2A{std::pow}>::Create("pow"),
+      FolderFactory<F2B, F2B{std::pow}>::Create("pow"),
+      FolderFactory<F, F{std::sin}>::Create("sin"),
+      FolderFactory<F, F{std::sinh}>::Create("sinh"),
+      FolderFactory<F, F{std::sqrt}>::Create("sqrt"),
+      FolderFactory<F, F{std::tan}>::Create("tan"),
+      FolderFactory<F, F{std::tanh}>::Create("tanh"),
+  };
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
 
+/// Define pgmath description
 #if LINK_WITH_LIBPGMATH
 // Only use libpgmath for folding if it is available.
 // First declare all libpgmaths functions
+#define PGMATH_LINKING
 #define PGMATH_DECLARE
 #include "../runtime/pgmath.h.inc"
 
-// Library versions: P for Precise, R for Relaxed, F for Fast
-enum class L { F, R, P };
-
-// Fill the function map used for folding with libpgmath symbols
-template <L Lib>
-static void AddLibpgmathFloatHostProcedures(
-    HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
-  if constexpr (Lib == L::F) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
+#define REAL_FOLDER(name, func) \
+  FolderFactory<decltype(&func), &func>::Create(#name)
+template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathFast> {
+  static constexpr HostRuntimeFunction table[]{
 #define PGMATH_FAST
-#define PGMATH_USE_S(name, function) {#name, function, true},
+#define PGMATH_USE_S(name, func) REAL_FOLDER(name, func),
 #include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  } else if constexpr (Lib == L::R) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_RELAXED
-#define PGMATH_USE_S(name, function) {#name, function, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  } else {
-    static_assert(Lib == L::P && "unexpected libpgmath version");
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_PRECISE
-#define PGMATH_USE_S(name, function) {#name, function, true},
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathFast> {
+  static constexpr HostRuntimeFunction table[]{
+#define PGMATH_FAST
+#define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
 #include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  }
-}
-
-template <L Lib>
-static void AddLibpgmathDoubleHostProcedures(
-    HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
-  if constexpr (Lib == L::F) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_FAST
-#define PGMATH_USE_D(name, function) {#name, function, true},
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
+template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathRelaxed> {
+  static constexpr HostRuntimeFunction table[]{
+#define PGMATH_RELAXED
+#define PGMATH_USE_S(name, func) REAL_FOLDER(name, func),
 #include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  } else if constexpr (Lib == L::R) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
+template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathRelaxed> {
+  static constexpr HostRuntimeFunction table[]{
 #define PGMATH_RELAXED
-#define PGMATH_USE_D(name, function) {#name, function, true},
+#define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
 #include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  } else {
-    static_assert(Lib == L::P && "unexpected libpgmath version");
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
+template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathPrecise> {
+  static constexpr HostRuntimeFunction table[]{
 #define PGMATH_PRECISE
-#define PGMATH_USE_D(name, function) {#name, function, true},
+#define PGMATH_USE_S(name, func) REAL_FOLDER(name, func),
 #include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  }
-}
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
+template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathPrecise> {
+  static constexpr HostRuntimeFunction table[]{
+#define PGMATH_PRECISE
+#define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
+#include "../runtime/pgmath.h.inc"
+  };
+  static constexpr HostRuntimeMap map{table};
+  static_assert(map.Verify(), "map must be sorted");
+};
 
-// Note: Lipgmath uses _Complex but the front-end use std::complex for folding.
-// std::complex and _Complex are layout compatible but are not guaranteed
-// to be linkage compatible. For instance, on i386, float _Complex is returned
-// by a pair of register but std::complex<float> is returned by structure
-// address. To fix the issue, wrapper around C _Complex functions are defined
-// below.
+// TODO: double _Complex/float _Complex have been removed from llvm flang
+// pgmath.h.inc because they caused warnings, they need to be added back
+// so that the complex pgmath versions can be used when requested.
 
-template <typename T> struct ToStdComplex {
-  using Type = T;
-  using AType = Type;
-};
+#endif /* LINK_WITH_LIBPGMATH */
 
-template <typename F, F func> struct CComplexFunc {};
-template <typename R, typename... A, FuncPointer<R, A...> func>
-struct CComplexFunc<FuncPointer<R, A...>, func> {
-  static typename ToStdComplex<R>::Type wrapper(
-      typename ToStdComplex<A>::AType... args) {
-    R res{func(*reinterpret_cast<A *>(&args)...)};
-    return *reinterpret_cast<typename ToStdComplex<R>::Type *>(&res);
+// Helper to check if a HostRuntimeLibrary specialization exists
+template <typename T, typename = void> struct IsAvailable : std::false_type {};
+template <typename T>
+struct IsAvailable<T, decltype((void)T::table, void())> : std::true_type {};
+// Define helpers to find host runtime library map according to desired version
+// and type.
+template <typename HostT, LibraryVersion version>
+static const HostRuntimeMap *GetHostRuntimeMapHelper(
+    [[maybe_unused]] DynamicType resultType) {
+  // A library must only be instantiated if LibraryVersion is
+  // available on the host and if HostT maps to a Fortran type.
+  // For instance, whenever long double and double are both 64-bits, double
+  // is mapped to Fortran 64bits real type, and long double will be left
+  // unmapped.
+  if constexpr (host::FortranTypeExists<HostT>()) {
+    using Lib = HostRuntimeLibrary<HostT, version>;
+    if constexpr (IsAvailable<Lib>::value) {
+      if (host::FortranType<HostT>{}.GetType() == resultType) {
+        return &Lib::map;
       }
-};
-
-template <L Lib>
-static void AddLibpgmathComplexHostProcedures(
-    HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
-  if constexpr (Lib == L::F) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_FAST
-#define PGMATH_USE_C(name, function) \
-  {#name, CComplexFunc<decltype(&function), &function>::wrapper, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
     }
-  } else if constexpr (Lib == L::R) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_RELAXED
-#define PGMATH_USE_C(name, function) \
-  {#name, CComplexFunc<decltype(&function), &function>::wrapper, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
-  } else {
-    static_assert(Lib == L::P && "unexpected libpgmath version");
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_PRECISE
-#define PGMATH_USE_C(name, function) \
-  {#name, CComplexFunc<decltype(&function), &function>::wrapper, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
+  return nullptr;
+template <LibraryVersion version>
+static const HostRuntimeMap *GetHostRuntimeMapVersion(DynamicType resultType) {
+  if (resultType.category() == TypeCategory::Real) {
+    if (const auto *map{GetHostRuntimeMapHelper<float, version>(resultType)}) {
+      return map;
     }
-
-  // cmath is used to complement pgmath when symbols are not available
-  using HostT = float;
-  using CHostT = std::complex<HostT>;
-  using CmathF = FuncPointer<CHostT, const CHostT &>;
-  hostIntrinsicLibrary.AddProcedure(
-      {"abs", FuncPointer<HostT, const CHostT &>{std::abs}, true});
-  hostIntrinsicLibrary.AddProcedure({"acosh", CmathF{std::acosh}, true});
-  hostIntrinsicLibrary.AddProcedure({"asinh", CmathF{std::asinh}, true});
-  hostIntrinsicLibrary.AddProcedure({"atanh", CmathF{std::atanh}, true});
-}
-
-template <L Lib>
-static void AddLibpgmathDoubleComplexHostProcedures(
-    HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
-  if constexpr (Lib == L::F) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_FAST
-#define PGMATH_USE_Z(name, function) \
-  {#name, CComplexFunc<decltype(&function), &function>::wrapper, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
+    if (const auto *map{GetHostRuntimeMapHelper<double, version>(resultType)}) {
+      return map;
     }
-  } else if constexpr (Lib == L::R) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_RELAXED
-#define PGMATH_USE_Z(name, function) \
-  {#name, CComplexFunc<decltype(&function), &function>::wrapper, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
+    if (const auto *map{
+            GetHostRuntimeMapHelper<long double, version>(resultType)}) {
+      return map;
     }
-  } else {
-    static_assert(Lib == L::P && "unexpected libpgmath version");
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_PRECISE
-#define PGMATH_USE_Z(name, function) \
-  {#name, CComplexFunc<decltype(&function), &function>::wrapper, true},
-#include "../runtime/pgmath.h.inc"
-    };
-    for (auto sym : pgmathSymbols) {
-      hostIntrinsicLibrary.AddProcedure(std::move(sym));
-    }
+  if (resultType.category() == TypeCategory::Complex) {
+    if (const auto *map{GetHostRuntimeMapHelper<std::complex<float>, version>(
+            resultType)}) {
+      return map;
     }
-
-  // cmath is used to complement pgmath when symbols are not available
-  using HostT = double;
-  using CHostT = std::complex<HostT>;
-  using CmathF = FuncPointer<CHostT, const CHostT &>;
-  hostIntrinsicLibrary.AddProcedure(
-      {"abs", FuncPointer<HostT, const CHostT &>{std::abs}, true});
-  hostIntrinsicLibrary.AddProcedure({"acosh", CmathF{std::acosh}, true});
-  hostIntrinsicLibrary.AddProcedure({"asinh", CmathF{std::asinh}, true});
-  hostIntrinsicLibrary.AddProcedure({"atanh", CmathF{std::atanh}, true});
+    if (const auto *map{GetHostRuntimeMapHelper<std::complex<double>, version>(
+            resultType)}) {
+      return map;
     }
-
-template <L Lib>
-static void InitHostIntrinsicLibraryWithLibpgmath(
-    HostIntrinsicProceduresLibrary &lib) {
-  if constexpr (host::FortranTypeExists<float>()) {
-    AddLibpgmathFloatHostProcedures<Lib>(lib);
+    if (const auto *map{
+            GetHostRuntimeMapHelper<std::complex<long double>, version>(
+                resultType)}) {
+      return map;
     }
-  if constexpr (host::FortranTypeExists<double>()) {
-    AddLibpgmathDoubleHostProcedures<Lib>(lib);
   }
-  if constexpr (host::FortranTypeExists<std::complex<float>>()) {
-    AddLibpgmathComplexHostProcedures<Lib>(lib);
-  }
-  if constexpr (host::FortranTypeExists<std::complex<double>>()) {
-    AddLibpgmathDoubleComplexHostProcedures<Lib>(lib);
+  return nullptr;
+static const HostRuntimeMap *GetHostRuntimeMap(
+    LibraryVersion version, DynamicType resultType) {
+  switch (version) {
+  case LibraryVersion::Libm:
+    return GetHostRuntimeMapVersion<LibraryVersion::Libm>(resultType);
+  case LibraryVersion::PgmathPrecise:
+    return GetHostRuntimeMapVersion<LibraryVersion::PgmathPrecise>(resultType);
+  case LibraryVersion::PgmathRelaxed:
+    return GetHostRuntimeMapVersion<LibraryVersion::PgmathRelaxed>(resultType);
+  case LibraryVersion::PgmathFast:
+    return GetHostRuntimeMapVersion<LibraryVersion::PgmathFast>(resultType);
   }
-  // No long double functions in libpgmath
-  if constexpr (host::FortranTypeExists<long double>()) {
-    AddLibmRealHostProcedures<long double>(lib);
-  }
-  if constexpr (host::FortranTypeExists<std::complex<long double>>()) {
-    AddLibmComplexHostProcedures<long double>(lib);
+  return nullptr;
+
+static const HostRuntimeFunction *SearchInHostRuntimeMap(
+    const HostRuntimeMap &map, const std::string &name, DynamicType resultType,
+    const std::vector<DynamicType> &argTypes) {
+  auto sameNameRange{map.equal_range(name)};
+  for (const auto *iter{sameNameRange.first}; iter != sameNameRange.second;
+       ++iter) {
+    if (iter->resultType == resultType && iter->argumentTypes == argTypes) {
+      return &*iter;
     }
   }
-#endif // LINK_WITH_LIBPGMATH
+  return nullptr;
+}
 
-// Define which host runtime functions will be used for folding
-HostIntrinsicProceduresLibrary::HostIntrinsicProceduresLibrary() {
+// Search host runtime libraries for an exact type match.
+static const HostRuntimeFunction *SearchHostRuntime(const std::string &name,
+    DynamicType resultType, const std::vector<DynamicType> &argTypes) {
   // TODO: When command line options regarding targeted numerical library is
   // available, this needs to be revisited to take it into account. So far,
   // default to libpgmath if F18 is built with it.
 #if LINK_WITH_LIBPGMATH
-  // This looks and is stupid for now (until TODO above), but it is needed
-  // to silence clang warnings on unused symbols if all declared pgmath
-  // symbols are not used somewhere.
-  if (true) {
-    InitHostIntrinsicLibraryWithLibpgmath<L::P>(*this);
-  } else if (false) {
-    InitHostIntrinsicLibraryWithLibpgmath<L::F>(*this);
-  } else {
-    InitHostIntrinsicLibraryWithLibpgmath<L::R>(*this);
+  if (const auto *map{
+          GetHostRuntimeMap(LibraryVersion::PgmathPrecise, resultType)}) {
+    if (const auto *hostFunction{
+            SearchInHostRuntimeMap(*map, name, resultType, argTypes)}) {
+      return hostFunction;
+    }
   }
-#else
-  InitHostIntrinsicLibraryWithLibm(*this);
+  // Default to libm if functions or types are not available in pgmath.
 #endif
+  if (const auto *map{GetHostRuntimeMap(LibraryVersion::Libm, resultType)}) {
+    if (const auto *hostFunction{
+            SearchInHostRuntimeMap(*map, name, resultType, argTypes)}) {
+      return hostFunction;
+    }
+  }
+  return nullptr;
+}
+
+// Return a DynamicType that can hold all values of a given type.
+// This is used to allow 16bit float to be folded with 32bits and
+// x87 float to be folded with IEEE 128bits.
+static DynamicType BiggerType(DynamicType type) {
+  if (type.category() == TypeCategory::Real ||
+      type.category() == TypeCategory::Complex) {
+    // 16 bits floats to IEEE 32 bits float
+    if (type.kind() == common::RealKindForPrecision(11) ||
+        type.kind() == common::RealKindForPrecision(8)) {
+      return {type.category(), common::RealKindForPrecision(24)};
+    }
+    // x87 float to IEEE 128 bits float
+    if (type.kind() == common::RealKindForPrecision(64)) {
+      return {type.category(), common::RealKindForPrecision(113)};
+    }
+  }
+  return type;
+}
+
+std::optional<HostRuntimeWrapper> GetHostRuntimeWrapper(const std::string &name,
+    DynamicType resultType, const std::vector<DynamicType> &argTypes) {
+  if (const auto *hostFunction{SearchHostRuntime(name, resultType, argTypes)}) {
+    return hostFunction->folder;
+  }
+  // If no exact match, search with "bigger" types and insert type
+  // conversions around the folder.
+  std::vector<evaluate::DynamicType> biggerArgTypes;
+  evaluate::DynamicType biggerResultType{BiggerType(resultType)};
+  for (auto type : argTypes) {
+    biggerArgTypes.emplace_back(BiggerType(type));
+  }
+  if (const auto *hostFunction{
+          SearchHostRuntime(name, biggerResultType, biggerArgTypes)}) {
+    return [hostFunction, resultType](
+               FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
+      auto nArgs{args.size()};
+      for (size_t i{0}; i < nArgs; ++i) {
+        args[i] = Fold(context,
+            ConvertToType(hostFunction->argumentTypes[i], std::move(args[i]))
+                .value());
+      }
+      return Fold(context,
+          ConvertToType(
+              resultType, hostFunction->folder(context, std::move(args)))
+              .value());
+    };
+  }
+  return std::nullopt;
 }
 } // namespace Fortran::evaluate

flang/lib/Lower/IntrinsicCall.cpp

 // dialect of MLIR, it makes extensive use of MLIR interfaces and MLIR's coding
 // style (https://mlir.llvm.org/getting_started/DeveloperGuide/) is used in this
 // module.
 //
 //===----------------------------------------------------------------------===//
 
 #include "flang/Lower/IntrinsicCall.h"
 #include "RTBuilder.h"
+#include "flang/Common/static-multimap-view.h"
 #include "flang/Lower/CharacterExpr.h"
 #include "flang/Lower/ComplexExpr.h"
 #include "flang/Lower/ConvertType.h"
 #include "flang/Lower/FIRBuilder.h"
 #include "flang/Lower/Mangler.h"
 #include "flang/Lower/Runtime.h"
 #include "llvm/Support/CommandLine.h"
 #include "llvm/Support/ErrorHandling.h"
 #include <algorithm>
+#include <string_view>
 #include <utility>
 
 #define PGMATH_DECLARE
 #include "../runtime/pgmath.h.inc"
 
 /// This file implements lowering of Fortran intrinsic procedures.
 /// Intrinsics are lowered to a mix of FIR and MLIR operations as
 /// well as call to runtime functions or LLVM intrinsics.
enum class ExtremumBehavior {
   // IEEE minNum/maxNum behavior (754-2008, section 5.3.1):
   // TODO: Not implemented.
   // It is the only behavior where the signaling/quiet aspect of a NaN argument
   // impacts if the result should be NaN or the argument that is a number.
   // LLVM/MLIR do not provide ways to observe this aspect, so it is not
   // possible to implement it without some target dependent runtime.
 };
 
-namespace {
-/// StaticMultimapView is a constexpr friendly multimap
-/// implementation over sorted constexpr arrays.
-/// As the View name suggests, it does not duplicate the
-/// sorted array but only brings range and search concepts
-/// over it. It provides compile time search and can also
-/// provide dynamic search (currently linear, can be improved to
-/// log(n) due to the sorted array property).
-
-// TODO: Find a better place for this if this is retained.
-// This is currently here because this was designed to provide
-// maps over runtime description without the burden of having to
-// instantiate these maps dynamically and to keep them somewhere.
-template <typename V>
-class StaticMultimapView {
-public:
-  using Key = typename V::Key;
-  struct Range {
-    using const_iterator = const V *;
-    constexpr const_iterator begin() const { return startPtr; }
-    constexpr const_iterator end() const { return endPtr; }
-    constexpr bool empty() const {
-      return startPtr == nullptr || endPtr == nullptr || endPtr <= startPtr;
-    }
-    constexpr std::size_t size() const {
-      return empty() ? 0 : static_cast<std::size_t>(endPtr - startPtr);
-    }
-    const V *startPtr{nullptr};
-    const V *endPtr{nullptr};
-  };
-  using const_iterator = typename Range::const_iterator;
-
-  template <std::size_t N>
-  constexpr StaticMultimapView(const V (&array)[N])
-      : range{&array[0], &array[0] + N} {}
-  template <typename Key>
-  constexpr bool verify() {
-    // TODO: sorted
-    // non empty increasing pointer direction
-    return !range.empty();
-  }
-  constexpr const_iterator begin() const { return range.begin(); }
-  constexpr const_iterator end() const { return range.end(); }
-
-  // Assume array is sorted.
-  // TODO make it a log(n) search based on sorted property
-  // std::equal_range will be constexpr in C++20 only.
-  constexpr Range getRange(const Key &key) const {
-    bool matched{false};
-    const V *start{nullptr}, *end{nullptr};
-    for (const auto &desc : range) {
-      if (desc.key == key) {
-        if (!matched) {
-          start = &desc;
-          matched = true;
-        }
-      } else if (matched) {
-        end = &desc;
-        matched = false;
-      }
-    }
-    if (matched) {
-      end = range.end();
-    }
-    return Range{start, end};
-  }
-
-  constexpr std::pair<const_iterator, const_iterator>
-  equal_range(const Key &key) const {
-    Range range{getRange(key)};
-    return {range.begin(), range.end()};
-  }
-
-  constexpr typename Range::const_iterator find(Key key) const {
-    const Range subRange{getRange(key)};
-    return subRange.size() == 1 ? subRange.begin() : end();
-  }
-
-private:
-  Range range{nullptr, nullptr};
-};
-} // namespace
-
 // TODO error handling -> return a code or directly emit messages ?
 struct IntrinsicLibrary {
 
   // Constructors.
   explicit IntrinsicLibrary(Fortran::lower::FirOpBuilder &builder,
                             mlir::Location loc)
       : builder{builder}, loc{loc} {}
   IntrinsicLibrary() = delete;
llvm::cl::values(
         clEnumValN(fastVersion, "fast", "use pgmath fast runtime"),
         clEnumValN(relaxedVersion, "relaxed", "use pgmath relaxed runtime"),
         clEnumValN(preciseVersion, "precise", "use pgmath precise runtime"),
         clEnumValN(llvmOnly, "llvm",
                    "only use LLVM intrinsics (may be incomplete)")),
     llvm::cl::init(fastVersion));
 
 struct RuntimeFunction {
-  using Key = llvm::StringRef;
-  Key key;
+  // llvm::StringRef comparison operator are not constexpr, so use string_view.
+  using Key = std::string_view;
+  // Needed for implicit compare with keys.
+  constexpr operator Key() const { return key; }
+  Key key; // intrinsic name
   llvm::StringRef symbol;
   Fortran::lower::FuncTypeBuilderFunc typeGenerator;
 };
 
 #define RUNTIME_STATIC_DESCRIPTION(name, func)                                 \
   {#name, #func,                                                               \
    Fortran::lower::RuntimeTableKey<decltype(func)>::getTypeModel()},
 static constexpr RuntimeFunction pgmathFast[] = {
auto function = builder.addNamedFunction(
       loc, runtime.symbol, runtime.typeGenerator(builder.getContext()));
   function.setAttr("fir.runtime", builder.getUnitAttr());
   return function;
 }
 
 /// Select runtime function that has the smallest distance to the intrinsic
 /// function type and that will not imply narrowing arguments or extending the
 /// result.
 /// If nothing is found, the mlir::FuncOp will contain a nullptr.
-template <std::size_t N>
-mlir::FuncOp searchFunctionInLibrary(mlir::Location loc,
-                                     Fortran::lower::FirOpBuilder &builder,
-                                     const RuntimeFunction (&lib)[N],
-                                     llvm::StringRef name,
-                                     mlir::FunctionType funcType,
-                                     const RuntimeFunction **bestNearMatch,
-                                     FunctionDistance &bestMatchDistance) {
-  auto map = StaticMultimapView(lib);
-  auto range = map.equal_range(name);
+mlir::FuncOp searchFunctionInLibrary(
+    mlir::Location loc, Fortran::lower::FirOpBuilder &builder,
+    const Fortran::common::StaticMultimapView<RuntimeFunction> &lib,
+    llvm::StringRef name, mlir::FunctionType funcType,
+  auto range = lib.equal_range(name);
   for (auto iter{range.first}; iter != range.second && iter; ++iter) {
     const auto &impl = *iter;
     auto implType = impl.typeGenerator(builder.getContext());
     if (funcType == implType) {
       return getFuncOp(loc, builder, impl); // exact match
     } else {
       FunctionDistance distance(funcType, implType);
       if (distance.isSmallerThan(bestMatchDistance)) {

 /// If nothing is found, the mlir::FuncOp will contain a nullptr.
 static mlir::FuncOp getRuntimeFunction(mlir::Location loc,
                                        Fortran::lower::FirOpBuilder &builder,
                                        llvm::StringRef name,
                                        mlir::FunctionType funcType) {
   const RuntimeFunction *bestNearMatch = nullptr;
   FunctionDistance bestMatchDistance{};
   mlir::FuncOp match;
+  using RtMap = Fortran::common::StaticMultimapView<RuntimeFunction>;
+  static constexpr RtMap pgmathF(pgmathFast);
+  static_assert(pgmathF.Verify() && "map must be sorted");
+  static constexpr RtMap pgmathR(pgmathRelaxed);
+  static_assert(pgmathR.Verify() && "map must be sorted");
+  static constexpr RtMap pgmathP(pgmathPrecise);
+  static_assert(pgmathP.Verify() && "map must be sorted");
   if (mathRuntimeVersion == fastVersion) {
-    match = searchFunctionInLibrary(loc, builder, pgmathFast, name, funcType,
+    match = searchFunctionInLibrary(loc, builder, pgmathF, name, funcType,
                                     &bestNearMatch, bestMatchDistance);
   } else if (mathRuntimeVersion == relaxedVersion) {
-    match = searchFunctionInLibrary(loc, builder, pgmathRelaxed, name, funcType,
+    match = searchFunctionInLibrary(loc, builder, pgmathR, name, funcType,
                                     &bestNearMatch, bestMatchDistance);
   } else if (mathRuntimeVersion == preciseVersion) {
-    match = searchFunctionInLibrary(loc, builder, pgmathPrecise, name, funcType,
+    match = searchFunctionInLibrary(loc, builder, pgmathP, name, funcType,
                                     &bestNearMatch, bestMatchDistance);
   } else {
     assert(mathRuntimeVersion == llvmOnly && "unknown math runtime");
   }
   if (match)
     return match;
 
   // Go through llvm intrinsics if not exact match in libpgmath or if
   // mathRuntimeVersion == llvmOnly
+  static constexpr RtMap llvmIntr(llvmIntrinsics);
+  static_assert(llvmIntr.Verify() && "map must be sorted");
   if (auto exactMatch =
-          searchFunctionInLibrary(loc, builder, llvmIntrinsics, name, funcType,
+          searchFunctionInLibrary(loc, builder, llvmIntr, name, funcType,
                                   &bestNearMatch, bestMatchDistance))
     return exactMatch;
 
   if (bestNearMatch != nullptr) {
     assert(!bestMatchDistance.isLosingPrecision() &&
            "runtime selection loses precision");
     return getFuncOp(loc, builder, *bestNearMatch);
   }

flang/runtime/pgmath.h.inc

 PGMATH_MTH_VERSION_REAL(acosh)
 PGMATH_ALL(asin)
 PGMATH_MTH_VERSION_REAL(asinh)
 PGMATH_ALL(atan)
 PGMATH_REAL2(atan2)
 PGMATH_MTH_VERSION_REAL(atanh)
 PGMATH_MTH_VERSION_REAL(bessel_j0)
 PGMATH_MTH_VERSION_REAL(bessel_j1)
-PGMATH_MTH_VERSION_REAL(bessel_y0)
-PGMATH_MTH_VERSION_REAL(bessel_y1)
 // bessel_jn and bessel_yn takes an int as first arg
 PGMATH_DECLARE(float __mth_i_bessel_jn(int, float))
 PGMATH_DECLARE(double __mth_i_dbessel_jn(int, double))
 PGMATH_USE_S(bessel_jn, __mth_i_bessel_jn)
 PGMATH_USE_D(bessel_jn, __mth_i_dbessel_jn)
+PGMATH_MTH_VERSION_REAL(bessel_y0)
+PGMATH_MTH_VERSION_REAL(bessel_y1)
 PGMATH_DECLARE(float __mth_i_bessel_yn(int, float))
 PGMATH_DECLARE(double __mth_i_dbessel_yn(int, double))
 PGMATH_USE_S(bessel_yn, __mth_i_bessel_yn)
 PGMATH_USE_D(bessel_yn, __mth_i_dbessel_yn)
 
 PGMATH_ALL(cos)
 PGMATH_ALL(cosh)
 PGMATH_MTH_VERSION_REAL(erf)

flang/test/Evaluate/folding02.f90

 ! RUN: %S/test_folding.sh %s %t %f18
 ! Check intrinsic function folding with host runtime library
 
 module m
+  real(2), parameter :: eps2 = 0.001_2
+  real(2), parameter :: eps3 = 0.001_3
   real(4), parameter :: eps4 = 0.000001_4
   real(8), parameter :: eps8 = 0.000000000000001_8
 
   ! These eps have been set so that constant folding of intrinsic
   ! functions that use host runtime can be tested independently of
   ! the underlying math library used.
   ! C++ <cmath> and libpgmath precise, relaxed and fast libraries pass the test.
   ! It may have to be relaxed to pass on all architectures.
   ! The purpose is to check that the "correct" runtime functions are mapped
   ! to intrinsic functions but not to test the stability between different libraries.
   ! The eps should just be low enough to give confidence that intrinsic
   ! functions are mapped to runtime functions implementing the same math
   ! function.
   ! Compared values were selected to be around 1 +/- 0.5 so that eps is meaningful.
   ! Expected values come from libpgmath-precise for Real(4) and Real(8) and
   ! were computed on X86_64.
 
-! Real scalar intrinsic function tests
-
-  #define TEST_R4(name, result, expected) \
-  real(kind=4), parameter :: res_##name##_r4 = result; \
-  real(kind=4), parameter :: exp_##name##_r4 = expected; \
-  logical, parameter :: test_##name##_r4 = abs(res_##name##_r4 - exp_##name##_r4).LE.(eps4)
-
   logical, parameter :: test_sign_i4 = sign(1_4,2_4) == 1_4 .and. sign(1_4,-3_4) == -1_4
   logical, parameter :: test_sign_i8 = sign(1_8,2_8) == 1_8 .and. sign(1_8,-3_8) == -1_8
 
+! Real scalar intrinsic function tests
+#define TEST_FLOATING(name, result, expected, t, k) \
+  t(kind = k), parameter ::res_##name##_##t##k = result; \
+  t(kind = k), parameter ::exp_##name##_##t##k = expected; \
+  logical, parameter ::test_##name##_##t##k = abs(res_##name##_##t##k - exp_##name##_##t##k).LE.(eps##k)
+
+#define TEST_R2(name, result, expected) TEST_FLOATING(name, result, expected, real, 2)
+#define TEST_R3(name, result, expected) TEST_FLOATING(name, result, expected, real, 3)
+#define TEST_R4(name, result, expected) TEST_FLOATING(name, result, expected, real, 4)
+#define TEST_R8(name, result, expected) TEST_FLOATING(name, result, expected, real, 8)
+#define TEST_C4(name, result, expected) TEST_FLOATING(name, result, expected, complex, 4)
+#define TEST_C8(name, result, expected) TEST_FLOATING(name, result, expected, complex, 8)
+
+! REAL(4) tests.
+
   logical, parameter :: test_abs_r4 = abs(-2._4).EQ.(2._4)
   TEST_R4(acos, acos(0.5_4), 1.0471975803375244140625_4)
   TEST_R4(acosh, acosh(1.5_4), 0.96242368221282958984375_4)
   logical, parameter :: test_aint1 = aint(2.783).EQ.(2.)
   logical, parameter :: test_anint1 = anint(2.783).EQ.(3.)
   logical, parameter :: test_floor1 = floor(-2.783).EQ.(-3.)
   logical, parameter :: test_floor2 = floor(2.783).EQ.(2.)
   logical, parameter :: test_ceiling1 = ceiling(-2.783).EQ.(-2.)
! REAL(4) tests.
   TEST_R4(real, real(z'3f800000'), 1._4)
   logical, parameter :: test_sign_r4 = sign(1._4,2._4) == 1._4 .and. sign(1._4,-2._4) == -1._4
   TEST_R4(sin, sin(1.6_4), 0.99957358837127685546875_4)
   TEST_R4(sinh, sinh(0.9_4), 1.0265166759490966796875_4)
   TEST_R4(sqrt, sqrt(1.1_4), 1.0488088130950927734375_4)
   TEST_R4(tan, tan(0.8_4), 1.0296385288238525390625_4)
   TEST_R4(tanh, tanh(3._4), 0.995054781436920166015625_4)
 
-! Real(kind=8) tests.
-
-  #define TEST_R8(name, result, expected) \
-  real(kind=8), parameter :: res_##name##_r8 = result; \
-  real(kind=8), parameter :: exp_##name##_r8 = expected; \
-  logical, parameter :: test_##name##_r8 = abs(res_##name##_r8 - exp_##name##_r8).LE.(eps8)
+! REAL(8) tests.
 
   logical, parameter :: test_abs_r8 = abs(-2._8).EQ.(2._8)
   TEST_R8(acos, acos(0.5_8), &
     1.047197551196597853362391106202267110347747802734375_8)
   TEST_R8(acosh, acosh(1.5_8), &
     0.9624236501192069415111518537742085754871368408203125_8)
   TEST_R8(asin, asin(0.9_8), &
     1.119769514998634196700777465594001114368438720703125_8)
TEST_R8(sinh, sinh(0.9_8), &
     1.0265167257081753149350333842448890209197998046875_8)
   TEST_R8(sqrt, sqrt(1.1_8), &
     1.048808848170151630796453900984488427639007568359375_8)
   TEST_R8(tan, tan(0.8_8), &
     1.0296385570503641115891468871268443763256072998046875_8)
   TEST_R8(tanh, tanh(3._8), &
     0.995054753686730464323773048818111419677734375_8)
 
-  #define TEST_C4(name, result, expected) \
-  complex(kind=4), parameter :: res_##name##_c4 = result; \
-  complex(kind=4), parameter :: exp_##name##_c4 = expected; \
-  logical, parameter :: test_##name##_c4 = abs(res_##name##_c4 - exp_##name##_c4).LE.(eps4)
+! COMPLEX(4) tests.
 
   logical, parameter :: test_abs_c4 = abs(abs((1.1_4, 0.1_4)) &
     - 1.10453617572784423828125_4).LE.(eps4)
   TEST_C4(acos, acos((0.7_4, 1.1_4)), &
     (1.11259567737579345703125_4, -1.03283786773681640625_4))
   TEST_C4(acosh, acosh((0.7_4, 1.1_4)), &
     (1.03283774852752685546875_4, 1.11259555816650390625_4))
   TEST_C4(asin, asin((1.4_4, 0.7_4)), &
TEST_C4(sinh, sinh((1.1_4, 0.7_4)), &
     (1.02155959606170654296875_4,1.07488918304443359375_4))
   TEST_C4(sqrt, sqrt((0.1_4, 2.1_4)), &
     (1.04937589168548583984375_4,1.00059473514556884765625_4))
   TEST_C4(tan, tan((1.1_4, 0.4_4)), &
     (1.07952976226806640625_4,1.1858270168304443359375_4))
   TEST_C4(tanh, tanh((0.4_4, 1.1_4)), &
     (1.1858270168304443359375_4,1.07952976226806640625_4))
 
-  #define TEST_C8(name, result, expected) \
-  complex(kind=8), parameter :: res_##name##_c8 = result; \
-  complex(kind=8), parameter :: exp_##name##_c8 = expected; \
-  logical, parameter :: test_##name##_c8 = abs(res_##name##_c8 - exp_##name##_c8).LE.(eps8)
+! COMPLEX(8) tests.
 
   logical, parameter :: test_abs_c8 = abs(abs((1.1_8, 0.1_8)) &
     - 1.1045361017187260710414875575224868953227996826171875_8).LE.(eps4)
   TEST_C8(acos, acos((0.7_8, 1.1_8)), &
     (1.1125956244800556671492586247040890157222747802734375_8, &
       (-1.032837729564676454430127705563791096210479736328125_8)))
   TEST_C8(acosh, acosh((0.7_8, 1.1_8)), &
     (1.0328377295646762323855227805324830114841461181640625_8, &
(1.04937591075907210580453465809114277362823486328125_8, &
       (1.0005947241922830059472419228351358112816260614863494993187487125396728515625_8)))
   TEST_C8(tan, tan((1.1_8, 0.4_8)), &
     (1.07952982287592025301137255155481398105621337890625_8, &
       (1.1858270353667335061942367246956564486026763916015625_8)))
   TEST_C8(tanh, tanh((0.4_8, 1.1_8)), &
     (1.1858270353667335061942367246956564486026763916015625_8, &
       (1.07952982287592025301137255155481398105621337890625_8)))
 
+
+  ! Only test a few REAL(2)/REAL(3) cases since they anyway use the real 4
+  ! runtime mapping.
+  TEST_R2(acos, acos(0.5_2), 1.046875_2)
+  TEST_R2(atan2, atan2(1.5_2, 1._2), 9.8291015625e-1_2)
+
+  TEST_R3(acos, acos(0.5_3), 1.046875_3)
+  TEST_R3(atan2, atan2(1.3_2, 1._3), 9.140625e-1_3)
+
 #ifdef TEST_LIBPGMATH
 ! Bessel functions and erfc_scaled can only be folded if libpgmath
 ! is used.
   TEST_R4(bessel_j0, bessel_j0(0.5_4), 0.938469827175140380859375_4)
   TEST_R4(bessel_j1, bessel_j1(1.8_4), 0.5815169811248779296875_4)
   TEST_R4(bessel_jn, bessel_jn(2, 3._4), 0.4860912859439849853515625_4)
   TEST_R4(bessel_y0, bessel_y0(2._4), 0.510375678539276123046875_4)
   TEST_R4(bessel_y1, bessel_y1(1._4), (-0.78121280670166015625_4))

flang/unittests/Evaluate/folding.cpp

 #include "testing.h"
 #include "../../lib/Evaluate/host.h"
-#include "../../lib/Evaluate/intrinsics-library-templates.h"
 #include "flang/Evaluate/call.h"
 #include "flang/Evaluate/expression.h"
 #include "flang/Evaluate/fold.h"
+#include "flang/Evaluate/intrinsics-library.h"
 #include "flang/Evaluate/intrinsics.h"
 #include "flang/Evaluate/tools.h"
 #include <tuple>
 
 using namespace Fortran::evaluate;
 
 // helper to call functions on all types from tuple
 template <typename... T> struct RunOnTypes {};
template <typename T> static void Run() {
     Expr<SomeType> exprSomeType{exprSomeKind};
     TEST(GetScalarConstantValue<T>(exprFullyTyped).has_value());
     TEST(GetScalarConstantValue<T>(exprSomeKind).has_value());
     TEST(GetScalarConstantValue<T>(exprSomeType).has_value());
   }
 };
 
 template <typename T>
-static FunctionRef<T> CreateIntrinsicElementalCall(
-    const std::string &name, const Expr<T> &arg) {
-  Fortran::semantics::Attrs attrs;
-  attrs.set(Fortran::semantics::Attr::ELEMENTAL);
-  ActualArguments args{ActualArgument{AsGenericExpr(arg)}};
-  ProcedureDesignator intrinsic{
-      SpecificIntrinsic{name, T::GetType(), 0, attrs}};
-  return FunctionRef<T>{std::move(intrinsic), std::move(args)};
-}
-
-// Test flushSubnormalsToZero when folding with host runtime.
-// Subnormal value flushing on host is handle in host.cpp
-// HostFloatingPointEnvironment::SetUpHostFloatingPointEnvironment
-
-// Dummy host runtime functions where subnormal flushing matters
-float SubnormalFlusher1(float f) { // given f is subnormal
-  return 2.3 * f; // returns 0 if subnormal arguments are flushed to zero
-}
-
-float SubnormalFlusher2(float f) { // given f/2 is subnormal
-  return f / 2.3; // returns 0 if subnormal
+Scalar<T> CallHostRt(
+    HostRuntimeWrapper func, FoldingContext &context, Scalar<T> x) {
+  return GetScalarConstantValue<T>(
+      func(context, {AsGenericExpr(Constant<T>{x})}))
+      .value();
 }
 
 void TestHostRuntimeSubnormalFlushing() {
   using R4 = Type<TypeCategory::Real, 4>;
   if constexpr (std::is_same_v<host::HostType<R4>, float>) {
     Fortran::parser::CharBlock src;
     Fortran::parser::ContextualMessages messages{src, nullptr};
     Fortran::common::IntrinsicTypeDefaultKinds defaults;
     auto intrinsics{Fortran::evaluate::IntrinsicProcTable::Configure(defaults)};
     FoldingContext flushingContext{
         messages, defaults, intrinsics, defaultRounding, true};
     FoldingContext noFlushingContext{
         messages, defaults, intrinsics, defaultRounding, false};
 
-    HostIntrinsicProceduresLibrary lib;
-    lib.AddProcedure(HostRuntimeIntrinsicProcedure{
-        "flusher_test1", SubnormalFlusher1, true});
-    lib.AddProcedure(HostRuntimeIntrinsicProcedure{
-        "flusher_test2", SubnormalFlusher2, true});
-
+    DynamicType r4{R4{}.GetType()};
     // Test subnormal argument flushing
-    if (auto callable{
-            lib.GetHostProcedureWrapper<Scalar, R4, R4>("flusher_test1")}) {
+    if (auto callable{GetHostRuntimeWrapper("log", r4, {r4})}) {
       // Biggest IEEE 32bits subnormal power of two
-      host::HostType<R4> input1{5.87747175411144e-39};
-      const Scalar<R4> x1{host::CastHostToFortran<R4>(input1)};
-      Scalar<R4> y1Flushing{callable.value()(flushingContext, x1)};
-      Scalar<R4> y1NoFlushing{callable.value()(noFlushingContext, x1)};
-      TEST(y1Flushing.IsZero());
-      TEST(!y1NoFlushing.IsZero());
-    } else {
-      TEST(false);
-    }
-    // Test subnormal result flushing
-    if (auto callable{
-            lib.GetHostProcedureWrapper<Scalar, R4, R4>("flusher_test2")}) {
-      // Smallest (positive) non-subnormal IEEE 32 bit float value
-      host::HostType<R4> input2{1.1754944e-38};
-      const Scalar<R4> x2{host::CastHostToFortran<R4>(input2)};
-      Scalar<R4> y2Flushing{callable.value()(flushingContext, x2)};
-      Scalar<R4> y2NoFlushing{callable.value()(noFlushingContext, x2)};
-      TEST(y2Flushing.IsZero());
-      TEST(!y2NoFlushing.IsZero());
+      const Scalar<R4> x1{Scalar<R4>::Word{0x00400000}};
+      Scalar<R4> y1Flushing{CallHostRt<R4>(*callable, flushingContext, x1)};
+      Scalar<R4> y1NoFlushing{CallHostRt<R4>(*callable, noFlushingContext, x1)};
+      // We would expect y1Flushing to be NaN, but some libc logf implementation
+      // "workaround" subnormal flushing by returning a constant negative
+      // results for all subnormal values (-1.03972076416015625e2_4). In case of
+      // flushing, the result should still be different than -88 +/- 2%.
+      TEST(y1Flushing.IsInfinite() ||
+          std::abs(host::CastFortranToHost<R4>(y1Flushing) + 88.) > 2);
+      TEST(!y1NoFlushing.IsInfinite() &&
+          std::abs(host::CastFortranToHost<R4>(y1NoFlushing) + 88.) < 2);
     } else {
       TEST(false);
     }
   } else {
     TEST(false); // Cannot run this test on the host
   }
 }
 
 int main() {
   RunOnTypes<TestGetScalarConstantValue, AllIntrinsicTypes>::Run();
   TestHostRuntimeSubnormalFlushing();
   return testing::Complete();
 }