diff --git a/flang/include/flang/Common/static-multimap-view.h b/flang/include/flang/Common/static-multimap-view.h
new file mode 100644
--- /dev/null
+++ b/flang/include/flang/Common/static-multimap-view.h
@@ -0,0 +1,65 @@
+//===-- 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 *;
+
+  template <typename T, const T &array>
+  static constexpr StaticMultimapView Create() {
+    static_assert(StaticMultimapView{array}.Verify(),
+        "array must be sorted in increasing order");
+    return StaticMultimapView{array};
+  }
+  constexpr const_iterator begin() const { return begin_; }
+  constexpr const_iterator end() const { return end_; }
+
+  // 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);
+  }
+
+private:
+  template <std::size_t N>
+  constexpr StaticMultimapView(const V (&array)[N])
+      : begin_{&array[0]}, end_{&array[0] + N} {}
+
+  // 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_};
+    for (const auto *x{begin_}; x != end_; ++x) {
+      if (x->key < lastSeen->key) {
+        return false;
+      }
+      lastSeen = x;
+    }
+    return true;
+  }
+  const_iterator begin_{nullptr};
+  const_iterator end_{nullptr};
+};
+} // namespace Fortran::common
+#endif // FORTRAN_COMMON_STATIC_MULTIMAP_VIEW_H_
diff --git a/flang/include/flang/Evaluate/common.h b/flang/include/flang/Evaluate/common.h
--- a/flang/include/flang/Evaluate/common.h
+++ b/flang/include/flang/Evaluate/common.h
@@ -9,7 +9,6 @@
 #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"
@@ -237,9 +236,6 @@
   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);
@@ -264,7 +260,6 @@
   bool bigEndian_{false};
   const semantics::DerivedTypeSpec *pdtInstance_{nullptr};
   std::map<parser::CharBlock, ConstantSubscript> impliedDos_;
-  HostIntrinsicProceduresLibrary hostIntrinsicsLibrary_;
 };
 
 void RealFlagWarnings(FoldingContext &, const RealFlags &, const char *op);
diff --git a/flang/include/flang/Evaluate/intrinsics-library.h b/flang/include/flang/Evaluate/intrinsics-library.h
--- a/flang/include/flang/Evaluate/intrinsics-library.h
+++ b/flang/include/flang/Evaluate/intrinsics-library.h
@@ -10,99 +10,36 @@
 #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_
diff --git a/flang/lib/Evaluate/fold-complex.cpp b/flang/lib/Evaluate/fold-complex.cpp
--- a/flang/lib/Evaluate/fold-complex.cpp
+++ b/flang/lib/Evaluate/fold-complex.cpp
@@ -23,8 +23,7 @@
       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 {
diff --git a/flang/lib/Evaluate/fold-implementation.h b/flang/lib/Evaluate/fold-implementation.h
--- a/flang/lib/Evaluate/fold-implementation.h
+++ b/flang/lib/Evaluate/fold-implementation.h
@@ -12,7 +12,6 @@
 #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"
@@ -22,6 +21,7 @@
 #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"
@@ -70,6 +70,24 @@
 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
@@ -1410,8 +1428,7 @@
       }
       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 {
diff --git a/flang/lib/Evaluate/fold-real.cpp b/flang/lib/Evaluate/fold-real.cpp
--- a/flang/lib/Evaluate/fold-real.cpp
+++ b/flang/lib/Evaluate/fold-real.cpp
@@ -29,8 +29,7 @@
       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 {
@@ -44,9 +43,7 @@
       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 {
@@ -58,9 +55,7 @@
     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 {
@@ -75,9 +70,7 @@
       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 {
diff --git a/flang/lib/Evaluate/host.h b/flang/lib/Evaluate/host.h
--- a/flang/lib/Evaluate/host.h
+++ b/flang/lib/Evaluate/host.h
@@ -95,22 +95,6 @@
   }
 }
 
-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;
@@ -139,21 +123,27 @@
 // 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 &&
@@ -161,7 +151,9 @@
       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 &&
@@ -211,49 +203,6 @@
   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
 
diff --git a/flang/lib/Evaluate/intrinsics-library-templates.h b/flang/lib/Evaluate/intrinsics-library-templates.h
deleted file mode 100644
--- a/flang/lib/Evaluate/intrinsics-library-templates.h
+++ /dev/null
@@ -1,209 +0,0 @@
-//===-- 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_
diff --git a/flang/lib/Evaluate/intrinsics-library.cpp b/flang/lib/Evaluate/intrinsics-library.cpp
--- a/flang/lib/Evaluate/intrinsics-library.cpp
+++ b/flang/lib/Evaluate/intrinsics-library.cpp
@@ -8,74 +8,231 @@
 
 // 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;
     }
   }
-  return false;
+  return true;
 }
 
-// Map numerical intrinsic to  <cmath>/<complex> functions
+// 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;
+};
 
-// Define which host runtime functions will be used for folding
+// 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);
+  }
+  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...>{});
+}
+
+// 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 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
 
@@ -88,313 +245,260 @@
   // 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 auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
 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},
+  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 auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
 
-  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);
-  }
-}
-
+/// 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},
-#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},
+#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 {
-    static_assert(Lib == L::P && "unexpected libpgmath version");
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
-#define PGMATH_PRECISE
-#define PGMATH_USE_S(name, function) {#name, function, true},
-#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[]{
+  };
+  static constexpr auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
+template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathFast> {
+  static constexpr HostRuntimeFunction table[]{
 #define PGMATH_FAST
-#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 if constexpr (Lib == L::R) {
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
+  };
+  static constexpr auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
+template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathRelaxed> {
+  static constexpr HostRuntimeFunction table[]{
 #define PGMATH_RELAXED
-#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));
-    }
-  } else {
-    static_assert(Lib == L::P && "unexpected libpgmath version");
-    HostRuntimeIntrinsicProcedure pgmathSymbols[]{
+  };
+  static constexpr auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
+template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathRelaxed> {
+  static constexpr HostRuntimeFunction table[]{
+#define PGMATH_RELAXED
+#define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
+#include "../runtime/pgmath.h.inc"
+  };
+  static constexpr auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
+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 auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
+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 auto map{HostRuntimeMap::Create<decltype(table), table>()};
+};
 
-// 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));
+  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;
     }
-  } 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));
+    if (const auto *map{GetHostRuntimeMapHelper<double, version>(resultType)}) {
+      return map;
     }
-  } 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));
+    if (const auto *map{
+            GetHostRuntimeMapHelper<long double, 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 (resultType.category() == TypeCategory::Complex) {
+    if (const auto *map{GetHostRuntimeMapHelper<std::complex<float>, 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<std::complex<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 (const auto *map{
+            GetHostRuntimeMapHelper<std::complex<long double>, 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});
+  return nullptr;
 }
-
-template <L Lib>
-static void InitHostIntrinsicLibraryWithLibpgmath(
-    HostIntrinsicProceduresLibrary &lib) {
-  if constexpr (host::FortranTypeExists<float>()) {
-    AddLibpgmathFloatHostProcedures<Lib>(lib);
-  }
-  if constexpr (host::FortranTypeExists<double>()) {
-    AddLibpgmathDoubleHostProcedures<Lib>(lib);
+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);
   }
-  if constexpr (host::FortranTypeExists<std::complex<float>>()) {
-    AddLibpgmathComplexHostProcedures<Lib>(lib);
-  }
-  if constexpr (host::FortranTypeExists<std::complex<double>>()) {
-    AddLibpgmathDoubleComplexHostProcedures<Lib>(lib);
-  }
-  // 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;
+    }
   }
+  return nullptr;
 }
-#endif // LINK_WITH_LIBPGMATH
 
-// 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
diff --git a/flang/lib/Lower/IntrinsicCall.cpp b/flang/lib/Lower/IntrinsicCall.cpp
--- a/flang/lib/Lower/IntrinsicCall.cpp
+++ b/flang/lib/Lower/IntrinsicCall.cpp
@@ -15,6 +15,7 @@
 
 #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"
@@ -24,6 +25,7 @@
 #include "llvm/Support/CommandLine.h"
 #include "llvm/Support/ErrorHandling.h"
 #include <algorithm>
+#include <string_view>
 #include <utility>
 
 #define PGMATH_DECLARE
@@ -85,89 +87,6 @@
   // 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 {
 
@@ -349,8 +268,11 @@
     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;
 };
@@ -583,16 +505,13 @@
 /// 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,
+    const RuntimeFunction **bestNearMatch,
+    FunctionDistance &bestMatchDistance) {
+  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());
@@ -620,14 +539,21 @@
   const RuntimeFunction *bestNearMatch = nullptr;
   FunctionDistance bestMatchDistance{};
   mlir::FuncOp match;
+  using RtMap = Fortran::common::StaticMultimapView<RuntimeFunction>;
+  static constexpr auto pgmathF =
+      RtMap::Create<decltype(pgmathFast), pgmathFast>();
+  static constexpr auto pgmathR =
+      RtMap::Create<decltype(pgmathRelaxed), pgmathRelaxed>();
+  static constexpr auto pgmathP =
+      RtMap::Create<decltype(pgmathPrecise), pgmathPrecise>();
   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");
@@ -637,8 +563,10 @@
 
   // Go through llvm intrinsics if not exact match in libpgmath or if
   // mathRuntimeVersion == llvmOnly
+  static constexpr auto llvmIntr =
+      RtMap::Create<decltype(llvmIntrinsics), llvmIntrinsics>();
   if (auto exactMatch =
-          searchFunctionInLibrary(loc, builder, llvmIntrinsics, name, funcType,
+          searchFunctionInLibrary(loc, builder, llvmIntr, name, funcType,
                                   &bestNearMatch, bestMatchDistance))
     return exactMatch;
 
diff --git a/flang/runtime/pgmath.h.inc b/flang/runtime/pgmath.h.inc
--- a/flang/runtime/pgmath.h.inc
+++ b/flang/runtime/pgmath.h.inc
@@ -167,13 +167,13 @@
 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)
diff --git a/flang/test/Evaluate/folding02.f90 b/flang/test/Evaluate/folding02.f90
--- a/flang/test/Evaluate/folding02.f90
+++ b/flang/test/Evaluate/folding02.f90
@@ -2,6 +2,8 @@
 ! 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
 
@@ -19,15 +21,23 @@
   ! Expected values come from libpgmath-precise for Real(4) and Real(8) and
   ! were computed on X86_64.
 
+  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_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)
+#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)
 
-  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(4) tests.
 
   logical, parameter :: test_abs_r4 = abs(-2._4).EQ.(2._4)
   TEST_R4(acos, acos(0.5_4), 1.0471975803375244140625_4)
@@ -63,12 +73,7 @@
   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), &
@@ -122,10 +127,7 @@
   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)
@@ -161,10 +163,7 @@
   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)
@@ -215,6 +214,15 @@
     (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.
diff --git a/flang/unittests/Evaluate/folding.cpp b/flang/unittests/Evaluate/folding.cpp
--- a/flang/unittests/Evaluate/folding.cpp
+++ b/flang/unittests/Evaluate/folding.cpp
@@ -1,9 +1,9 @@
 #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>
@@ -30,27 +30,11 @@
 };
 
 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() {
@@ -65,35 +49,21 @@
     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);
     }