flang/lib/Evaluate/intrinsics-library.cpp

   1 //===-- lib/Evaluate/intrinsics-library.cpp -------------------------------===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 // This file defines host runtime functions that can be used for folding
  10 // intrinsic functions.
  11 // The default host runtime folders are built with <cmath> and
  12 // <complex> functions that are guaranteed to exist from the C++ standard.
  13
  14 #include "flang/Evaluate/intrinsics-library.h"
  15 #include "fold-implementation.h"
  16 #include "host.h"
  17 #include "flang/Common/erfc-scaled.h"
  18 #include "flang/Common/idioms.h"
  19 #include "flang/Common/static-multimap-view.h"
  20 #include "flang/Evaluate/expression.h"
  21 #include <cfloat>
  22 #include <cmath>
  23 #include <complex>
  24 #include <functional>
  25 #if HAS_QUADMATHLIB
  26 #include "quadmath.h"
  27 #endif
  28 #include "flang/Common/float128.h"
  29 #include "flang/Common/float80.h"
  30 #include <type_traits>
  31
  32 namespace Fortran::evaluate {
  33
  34 // Define a vector like class that can hold an arbitrary number of
  35 // Dynamic type and be built at compile time. This is like a
  36 // std::vector<DynamicType>, but constexpr only.
  37 template <typename... FortranType> struct TypeVectorStorage {
  38   static constexpr DynamicType values[]{FortranType{}.GetType()...};
  39   static constexpr const DynamicType *start{&values[0]};
  40   static constexpr const DynamicType *end{start + sizeof...(FortranType)};
  41 };
  42 template <> struct TypeVectorStorage<> {
  43   static constexpr const DynamicType *start{nullptr}, *end{nullptr};
  44 };
  45 struct TypeVector {
  46   template <typename... FortranType> static constexpr TypeVector Create() {
  47     using storage = TypeVectorStorage<FortranType...>;
  48     return TypeVector{storage::start, storage::end, sizeof...(FortranType)};
  49   }
  50   constexpr size_t size() const { return size_; };
  51   using const_iterator = const DynamicType *;
  52   constexpr const_iterator begin() const { return startPtr; }
  53   constexpr const_iterator end() const { return endPtr; }
  54   const DynamicType &operator[](size_t i) const { return *(startPtr + i); }
  55
  56   const DynamicType *startPtr{nullptr};
  57   const DynamicType *endPtr{nullptr};
  58   const size_t size_;
  59 };
  60 inline bool operator==(
  61     const TypeVector &lhs, const std::vector<DynamicType> &rhs) {
  62   if (lhs.size() != rhs.size()) {
  63     return false;
  64   }
  65   for (size_t i{0}; i < lhs.size(); ++i) {
  66     if (lhs[i] != rhs[i]) {
  67       return false;
  68     }
  69   }
  70   return true;
  71 }
  72
  73 // HostRuntimeFunction holds a pointer to a Folder function that can fold
  74 // a Fortran scalar intrinsic using host runtime functions (e.g libm).
  75 // The folder take care of all conversions between Fortran types and the related
  76 // host types as well as setting and cleaning-up the floating point environment.
  77 // HostRuntimeFunction are intended to be built at compile time (members are all
  78 // constexpr constructible) so that they can be stored in a compile time static
  79 // map.
  80 struct HostRuntimeFunction {
  81   using Folder = Expr<SomeType> (*)(
  82       FoldingContext &, std::vector<Expr<SomeType>> &&);
  83   using Key = std::string_view;
  84   // Needed for implicit compare with keys.
  85   constexpr operator Key() const { return key; }
  86   // Name of the related Fortran intrinsic.
  87   Key key;
  88   // DynamicType of the Expr<SomeType> returns by folder.
  89   DynamicType resultType;
  90   // DynamicTypes expected for the Expr<SomeType> arguments of the folder.
  91   // The folder will crash if provided arguments of different types.
  92   TypeVector argumentTypes;
  93   // Folder to be called to fold the intrinsic with host runtime. The provided
  94   // Expr<SomeType> arguments must wrap scalar constants of the type described
  95   // in argumentTypes, otherwise folder will crash. Any floating point issue
  96   // raised while executing the host runtime will be reported in FoldingContext
  97   // messages.
  98   Folder folder;
  99 };
 100
 101 // Translate a host function type signature (template arguments) into a
 102 // constexpr data representation based on Fortran DynamicType that can be
 103 // stored.
 104 template <typename TR, typename... TA> using FuncPointer = TR (*)(TA...);
 105 template <typename T> struct FuncTypeAnalyzer {};
 106 template <typename HostTR, typename... HostTA>
 107 struct FuncTypeAnalyzer<FuncPointer<HostTR, HostTA...>> {
 108   static constexpr DynamicType result{host::FortranType<HostTR>{}.GetType()};
 109   static constexpr TypeVector arguments{
 110       TypeVector::Create<host::FortranType<HostTA>...>()};
 111 };
 112
 113 // Define helpers to deal with host floating environment.
 114 template <typename TR>
 115 static void CheckFloatingPointIssues(
 116     host::HostFloatingPointEnvironment &hostFPE, const Scalar<TR> &x) {
 117   if constexpr (TR::category == TypeCategory::Complex ||
 118       TR::category == TypeCategory::Real) {
 119     if (x.IsNotANumber()) {
 120       hostFPE.SetFlag(RealFlag::InvalidArgument);
 121     } else if (x.IsInfinite()) {
 122       hostFPE.SetFlag(RealFlag::Overflow);
 123     }
 124   }
 125 }
 126 // Software Subnormal Flushing helper.
 127 // Only flush floating-points. Forward other scalars untouched.
 128 // Software flushing is only performed if hardware flushing is not available
 129 // because it may not result in the same behavior as hardware flushing.
 130 // Some runtime implementations are "working around" subnormal flushing to
 131 // return results that they deem better than returning the result they would
 132 // with a null argument. An example is logf that should return -inf if arguments
 133 // are flushed to zero, but some implementations return -1.03972076416015625e2_4
 134 // for all subnormal values instead. It is impossible to reproduce this with the
 135 // simple software flushing below.
 136 template <typename T>
 137 static constexpr inline const Scalar<T> FlushSubnormals(Scalar<T> &&x) {
 138   if constexpr (T::category == TypeCategory::Real ||
 139       T::category == TypeCategory::Complex) {
 140     return x.FlushSubnormalToZero();
 141   }
 142   return x;
 143 }
 144
 145 // This is the kernel called by all HostRuntimeFunction folders, it convert the
 146 // Fortran Expr<SomeType> to the host runtime function argument types, calls
 147 // the runtime function, and wrap back the result into an Expr<SomeType>.
 148 // It deals with host floating point environment set-up and clean-up.
 149 template <typename FuncType, typename TR, typename... TA, size_t... I>
 150 static Expr<SomeType> ApplyHostFunctionHelper(FuncType func,
 151     FoldingContext &context, std::vector<Expr<SomeType>> &&args,
 152     std::index_sequence<I...>) {
 153   host::HostFloatingPointEnvironment hostFPE;
 154   hostFPE.SetUpHostFloatingPointEnvironment(context);
 155   host::HostType<TR> hostResult{};
 156   Scalar<TR> result{};
 157   std::tuple<Scalar<TA>...> scalarArgs{
 158       GetScalarConstantValue<TA>(args[I]).value()...};
 159   if (context.targetCharacteristics().areSubnormalsFlushedToZero() &&
 160       !hostFPE.hasSubnormalFlushingHardwareControl()) {
 161     hostResult = func(host::CastFortranToHost<TA>(
 162         FlushSubnormals<TA>(std::move(std::get<I>(scalarArgs))))...);
 163     result = FlushSubnormals<TR>(host::CastHostToFortran<TR>(hostResult));
 164   } else {
 165     hostResult = func(host::CastFortranToHost<TA>(std::get<I>(scalarArgs))...);
 166     result = host::CastHostToFortran<TR>(hostResult);
 167   }
 168   if (!hostFPE.hardwareFlagsAreReliable()) {
 169     CheckFloatingPointIssues<TR>(hostFPE, result);
 170   }
 171   hostFPE.CheckAndRestoreFloatingPointEnvironment(context);
 172   return AsGenericExpr(Constant<TR>(std::move(result)));
 173 }
 174 template <typename HostTR, typename... HostTA>
 175 Expr<SomeType> ApplyHostFunction(FuncPointer<HostTR, HostTA...> func,
 176     FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
 177   return ApplyHostFunctionHelper<decltype(func), host::FortranType<HostTR>,
 178       host::FortranType<HostTA>...>(
 179       func, context, std::move(args), std::index_sequence_for<HostTA...>{});
 180 }
 181
 182 // FolderFactory builds a HostRuntimeFunction for the host runtime function
 183 // passed as a template argument.
 184 // Its static member function "fold" is the resulting folder. It captures the
 185 // host runtime function pointer and pass it to the host runtime function folder
 186 // kernel.
 187 template <typename HostFuncType, HostFuncType func> class FolderFactory {
 188 public:
 189   static constexpr HostRuntimeFunction Create(const std::string_view &name) {
 190     return HostRuntimeFunction{name, FuncTypeAnalyzer<HostFuncType>::result,
 191         FuncTypeAnalyzer<HostFuncType>::arguments, &Fold};
 192   }
 193
 194 private:
 195   static Expr<SomeType> Fold(
 196       FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
 197     return ApplyHostFunction(func, context, std::move(args));
 198   }
 199 };
 200
 201 // Define host runtime libraries that can be used for folding and
 202 // fill their description if they are available.
 203 enum class LibraryVersion {
 204   Libm,
 205   LibmExtensions,
 206   PgmathFast,
 207   PgmathRelaxed,
 208   PgmathPrecise
 209 };
 210 template <typename HostT, LibraryVersion> struct HostRuntimeLibrary {
 211   // When specialized, this class holds a static constexpr table containing
 212   // all the HostRuntimeLibrary for functions of library LibraryVersion
 213   // that returns a value of type HostT.
 214 };
 215
 216 using HostRuntimeMap = common::StaticMultimapView<HostRuntimeFunction>;
 217
 218 // Map numerical intrinsic to  <cmath>/<complex> functions
 219 // (Note: ABS() is folded in fold-real.cpp.)
 220 template <typename HostT>
 221 struct HostRuntimeLibrary<HostT, LibraryVersion::Libm> {
 222   using F = FuncPointer<HostT, HostT>;
 223   using F2 = FuncPointer<HostT, HostT, HostT>;
 224   static constexpr HostRuntimeFunction table[]{
 225       FolderFactory<F, F{std::acos}>::Create("acos"),
 226       FolderFactory<F, F{std::acosh}>::Create("acosh"),
 227       FolderFactory<F, F{std::asin}>::Create("asin"),
 228       FolderFactory<F, F{std::asinh}>::Create("asinh"),
 229       FolderFactory<F, F{std::atan}>::Create("atan"),
 230       FolderFactory<F2, F2{std::atan2}>::Create("atan2"),
 231       FolderFactory<F, F{std::atanh}>::Create("atanh"),
 232       FolderFactory<F, F{std::cos}>::Create("cos"),
 233       FolderFactory<F, F{std::cosh}>::Create("cosh"),
 234       FolderFactory<F, F{std::erf}>::Create("erf"),
 235       FolderFactory<F, F{std::erfc}>::Create("erfc"),
 236       FolderFactory<F, F{common::ErfcScaled}>::Create("erfc_scaled"),
 237       FolderFactory<F, F{std::exp}>::Create("exp"),
 238       FolderFactory<F, F{std::tgamma}>::Create("gamma"),
 239       FolderFactory<F, F{std::log}>::Create("log"),
 240       FolderFactory<F, F{std::log10}>::Create("log10"),
 241       FolderFactory<F, F{std::lgamma}>::Create("log_gamma"),
 242       FolderFactory<F2, F2{std::pow}>::Create("pow"),
 243       FolderFactory<F, F{std::sin}>::Create("sin"),
 244       FolderFactory<F, F{std::sinh}>::Create("sinh"),
 245       FolderFactory<F, F{std::tan}>::Create("tan"),
 246       FolderFactory<F, F{std::tanh}>::Create("tanh"),
 247   };
 248   // Note: cmath does not have modulo and erfc_scaled equivalent
 249
 250   // Note regarding  lack of bessel function support:
 251   // C++17 defined standard Bessel math functions std::cyl_bessel_j
 252   // and std::cyl_neumann that can be used for Fortran j and y
 253   // bessel functions. However, they are not yet implemented in
 254   // clang libc++ (ok in GNU libstdc++). C maths functions j0...
 255   // are not C standard but a GNU extension so they are not used
 256   // to avoid introducing incompatibilities.
 257   // Use libpgmath to get bessel function folding support.
 258   // TODO:  Add Bessel functions when possible.
 259   static constexpr HostRuntimeMap map{table};
 260   static_assert(map.Verify(), "map must be sorted");
 261 };
 262
 263 // Helpers to map complex std::pow whose resolution in F2{std::pow} is
 264 // ambiguous as of clang++ 20.
 265 template <typename HostT>
 266 static std::complex<HostT> StdPowF2(
 267     const std::complex<HostT> &x, const std::complex<HostT> &y) {
 268   return std::pow(x, y);
 269 }
 270 template <typename HostT>
 271 static std::complex<HostT> StdPowF2A(
 272     const HostT &x, const std::complex<HostT> &y) {
 273   return std::pow(x, y);
 274 }
 275 template <typename HostT>
 276 static std::complex<HostT> StdPowF2B(
 277     const std::complex<HostT> &x, const HostT &y) {
 278   return std::pow(x, y);
 279 }
 280
 281 #ifdef _AIX
 282 #ifdef __clang_major__
 283 #pragma clang diagnostic ignored "-Wc99-extensions"
 284 #endif
 285
 286 extern "C" {
 287 float _Complex cacosf(float _Complex);
 288 double _Complex cacos(double _Complex);
 289 float _Complex csqrtf(float _Complex);
 290 double _Complex csqrt(double _Complex);
 291 }
 292
 293 enum CRI { Real, Imag };
 294 template <typename TR, typename TA> static TR &reIm(TA &x, CRI n) {
 295   return reinterpret_cast<TR(&)[2]>(x)[n];
 296 }
 297 template <typename TR, typename T> static TR CppToC(const std::complex<T> &x) {
 298   TR r;
 299   reIm<T, TR>(r, CRI::Real) = x.real();
 300   reIm<T, TR>(r, CRI::Imag) = x.imag();
 301   return r;
 302 }
 303 template <typename T, typename TA> static std::complex<T> CToCpp(const TA &x) {
 304   TA &z{const_cast<TA &>(x)};
 305   return std::complex<T>(reIm<T, TA>(z, CRI::Real), reIm<T, TA>(z, CRI::Imag));
 306 }
 307 #endif
 308
 309 template <typename HostT>
 310 static std::complex<HostT> CSqrt(const std::complex<HostT> &x) {
 311   std::complex<HostT> res;
 312 #ifdef _AIX
 313   // On AIX, the implementation of csqrt[f] and std::sqrt is different,
 314   // use csqrt[f] in folding.
 315   if constexpr (std::is_same_v<HostT, float>) {
 316     float _Complex r{csqrtf(CppToC<float _Complex, float>(x))};
 317     res = CToCpp<float, float _Complex>(r);
 318   } else if constexpr (std::is_same_v<HostT, double>) {
 319     double _Complex r{csqrt(CppToC<double _Complex, double>(x))};
 320     res = CToCpp<double, double _Complex>(r);
 321   } else {
 322     DIE("bad complex component type");
 323   }
 324 #else
 325   res = std::sqrt(x);
 326 #endif
 327   return res;
 328 }
 329
 330 template <typename HostT>
 331 static std::complex<HostT> CAcos(const std::complex<HostT> &x) {
 332   std::complex<HostT> res;
 333 #ifdef _AIX
 334   // On AIX, the implementation of cacos[f] and std::acos is different,
 335   // use cacos[f] in folding.
 336   if constexpr (std::is_same_v<HostT, float>) {
 337     float _Complex r{cacosf(CppToC<float _Complex, float>(x))};
 338     res = CToCpp<float, float _Complex>(r);
 339   } else if constexpr (std::is_same_v<HostT, double>) {
 340     double _Complex r{cacos(CppToC<double _Complex, double>(x))};
 341     res = CToCpp<double, double _Complex>(r);
 342   } else {
 343     DIE("bad complex component type");
 344   }
 345 #else
 346   res = std::acos(x);
 347 #endif
 348   return res;
 349 }
 350
 351 template <typename HostT>
 352 struct HostRuntimeLibrary<std::complex<HostT>, LibraryVersion::Libm> {
 353   using F = FuncPointer<std::complex<HostT>, const std::complex<HostT> &>;
 354   using F2 = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
 355       const std::complex<HostT> &>;
 356   using F2A = FuncPointer<std::complex<HostT>, const HostT &,
 357       const std::complex<HostT> &>;
 358   using F2B = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
 359       const HostT &>;
 360   static constexpr HostRuntimeFunction table[]{
 361       FolderFactory<F, F{CAcos}>::Create("acos"),
 362       FolderFactory<F, F{std::acosh}>::Create("acosh"),
 363       FolderFactory<F, F{std::asin}>::Create("asin"),
 364       FolderFactory<F, F{std::asinh}>::Create("asinh"),
 365       FolderFactory<F, F{std::atan}>::Create("atan"),
 366       FolderFactory<F, F{std::atanh}>::Create("atanh"),
 367       FolderFactory<F, F{std::cos}>::Create("cos"),
 368       FolderFactory<F, F{std::cosh}>::Create("cosh"),
 369       FolderFactory<F, F{std::exp}>::Create("exp"),
 370       FolderFactory<F, F{std::log}>::Create("log"),
 371       FolderFactory<F2, F2{StdPowF2}>::Create("pow"),
 372       FolderFactory<F2A, F2A{StdPowF2A}>::Create("pow"),
 373       FolderFactory<F2B, F2B{StdPowF2B}>::Create("pow"),
 374       FolderFactory<F, F{std::sin}>::Create("sin"),
 375       FolderFactory<F, F{std::sinh}>::Create("sinh"),
 376       FolderFactory<F, F{CSqrt}>::Create("sqrt"),
 377       FolderFactory<F, F{std::tan}>::Create("tan"),
 378       FolderFactory<F, F{std::tanh}>::Create("tanh"),
 379   };
 380   static constexpr HostRuntimeMap map{table};
 381   static_assert(map.Verify(), "map must be sorted");
 382 };
 383 // Note regarding cmath:
 384 //  - cmath does not have modulo and erfc_scaled equivalent
 385 //  - C++17 defined standard Bessel math functions std::cyl_bessel_j
 386 //    and std::cyl_neumann that can be used for Fortran j and y
 387 //    bessel functions. However, they are not yet implemented in
 388 //    clang libc++ (ok in GNU libstdc++). Instead, the Posix libm
 389 //    extensions are used when available below.
 390
 391 #if _POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600
 392 /// Define libm extensions
 393 /// Bessel functions are defined in POSIX.1-2001.
 394
 395 // Remove float bessel functions for AIX and Darwin as they are not supported
 396 #if !defined(_AIX) && !defined(__APPLE__)
 397 template <> struct HostRuntimeLibrary<float, LibraryVersion::LibmExtensions> {
 398   using F = FuncPointer<float, float>;
 399   using FN = FuncPointer<float, int, float>;
 400   static constexpr HostRuntimeFunction table[]{
 401       FolderFactory<F, F{::j0f}>::Create("bessel_j0"),
 402       FolderFactory<F, F{::j1f}>::Create("bessel_j1"),
 403       FolderFactory<FN, FN{::jnf}>::Create("bessel_jn"),
 404       FolderFactory<F, F{::y0f}>::Create("bessel_y0"),
 405       FolderFactory<F, F{::y1f}>::Create("bessel_y1"),
 406       FolderFactory<FN, FN{::ynf}>::Create("bessel_yn"),
 407   };
 408   static constexpr HostRuntimeMap map{table};
 409   static_assert(map.Verify(), "map must be sorted");
 410 };
 411 #endif
 412
 413 #if HAS_QUADMATHLIB
 414 template <> struct HostRuntimeLibrary<__float128, LibraryVersion::Libm> {
 415   using F = FuncPointer<__float128, __float128>;
 416   using F2 = FuncPointer<__float128, __float128, __float128>;
 417   using FN = FuncPointer<__float128, int, __float128>;
 418   static constexpr HostRuntimeFunction table[]{
 419       FolderFactory<F, F{::acosq}>::Create("acos"),
 420       FolderFactory<F, F{::acoshq}>::Create("acosh"),
 421       FolderFactory<F, F{::asinq}>::Create("asin"),
 422       FolderFactory<F, F{::asinhq}>::Create("asinh"),
 423       FolderFactory<F, F{::atanq}>::Create("atan"),
 424       FolderFactory<F2, F2{::atan2q}>::Create("atan2"),
 425       FolderFactory<F, F{::atanhq}>::Create("atanh"),
 426       FolderFactory<F, F{::j0q}>::Create("bessel_j0"),
 427       FolderFactory<F, F{::j1q}>::Create("bessel_j1"),
 428       FolderFactory<FN, FN{::jnq}>::Create("bessel_jn"),
 429       FolderFactory<F, F{::y0q}>::Create("bessel_y0"),
 430       FolderFactory<F, F{::y1q}>::Create("bessel_y1"),
 431       FolderFactory<FN, FN{::ynq}>::Create("bessel_yn"),
 432       FolderFactory<F, F{::cosq}>::Create("cos"),
 433       FolderFactory<F, F{::coshq}>::Create("cosh"),
 434       FolderFactory<F, F{::erfq}>::Create("erf"),
 435       FolderFactory<F, F{::erfcq}>::Create("erfc"),
 436       FolderFactory<F, F{::expq}>::Create("exp"),
 437       FolderFactory<F, F{::tgammaq}>::Create("gamma"),
 438       FolderFactory<F, F{::logq}>::Create("log"),
 439       FolderFactory<F, F{::log10q}>::Create("log10"),
 440       FolderFactory<F, F{::lgammaq}>::Create("log_gamma"),
 441       FolderFactory<F2, F2{::powq}>::Create("pow"),
 442       FolderFactory<F, F{::sinq}>::Create("sin"),
 443       FolderFactory<F, F{::sinhq}>::Create("sinh"),
 444       FolderFactory<F, F{::tanq}>::Create("tan"),
 445       FolderFactory<F, F{::tanhq}>::Create("tanh"),
 446   };
 447   static constexpr HostRuntimeMap map{table};
 448   static_assert(map.Verify(), "map must be sorted");
 449 };
 450 template <> struct HostRuntimeLibrary<__complex128, LibraryVersion::Libm> {
 451   using F = FuncPointer<__complex128, __complex128>;
 452   using F2 = FuncPointer<__complex128, __complex128, __complex128>;
 453   static constexpr HostRuntimeFunction table[]{
 454       FolderFactory<F, F{::cacosq}>::Create("acos"),
 455       FolderFactory<F, F{::cacoshq}>::Create("acosh"),
 456       FolderFactory<F, F{::casinq}>::Create("asin"),
 457       FolderFactory<F, F{::casinhq}>::Create("asinh"),
 458       FolderFactory<F, F{::catanq}>::Create("atan"),
 459       FolderFactory<F, F{::catanhq}>::Create("atanh"),
 460       FolderFactory<F, F{::ccosq}>::Create("cos"),
 461       FolderFactory<F, F{::ccoshq}>::Create("cosh"),
 462       FolderFactory<F, F{::cexpq}>::Create("exp"),
 463       FolderFactory<F, F{::clogq}>::Create("log"),
 464       FolderFactory<F2, F2{::cpowq}>::Create("pow"),
 465       FolderFactory<F, F{::csinq}>::Create("sin"),
 466       FolderFactory<F, F{::csinhq}>::Create("sinh"),
 467       FolderFactory<F, F{::csqrtq}>::Create("sqrt"),
 468       FolderFactory<F, F{::ctanq}>::Create("tan"),
 469       FolderFactory<F, F{::ctanhq}>::Create("tanh"),
 470   };
 471   static constexpr HostRuntimeMap map{table};
 472   static_assert(map.Verify(), "map must be sorted");
 473 };
 474 #endif
 475
 476 template <> struct HostRuntimeLibrary<double, LibraryVersion::LibmExtensions> {
 477   using F = FuncPointer<double, double>;
 478   using FN = FuncPointer<double, int, double>;
 479   static constexpr HostRuntimeFunction table[]{
 480       FolderFactory<F, F{::j0}>::Create("bessel_j0"),
 481       FolderFactory<F, F{::j1}>::Create("bessel_j1"),
 482       FolderFactory<FN, FN{::jn}>::Create("bessel_jn"),
 483       FolderFactory<F, F{::y0}>::Create("bessel_y0"),
 484       FolderFactory<F, F{::y1}>::Create("bessel_y1"),
 485       FolderFactory<FN, FN{::yn}>::Create("bessel_yn"),
 486   };
 487   static constexpr HostRuntimeMap map{table};
 488   static_assert(map.Verify(), "map must be sorted");
 489 };
 490
 491 #if defined(__GLIBC__) && (HAS_FLOAT80 || HAS_LDBL128)
 492 template <>
 493 struct HostRuntimeLibrary<long double, LibraryVersion::LibmExtensions> {
 494   using F = FuncPointer<long double, long double>;
 495   using FN = FuncPointer<long double, int, long double>;
 496   static constexpr HostRuntimeFunction table[]{
 497       FolderFactory<F, F{::j0l}>::Create("bessel_j0"),
 498       FolderFactory<F, F{::j1l}>::Create("bessel_j1"),
 499       FolderFactory<FN, FN{::jnl}>::Create("bessel_jn"),
 500       FolderFactory<F, F{::y0l}>::Create("bessel_y0"),
 501       FolderFactory<F, F{::y1l}>::Create("bessel_y1"),
 502       FolderFactory<FN, FN{::ynl}>::Create("bessel_yn"),
 503   };
 504   static constexpr HostRuntimeMap map{table};
 505   static_assert(map.Verify(), "map must be sorted");
 506 };
 507 #endif // HAS_FLOAT80 || HAS_LDBL128
 508 #endif //_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600
 509
 510 /// Define pgmath description
 511 #if LINK_WITH_LIBPGMATH
 512 // Only use libpgmath for folding if it is available.
 513 // First declare all libpgmaths functions
 514 #define PGMATH_LINKING
 515 #define PGMATH_DECLARE
 516 #include "flang/Evaluate/pgmath.h.inc"
 517
 518 #define REAL_FOLDER(name, func) \
 519   FolderFactory<decltype(&func), &func>::Create(#name)
 520 template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathFast> {
 521   static constexpr HostRuntimeFunction table[]{
 522 #define PGMATH_FAST
 523 #define PGMATH_USE_S(name, func) REAL_FOLDER(name, func),
 524 #include "flang/Evaluate/pgmath.h.inc"
 525   };
 526   static constexpr HostRuntimeMap map{table};
 527   static_assert(map.Verify(), "map must be sorted");
 528 };
 529 template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathFast> {
 530   static constexpr HostRuntimeFunction table[]{
 531 #define PGMATH_FAST
 532 #define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
 533 #include "flang/Evaluate/pgmath.h.inc"
 534   };
 535   static constexpr HostRuntimeMap map{table};
 536   static_assert(map.Verify(), "map must be sorted");
 537 };
 538 template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathRelaxed> {
 539   static constexpr HostRuntimeFunction table[]{
 540 #define PGMATH_RELAXED
 541 #define PGMATH_USE_S(name, func) REAL_FOLDER(name, func),
 542 #include "flang/Evaluate/pgmath.h.inc"
 543   };
 544   static constexpr HostRuntimeMap map{table};
 545   static_assert(map.Verify(), "map must be sorted");
 546 };
 547 template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathRelaxed> {
 548   static constexpr HostRuntimeFunction table[]{
 549 #define PGMATH_RELAXED
 550 #define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
 551 #include "flang/Evaluate/pgmath.h.inc"
 552   };
 553   static constexpr HostRuntimeMap map{table};
 554   static_assert(map.Verify(), "map must be sorted");
 555 };
 556 template <> struct HostRuntimeLibrary<float, LibraryVersion::PgmathPrecise> {
 557   static constexpr HostRuntimeFunction table[]{
 558 #define PGMATH_PRECISE
 559 #define PGMATH_USE_S(name, func) REAL_FOLDER(name, func),
 560 #include "flang/Evaluate/pgmath.h.inc"
 561   };
 562   static constexpr HostRuntimeMap map{table};
 563   static_assert(map.Verify(), "map must be sorted");
 564 };
 565 template <> struct HostRuntimeLibrary<double, LibraryVersion::PgmathPrecise> {
 566   static constexpr HostRuntimeFunction table[]{
 567 #define PGMATH_PRECISE
 568 #define PGMATH_USE_D(name, func) REAL_FOLDER(name, func),
 569 #include "flang/Evaluate/pgmath.h.inc"
 570   };
 571   static constexpr HostRuntimeMap map{table};
 572   static_assert(map.Verify(), "map must be sorted");
 573 };
 574
 575 // TODO: double _Complex/float _Complex have been removed from llvm flang
 576 // pgmath.h.inc because they caused warnings, they need to be added back
 577 // so that the complex pgmath versions can be used when requested.
 578
 579 #endif /* LINK_WITH_LIBPGMATH */
 580
 581 // Helper to check if a HostRuntimeLibrary specialization exists
 582 template <typename T, typename = void> struct IsAvailable : std::false_type {};
 583 template <typename T>
 584 struct IsAvailable<T, decltype((void)T::table, void())> : std::true_type {};
 585 // Define helpers to find host runtime library map according to desired version
 586 // and type.
 587 template <typename HostT, LibraryVersion version>
 588 static const HostRuntimeMap *GetHostRuntimeMapHelper(
 589     [[maybe_unused]] DynamicType resultType) {
 590   // A library must only be instantiated if LibraryVersion is
 591   // available on the host and if HostT maps to a Fortran type.
 592   // For instance, whenever long double and double are both 64-bits, double
 593   // is mapped to Fortran 64bits real type, and long double will be left
 594   // unmapped.
 595   if constexpr (host::FortranTypeExists<HostT>()) {
 596     using Lib = HostRuntimeLibrary<HostT, version>;
 597     if constexpr (IsAvailable<Lib>::value) {
 598       if (host::FortranType<HostT>{}.GetType() == resultType) {
 599         return &Lib::map;
 600       }
 601     }
 602   }
 603   return nullptr;
 604 }
 605 template <LibraryVersion version>
 606 static const HostRuntimeMap *GetHostRuntimeMapVersion(DynamicType resultType) {
 607   if (resultType.category() == TypeCategory::Real) {
 608     if (const auto *map{GetHostRuntimeMapHelper<float, version>(resultType)}) {
 609       return map;
 610     }
 611     if (const auto *map{GetHostRuntimeMapHelper<double, version>(resultType)}) {
 612       return map;
 613     }
 614     if (const auto *map{
 615             GetHostRuntimeMapHelper<long double, version>(resultType)}) {
 616       return map;
 617     }
 618 #if HAS_QUADMATHLIB
 619     if (const auto *map{
 620             GetHostRuntimeMapHelper<__float128, version>(resultType)}) {
 621       return map;
 622     }
 623 #endif
 624   }
 625   if (resultType.category() == TypeCategory::Complex) {
 626     if (const auto *map{GetHostRuntimeMapHelper<std::complex<float>, version>(
 627             resultType)}) {
 628       return map;
 629     }
 630     if (const auto *map{GetHostRuntimeMapHelper<std::complex<double>, version>(
 631             resultType)}) {
 632       return map;
 633     }
 634     if (const auto *map{
 635             GetHostRuntimeMapHelper<std::complex<long double>, version>(
 636                 resultType)}) {
 637       return map;
 638     }
 639 #if HAS_QUADMATHLIB
 640     if (const auto *map{
 641             GetHostRuntimeMapHelper<__complex128, version>(resultType)}) {
 642       return map;
 643     }
 644 #endif
 645   }
 646   return nullptr;
 647 }
 648 static const HostRuntimeMap *GetHostRuntimeMap(
 649     LibraryVersion version, DynamicType resultType) {
 650   switch (version) {
 651   case LibraryVersion::Libm:
 652     return GetHostRuntimeMapVersion<LibraryVersion::Libm>(resultType);
 653   case LibraryVersion::LibmExtensions:
 654     return GetHostRuntimeMapVersion<LibraryVersion::LibmExtensions>(resultType);
 655   case LibraryVersion::PgmathPrecise:
 656     return GetHostRuntimeMapVersion<LibraryVersion::PgmathPrecise>(resultType);
 657   case LibraryVersion::PgmathRelaxed:
 658     return GetHostRuntimeMapVersion<LibraryVersion::PgmathRelaxed>(resultType);
 659   case LibraryVersion::PgmathFast:
 660     return GetHostRuntimeMapVersion<LibraryVersion::PgmathFast>(resultType);
 661   }
 662   return nullptr;
 663 }
 664
 665 static const HostRuntimeFunction *SearchInHostRuntimeMap(
 666     const HostRuntimeMap &map, const std::string &name, DynamicType resultType,
 667     const std::vector<DynamicType> &argTypes) {
 668   auto sameNameRange{map.equal_range(name)};
 669   for (const auto *iter{sameNameRange.first}; iter != sameNameRange.second;
 670        ++iter) {
 671     if (iter->resultType == resultType && iter->argumentTypes == argTypes) {
 672       return &*iter;
 673     }
 674   }
 675   return nullptr;
 676 }
 677
 678 // Search host runtime libraries for an exact type match.
 679 static const HostRuntimeFunction *SearchHostRuntime(const std::string &name,
 680     DynamicType resultType, const std::vector<DynamicType> &argTypes) {
 681   // TODO: When command line options regarding targeted numerical library is
 682   // available, this needs to be revisited to take it into account. So far,
 683   // default to libpgmath if F18 is built with it.
 684 #if LINK_WITH_LIBPGMATH
 685   if (const auto *map{
 686           GetHostRuntimeMap(LibraryVersion::PgmathPrecise, resultType)}) {
 687     if (const auto *hostFunction{
 688             SearchInHostRuntimeMap(*map, name, resultType, argTypes)}) {
 689       return hostFunction;
 690     }
 691   }
 692   // Default to libm if functions or types are not available in pgmath.
 693 #endif
 694   if (const auto *map{GetHostRuntimeMap(LibraryVersion::Libm, resultType)}) {
 695     if (const auto *hostFunction{
 696             SearchInHostRuntimeMap(*map, name, resultType, argTypes)}) {
 697       return hostFunction;
 698     }
 699   }
 700   if (const auto *map{
 701           GetHostRuntimeMap(LibraryVersion::LibmExtensions, resultType)}) {
 702     if (const auto *hostFunction{
 703             SearchInHostRuntimeMap(*map, name, resultType, argTypes)}) {
 704       return hostFunction;
 705     }
 706   }
 707   return nullptr;
 708 }
 709
 710 // Return a DynamicType that can hold all values of a given type.
 711 // This is used to allow 16bit float to be folded with 32bits and
 712 // x87 float to be folded with IEEE 128bits.
 713 static DynamicType BiggerType(DynamicType type) {
 714   if (type.category() == TypeCategory::Real ||
 715       type.category() == TypeCategory::Complex) {
 716     // 16 bits floats to IEEE 32 bits float
 717     if (type.kind() == common::RealKindForPrecision(11) ||
 718         type.kind() == common::RealKindForPrecision(8)) {
 719       return {type.category(), common::RealKindForPrecision(24)};
 720     }
 721     // x87 float to IEEE 128 bits float
 722     if (type.kind() == common::RealKindForPrecision(64)) {
 723       return {type.category(), common::RealKindForPrecision(113)};
 724     }
 725   }
 726   return type;
 727 }
 728
 729 /// Structure to register intrinsic argument checks that must be performed.
 730 using ArgumentVerifierFunc = bool (*)(
 731     const std::vector<Expr<SomeType>> &, FoldingContext &);
 732 struct ArgumentVerifier {
 733   using Key = std::string_view;
 734   // Needed for implicit compare with keys.
 735   constexpr operator Key() const { return key; }
 736   Key key;
 737   ArgumentVerifierFunc verifier;
 738 };
 739
 740 static constexpr int lastArg{-1};
 741 static constexpr int firstArg{0};
 742
 743 static const Expr<SomeType> &GetArg(
 744     int position, const std::vector<Expr<SomeType>> &args) {
 745   if (position == lastArg) {
 746     CHECK(!args.empty());
 747     return args.back();
 748   }
 749   CHECK(position >= 0 && static_cast<std::size_t>(position) < args.size());
 750   return args[position];
 751 }
 752
 753 template <typename T>
 754 static bool IsInRange(const Expr<T> &expr, int lb, int ub) {
 755   if (auto scalar{GetScalarConstantValue<T>(expr)}) {
 756     auto lbValue{Scalar<T>::FromInteger(value::Integer<8>{lb}).value};
 757     auto ubValue{Scalar<T>::FromInteger(value::Integer<8>{ub}).value};
 758     return Satisfies(RelationalOperator::LE, lbValue.Compare(*scalar)) &&
 759         Satisfies(RelationalOperator::LE, scalar->Compare(ubValue));
 760   }
 761   return true;
 762 }
 763
 764 /// Verify that the argument in an intrinsic call belongs to [lb, ub] if is
 765 /// real.
 766 template <int lb, int ub>
 767 static bool VerifyInRangeIfReal(
 768     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 769   if (const auto *someReal{
 770           std::get_if<Expr<SomeReal>>(&GetArg(firstArg, args).u)}) {
 771     bool isInRange{
 772         std::visit([&](const auto &x) -> bool { return IsInRange(x, lb, ub); },
 773             someReal->u)};
 774     if (!isInRange) {
 775       context.messages().Say(
 776           "argument is out of range [%d., %d.]"_warn_en_US, lb, ub);
 777     }
 778     return isInRange;
 779   }
 780   return true;
 781 }
 782
 783 template <int argPosition, const char *argName>
 784 static bool VerifyStrictlyPositiveIfReal(
 785     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 786   if (const auto *someReal =
 787           std::get_if<Expr<SomeReal>>(&GetArg(argPosition, args).u)) {
 788     const bool isStrictlyPositive{std::visit(
 789         [&](const auto &x) -> bool {
 790           using T = typename std::decay_t<decltype(x)>::Result;
 791           auto scalar{GetScalarConstantValue<T>(x)};
 792           return Satisfies(
 793               RelationalOperator::LT, Scalar<T>{}.Compare(*scalar));
 794         },
 795         someReal->u)};
 796     if (!isStrictlyPositive) {
 797       context.messages().Say(
 798           "argument '%s' must be strictly positive"_warn_en_US, argName);
 799     }
 800     return isStrictlyPositive;
 801   }
 802   return true;
 803 }
 804
 805 /// Verify that an intrinsic call argument is not zero if it is real.
 806 template <int argPosition, const char *argName>
 807 static bool VerifyNotZeroIfReal(
 808     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 809   if (const auto *someReal =
 810           std::get_if<Expr<SomeReal>>(&GetArg(argPosition, args).u)) {
 811     const bool isNotZero{std::visit(
 812         [&](const auto &x) -> bool {
 813           using T = typename std::decay_t<decltype(x)>::Result;
 814           auto scalar{GetScalarConstantValue<T>(x)};
 815           return !scalar || !scalar->IsZero();
 816         },
 817         someReal->u)};
 818     if (!isNotZero) {
 819       context.messages().Say(
 820           "argument '%s' must be different from zero"_warn_en_US, argName);
 821     }
 822     return isNotZero;
 823   }
 824   return true;
 825 }
 826
 827 /// Verify that the argument in an intrinsic call is not zero if is complex.
 828 static bool VerifyNotZeroIfComplex(
 829     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 830   if (const auto *someComplex =
 831           std::get_if<Expr<SomeComplex>>(&GetArg(firstArg, args).u)) {
 832     const bool isNotZero{std::visit(
 833         [&](const auto &z) -> bool {
 834           using T = typename std::decay_t<decltype(z)>::Result;
 835           auto scalar{GetScalarConstantValue<T>(z)};
 836           return !scalar || !scalar->IsZero();
 837         },
 838         someComplex->u)};
 839     if (!isNotZero) {
 840       context.messages().Say(
 841           "complex argument must be different from zero"_warn_en_US);
 842     }
 843     return isNotZero;
 844   }
 845   return true;
 846 }
 847
 848 // Verify that the argument in an intrinsic call is not zero and not a negative
 849 // integer.
 850 static bool VerifyGammaLikeArgument(
 851     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 852   if (const auto *someReal =
 853           std::get_if<Expr<SomeReal>>(&GetArg(firstArg, args).u)) {
 854     const bool isValid{std::visit(
 855         [&](const auto &x) -> bool {
 856           using T = typename std::decay_t<decltype(x)>::Result;
 857           auto scalar{GetScalarConstantValue<T>(x)};
 858           if (scalar) {
 859             return !scalar->IsZero() &&
 860                 !(scalar->IsNegative() &&
 861                     scalar->ToWholeNumber().value == scalar);
 862           }
 863           return true;
 864         },
 865         someReal->u)};
 866     if (!isValid) {
 867       context.messages().Say(
 868           "argument must not be a negative integer or zero"_warn_en_US);
 869     }
 870     return isValid;
 871   }
 872   return true;
 873 }
 874
 875 // Verify that two real arguments are not both zero.
 876 static bool VerifyAtan2LikeArguments(
 877     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 878   if (const auto *someReal =
 879           std::get_if<Expr<SomeReal>>(&GetArg(firstArg, args).u)) {
 880     const bool isValid{std::visit(
 881         [&](const auto &typedExpr) -> bool {
 882           using T = typename std::decay_t<decltype(typedExpr)>::Result;
 883           auto x{GetScalarConstantValue<T>(typedExpr)};
 884           auto y{GetScalarConstantValue<T>(GetArg(lastArg, args))};
 885           if (x && y) {
 886             return !(x->IsZero() && y->IsZero());
 887           }
 888           return true;
 889         },
 890         someReal->u)};
 891     if (!isValid) {
 892       context.messages().Say(
 893           "'x' and 'y' arguments must not be both zero"_warn_en_US);
 894     }
 895     return isValid;
 896   }
 897   return true;
 898 }
 899
 900 template <ArgumentVerifierFunc... F>
 901 static bool CombineVerifiers(
 902     const std::vector<Expr<SomeType>> &args, FoldingContext &context) {
 903   return (... && F(args, context));
 904 }
 905
 906 /// Define argument names to be used error messages when the intrinsic have
 907 /// several arguments.
 908 static constexpr char xName[]{"x"};
 909 static constexpr char pName[]{"p"};
 910
 911 /// Register argument verifiers for all intrinsics folded with runtime.
 912 static constexpr ArgumentVerifier intrinsicArgumentVerifiers[]{
 913     {"acos", VerifyInRangeIfReal<-1, 1>},
 914     {"asin", VerifyInRangeIfReal<-1, 1>},
 915     {"atan2", VerifyAtan2LikeArguments},
 916     {"bessel_y0", VerifyStrictlyPositiveIfReal<firstArg, xName>},
 917     {"bessel_y1", VerifyStrictlyPositiveIfReal<firstArg, xName>},
 918     {"bessel_yn", VerifyStrictlyPositiveIfReal<lastArg, xName>},
 919     {"gamma", VerifyGammaLikeArgument},
 920     {"log",
 921         CombineVerifiers<VerifyStrictlyPositiveIfReal<firstArg, xName>,
 922             VerifyNotZeroIfComplex>},
 923     {"log10", VerifyStrictlyPositiveIfReal<firstArg, xName>},
 924     {"log_gamma", VerifyGammaLikeArgument},
 925     {"mod", VerifyNotZeroIfReal<lastArg, pName>},
 926 };
 927
 928 const ArgumentVerifierFunc *findVerifier(const std::string &intrinsicName) {
 929   static constexpr Fortran::common::StaticMultimapView<ArgumentVerifier>
 930       verifiers(intrinsicArgumentVerifiers);
 931   static_assert(verifiers.Verify(), "map must be sorted");
 932   auto range{verifiers.equal_range(intrinsicName)};
 933   if (range.first != range.second) {
 934     return &range.first->verifier;
 935   }
 936   return nullptr;
 937 }
 938
 939 /// Ensure argument verifiers, if any, are run before calling the runtime
 940 /// wrapper to fold an intrinsic.
 941 static HostRuntimeWrapper AddArgumentVerifierIfAny(
 942     const std::string &intrinsicName, const HostRuntimeFunction &hostFunction) {
 943   if (const auto *verifier{findVerifier(intrinsicName)}) {
 944     const HostRuntimeFunction *hostFunctionPtr = &hostFunction;
 945     return [hostFunctionPtr, verifier](
 946                FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
 947       const bool validArguments{(*verifier)(args, context)};
 948       if (!validArguments) {
 949         // Silence fp signal warnings since a more detailed warning about
 950         // invalid arguments was already emitted.
 951         parser::Messages localBuffer;
 952         parser::ContextualMessages localMessages{&localBuffer};
 953         FoldingContext localContext{context, localMessages};
 954         return hostFunctionPtr->folder(localContext, std::move(args));
 955       }
 956       return hostFunctionPtr->folder(context, std::move(args));
 957     };
 958   }
 959   return hostFunction.folder;
 960 }
 961
 962 std::optional<HostRuntimeWrapper> GetHostRuntimeWrapper(const std::string &name,
 963     DynamicType resultType, const std::vector<DynamicType> &argTypes) {
 964   if (const auto *hostFunction{SearchHostRuntime(name, resultType, argTypes)}) {
 965     return AddArgumentVerifierIfAny(name, *hostFunction);
 966   }
 967   // If no exact match, search with "bigger" types and insert type
 968   // conversions around the folder.
 969   std::vector<evaluate::DynamicType> biggerArgTypes;
 970   evaluate::DynamicType biggerResultType{BiggerType(resultType)};
 971   for (auto type : argTypes) {
 972     biggerArgTypes.emplace_back(BiggerType(type));
 973   }
 974   if (const auto *hostFunction{
 975           SearchHostRuntime(name, biggerResultType, biggerArgTypes)}) {
 976     auto hostFolderWithChecks{AddArgumentVerifierIfAny(name, *hostFunction)};
 977     return [hostFunction, resultType, hostFolderWithChecks](
 978                FoldingContext &context, std::vector<Expr<SomeType>> &&args) {
 979       auto nArgs{args.size()};
 980       for (size_t i{0}; i < nArgs; ++i) {
 981         args[i] = Fold(context,
 982             ConvertToType(hostFunction->argumentTypes[i], std::move(args[i]))
 983                 .value());
 984       }
 985       return Fold(context,
 986           ConvertToType(
 987               resultType, hostFolderWithChecks(context, std::move(args)))
 988               .value());
 989     };
 990   }
 991   return std::nullopt;
 992 }
 993 } // namespace Fortran::evaluate