xref: /llvm-project/flang/lib/Evaluate/intrinsics-library.cpp (revision 1f8790050b0e99e7b46cc69518aa84f46f50738e) 
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 HostIntrinsicProceduresLibrary is built with <cmath> and
12 // <complex> functions that are guaranteed to exist from the C++ standard.
13 
14 #include "intrinsics-library-templates.h"
15 #include <cmath>
16 #include <complex>
17 
18 namespace Fortran::evaluate {
19 
20 // Note: argument passing is ignored in equivalence
21 bool HostIntrinsicProceduresLibrary::HasEquivalentProcedure(
22     const IntrinsicProcedureRuntimeDescription &sym) const {
23   const auto rteProcRange{procedures_.equal_range(sym.name)};
24   const size_t nargs{sym.argumentsType.size()};
25   for (auto iter{rteProcRange.first}; iter != rteProcRange.second; ++iter) {
26     if (nargs == iter->second.argumentsType.size() &&
27         sym.returnType == iter->second.returnType &&
28         (sym.isElemental || iter->second.isElemental)) {
29       bool match{true};
30       int pos{0};
31       for (const auto &type : sym.argumentsType) {
32         if (type != iter->second.argumentsType[pos++]) {
33           match = false;
34           break;
35         }
36       }
37       if (match) {
38         return true;
39       }
40     }
41   }
42   return false;
43 }
44 
45 // Map numerical intrinsic to  <cmath>/<complex> functions
46 
47 // Define which host runtime functions will be used for folding
48 
49 template <typename HostT>
50 static void AddLibmRealHostProcedures(
51     HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
52   using F = FuncPointer<HostT, HostT>;
53   using F2 = FuncPointer<HostT, HostT, HostT>;
54   HostRuntimeIntrinsicProcedure libmSymbols[]{
55       {"acos", F{std::acos}, true},
56       {"acosh", F{std::acosh}, true},
57       {"asin", F{std::asin}, true},
58       {"asinh", F{std::asinh}, true},
59       {"atan", F{std::atan}, true},
60       {"atan", F2{std::atan2}, true},
61       {"atanh", F{std::atanh}, true},
62       {"cos", F{std::cos}, true},
63       {"cosh", F{std::cosh}, true},
64       {"erf", F{std::erf}, true},
65       {"erfc", F{std::erfc}, true},
66       {"exp", F{std::exp}, true},
67       {"gamma", F{std::tgamma}, true},
68       {"hypot", F2{std::hypot}, true},
69       {"log", F{std::log}, true},
70       {"log10", F{std::log10}, true},
71       {"log_gamma", F{std::lgamma}, true},
72       {"mod", F2{std::fmod}, true},
73       {"pow", F2{std::pow}, true},
74       {"sin", F{std::sin}, true},
75       {"sinh", F{std::sinh}, true},
76       {"sqrt", F{std::sqrt}, true},
77       {"tan", F{std::tan}, true},
78       {"tanh", F{std::tanh}, true},
79   };
80   // Note: cmath does not have modulo and erfc_scaled equivalent
81 
82   // Note regarding  lack of bessel function support:
83   // C++17 defined standard Bessel math functions std::cyl_bessel_j
84   // and std::cyl_neumann that can be used for Fortran j and y
85   // bessel functions. However, they are not yet implemented in
86   // clang libc++ (ok in GNU libstdc++). C maths functions j0...
87   // are not C standard but a GNU extension so they are not used
88   // to avoid introducing incompatibilities.
89   // Use libpgmath to get bessel function folding support.
90   // TODO:  Add Bessel functions when possible.
91 
92   for (auto sym : libmSymbols) {
93     if (!hostIntrinsicLibrary.HasEquivalentProcedure(sym)) {
94       hostIntrinsicLibrary.AddProcedure(std::move(sym));
95     }
96   }
97 }
98 
99 template <typename HostT>
100 static void AddLibmComplexHostProcedures(
101     HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
102   using F = FuncPointer<std::complex<HostT>, const std::complex<HostT> &>;
103   using F2 = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
104       const std::complex<HostT> &>;
105   using F2a = FuncPointer<std::complex<HostT>, const HostT &,
106       const std::complex<HostT> &>;
107   using F2b = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
108       const HostT &>;
109   HostRuntimeIntrinsicProcedure libmSymbols[]{
110       {"abs", FuncPointer<HostT, const std::complex<HostT> &>{std::abs}, true},
111       {"acos", F{std::acos}, true},
112       {"acosh", F{std::acosh}, true},
113       {"asin", F{std::asin}, true},
114       {"asinh", F{std::asinh}, true},
115       {"atan", F{std::atan}, true},
116       {"atanh", F{std::atanh}, true},
117       {"cos", F{std::cos}, true},
118       {"cosh", F{std::cosh}, true},
119       {"exp", F{std::exp}, true},
120       {"log", F{std::log}, true},
121       {"pow", F2{std::pow}, true},
122       {"pow", F2a{std::pow}, true},
123       {"pow", F2b{std::pow}, true},
124       {"sin", F{std::sin}, true},
125       {"sinh", F{std::sinh}, true},
126       {"sqrt", F{std::sqrt}, true},
127       {"tan", F{std::tan}, true},
128       {"tanh", F{std::tanh}, true},
129   };
130 
131   for (auto sym : libmSymbols) {
132     if (!hostIntrinsicLibrary.HasEquivalentProcedure(sym)) {
133       hostIntrinsicLibrary.AddProcedure(std::move(sym));
134     }
135   }
136 }
137 
138 void InitHostIntrinsicLibraryWithLibm(HostIntrinsicProceduresLibrary &lib) {
139   if constexpr (host::FortranTypeExists<float>()) {
140     AddLibmRealHostProcedures<float>(lib);
141   }
142   if constexpr (host::FortranTypeExists<double>()) {
143     AddLibmRealHostProcedures<double>(lib);
144   }
145   if constexpr (host::FortranTypeExists<long double>()) {
146     AddLibmRealHostProcedures<long double>(lib);
147   }
148 
149   if constexpr (host::FortranTypeExists<std::complex<float>>()) {
150     AddLibmComplexHostProcedures<float>(lib);
151   }
152   if constexpr (host::FortranTypeExists<std::complex<double>>()) {
153     AddLibmComplexHostProcedures<double>(lib);
154   }
155   if constexpr (host::FortranTypeExists<std::complex<long double>>()) {
156     AddLibmComplexHostProcedures<long double>(lib);
157   }
158 }
159 
160 #if LINK_WITH_LIBPGMATH
161 namespace pgmath {
162 // Define mapping between numerical intrinsics and libpgmath symbols
163 // namespace is used to have shorter names on repeated patterns.
164 // A class would be better to hold all these defs, but GCC does not
165 // support specialization of template variables inside class even
166 // if it is C++14 standard compliant here because there are only full
167 // specializations.
168 
169 // List of intrinsics that have libpgmath implementations that can be used for
170 // constant folding The tag names must match the name used inside libpgmath name
171 // so that the macro below work.
172 enum class I {
173   acos,
174   acosh,
175   asin,
176   asinh,
177   atan,
178   atan2,
179   atanh,
180   bessel_j0,
181   bessel_j1,
182   bessel_jn,
183   bessel_y0,
184   bessel_y1,
185   bessel_yn,
186   cos,
187   cosh,
188   erf,
189   erfc,
190   erfc_scaled,
191   exp,
192   gamma,
193   hypot,
194   log,
195   log10,
196   log_gamma,
197   mod,
198   pow,
199   sin,
200   sinh,
201   sqrt,
202   tan,
203   tanh
204 };
205 
206 // Library versions: P for Precise, R for Relaxed, F for Fast
207 enum class L { F, R, P };
208 
209 struct NoSuchRuntimeSymbol {};
210 template <L, I, typename> constexpr auto Sym{NoSuchRuntimeSymbol{}};
211 
212 // Macros to declare fast/relaxed/precise libpgmath variants.
213 #define DECLARE_PGMATH_FAST_REAL(func) \
214   extern "C" float __fs_##func##_1(float); \
215   extern "C" double __fd_##func##_1(double); \
216   template <> constexpr auto Sym<L::F, I::func, float>{__fs_##func##_1}; \
217   template <> constexpr auto Sym<L::F, I::func, double>{__fd_##func##_1};
218 
219 #define DECLARE_PGMATH_FAST_COMPLEX(func) \
220   extern "C" float _Complex __fc_##func##_1(float _Complex); \
221   extern "C" double _Complex __fz_##func##_1(double _Complex); \
222   template <> \
223   constexpr auto Sym<L::F, I::func, std::complex<float>>{__fc_##func##_1}; \
224   template <> \
225   constexpr auto Sym<L::F, I::func, std::complex<double>>{__fz_##func##_1};
226 
227 #define DECLARE_PGMATH_FAST_ALL_FP(func) \
228   DECLARE_PGMATH_FAST_REAL(func) \
229   DECLARE_PGMATH_FAST_COMPLEX(func)
230 
231 #define DECLARE_PGMATH_PRECISE_REAL(func) \
232   extern "C" float __ps_##func##_1(float); \
233   extern "C" double __pd_##func##_1(double); \
234   template <> constexpr auto Sym<L::P, I::func, float>{__ps_##func##_1}; \
235   template <> constexpr auto Sym<L::P, I::func, double>{__pd_##func##_1};
236 
237 #define DECLARE_PGMATH_PRECISE_COMPLEX(func) \
238   extern "C" float _Complex __pc_##func##_1(float _Complex); \
239   extern "C" double _Complex __pz_##func##_1(double _Complex); \
240   template <> \
241   constexpr auto Sym<L::P, I::func, std::complex<float>>{__pc_##func##_1}; \
242   template <> \
243   constexpr auto Sym<L::P, I::func, std::complex<double>>{__pz_##func##_1};
244 
245 #define DECLARE_PGMATH_PRECISE_ALL_FP(func) \
246   DECLARE_PGMATH_PRECISE_REAL(func) \
247   DECLARE_PGMATH_PRECISE_COMPLEX(func)
248 
249 #define DECLARE_PGMATH_RELAXED_REAL(func) \
250   extern "C" float __rs_##func##_1(float); \
251   extern "C" double __rd_##func##_1(double); \
252   template <> constexpr auto Sym<L::R, I::func, float>{__rs_##func##_1}; \
253   template <> constexpr auto Sym<L::R, I::func, double>{__rd_##func##_1};
254 
255 #define DECLARE_PGMATH_RELAXED_COMPLEX(func) \
256   extern "C" float _Complex __rc_##func##_1(float _Complex); \
257   extern "C" double _Complex __rz_##func##_1(double _Complex); \
258   template <> \
259   constexpr auto Sym<L::R, I::func, std::complex<float>>{__rc_##func##_1}; \
260   template <> \
261   constexpr auto Sym<L::R, I::func, std::complex<double>>{__rz_##func##_1};
262 
263 #define DECLARE_PGMATH_RELAXED_ALL_FP(func) \
264   DECLARE_PGMATH_RELAXED_REAL(func) \
265   DECLARE_PGMATH_RELAXED_COMPLEX(func)
266 
267 #define DECLARE_PGMATH_REAL(func) \
268   DECLARE_PGMATH_FAST_REAL(func) \
269   DECLARE_PGMATH_PRECISE_REAL(func) \
270   DECLARE_PGMATH_RELAXED_REAL(func)
271 
272 #define DECLARE_PGMATH_COMPLEX(func) \
273   DECLARE_PGMATH_FAST_COMPLEX(func) \
274   DECLARE_PGMATH_PRECISE_COMPLEX(func) \
275   DECLARE_PGMATH_RELAXED_COMPLEX(func)
276 
277 #define DECLARE_PGMATH_ALL(func) \
278   DECLARE_PGMATH_REAL(func) \
279   DECLARE_PGMATH_COMPLEX(func)
280 
281 // Macros to declare fast/relaxed/precise libpgmath variants with two arguments.
282 #define DECLARE_PGMATH_FAST_REAL2(func) \
283   extern "C" float __fs_##func##_1(float, float); \
284   extern "C" double __fd_##func##_1(double, double); \
285   template <> constexpr auto Sym<L::F, I::func, float>{__fs_##func##_1}; \
286   template <> constexpr auto Sym<L::F, I::func, double>{__fd_##func##_1};
287 
288 #define DECLARE_PGMATH_FAST_COMPLEX2(func) \
289   extern "C" float _Complex __fc_##func##_1(float _Complex, float _Complex); \
290   extern "C" double _Complex __fz_##func##_1( \
291       double _Complex, double _Complex); \
292   template <> \
293   constexpr auto Sym<L::F, I::func, std::complex<float>>{__fc_##func##_1}; \
294   template <> \
295   constexpr auto Sym<L::F, I::func, std::complex<double>>{__fz_##func##_1};
296 
297 #define DECLARE_PGMATH_FAST_ALL_FP2(func) \
298   DECLARE_PGMATH_FAST_REAL2(func) \
299   DECLARE_PGMATH_FAST_COMPLEX2(func)
300 
301 #define DECLARE_PGMATH_PRECISE_REAL2(func) \
302   extern "C" float __ps_##func##_1(float, float); \
303   extern "C" double __pd_##func##_1(double, double); \
304   template <> constexpr auto Sym<L::P, I::func, float>{__ps_##func##_1}; \
305   template <> constexpr auto Sym<L::P, I::func, double>{__pd_##func##_1};
306 
307 #define DECLARE_PGMATH_PRECISE_COMPLEX2(func) \
308   extern "C" float _Complex __pc_##func##_1(float _Complex, float _Complex); \
309   extern "C" double _Complex __pz_##func##_1( \
310       double _Complex, double _Complex); \
311   template <> \
312   constexpr auto Sym<L::P, I::func, std::complex<float>>{__pc_##func##_1}; \
313   template <> \
314   constexpr auto Sym<L::P, I::func, std::complex<double>>{__pz_##func##_1};
315 
316 #define DECLARE_PGMATH_PRECISE_ALL_FP2(func) \
317   DECLARE_PGMATH_PRECISE_REAL2(func) \
318   DECLARE_PGMATH_PRECISE_COMPLEX2(func)
319 
320 #define DECLARE_PGMATH_RELAXED_REAL2(func) \
321   extern "C" float __rs_##func##_1(float, float); \
322   extern "C" double __rd_##func##_1(double, double); \
323   template <> constexpr auto Sym<L::R, I::func, float>{__rs_##func##_1}; \
324   template <> constexpr auto Sym<L::R, I::func, double>{__rd_##func##_1};
325 
326 #define DECLARE_PGMATH_RELAXED_COMPLEX2(func) \
327   extern "C" float _Complex __rc_##func##_1(float _Complex, float _Complex); \
328   extern "C" double _Complex __rz_##func##_1( \
329       double _Complex, double _Complex); \
330   template <> \
331   constexpr auto Sym<L::R, I::func, std::complex<float>>{__rc_##func##_1}; \
332   template <> \
333   constexpr auto Sym<L::R, I::func, std::complex<double>>{__rz_##func##_1};
334 
335 #define DECLARE_PGMATH_RELAXED_ALL_FP2(func) \
336   DECLARE_PGMATH_RELAXED_REAL2(func) \
337   DECLARE_PGMATH_RELAXED_COMPLEX2(func)
338 
339 #define DECLARE_PGMATH_REAL2(func) \
340   DECLARE_PGMATH_FAST_REAL2(func) \
341   DECLARE_PGMATH_PRECISE_REAL2(func) \
342   DECLARE_PGMATH_RELAXED_REAL2(func)
343 
344 #define DECLARE_PGMATH_COMPLEX2(func) \
345   DECLARE_PGMATH_FAST_COMPLEX2(func) \
346   DECLARE_PGMATH_PRECISE_COMPLEX2(func) \
347   DECLARE_PGMATH_RELAXED_COMPLEX2(func)
348 
349 #define DECLARE_PGMATH_ALL2(func) \
350   DECLARE_PGMATH_REAL2(func) \
351   DECLARE_PGMATH_COMPLEX2(func)
352 
353 // Marcos to declare __mth_i libpgmath variants
354 #define DECLARE_PGMATH_MTH_VERSION_REAL(func) \
355   extern "C" float __mth_i_##func(float); \
356   extern "C" double __mth_i_d##func(double); \
357   template <> constexpr auto Sym<L::F, I::func, float>{__mth_i_##func}; \
358   template <> constexpr auto Sym<L::F, I::func, double>{__mth_i_d##func}; \
359   template <> constexpr auto Sym<L::P, I::func, float>{__mth_i_##func}; \
360   template <> constexpr auto Sym<L::P, I::func, double>{__mth_i_d##func}; \
361   template <> constexpr auto Sym<L::R, I::func, float>{__mth_i_##func}; \
362   template <> constexpr auto Sym<L::R, I::func, double>{__mth_i_d##func};
363 
364 // Actual libpgmath declarations
365 DECLARE_PGMATH_ALL(acos)
366 DECLARE_PGMATH_MTH_VERSION_REAL(acosh)
367 DECLARE_PGMATH_ALL(asin)
368 DECLARE_PGMATH_MTH_VERSION_REAL(asinh)
369 DECLARE_PGMATH_ALL(atan)
370 DECLARE_PGMATH_REAL2(atan2)
371 DECLARE_PGMATH_MTH_VERSION_REAL(atanh)
372 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_j0)
373 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_j1)
374 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_y0)
375 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_y1)
376 // bessel_jn and bessel_yn takes an int as first arg
377 extern "C" float __mth_i_bessel_jn(int, float);
378 extern "C" double __mth_i_dbessel_jn(int, double);
379 template <> constexpr auto Sym<L::F, I::bessel_jn, float>{__mth_i_bessel_jn};
380 template <> constexpr auto Sym<L::F, I::bessel_jn, double>{__mth_i_dbessel_jn};
381 template <> constexpr auto Sym<L::P, I::bessel_jn, float>{__mth_i_bessel_jn};
382 template <> constexpr auto Sym<L::P, I::bessel_jn, double>{__mth_i_dbessel_jn};
383 template <> constexpr auto Sym<L::R, I::bessel_jn, float>{__mth_i_bessel_jn};
384 template <> constexpr auto Sym<L::R, I::bessel_jn, double>{__mth_i_dbessel_jn};
385 extern "C" float __mth_i_bessel_yn(int, float);
386 extern "C" double __mth_i_dbessel_yn(int, double);
387 template <> constexpr auto Sym<L::F, I::bessel_yn, float>{__mth_i_bessel_yn};
388 template <> constexpr auto Sym<L::F, I::bessel_yn, double>{__mth_i_dbessel_yn};
389 template <> constexpr auto Sym<L::P, I::bessel_yn, float>{__mth_i_bessel_yn};
390 template <> constexpr auto Sym<L::P, I::bessel_yn, double>{__mth_i_dbessel_yn};
391 template <> constexpr auto Sym<L::R, I::bessel_yn, float>{__mth_i_bessel_yn};
392 template <> constexpr auto Sym<L::R, I::bessel_yn, double>{__mth_i_dbessel_yn};
393 DECLARE_PGMATH_ALL(cos)
394 DECLARE_PGMATH_ALL(cosh)
395 DECLARE_PGMATH_MTH_VERSION_REAL(erf)
396 DECLARE_PGMATH_MTH_VERSION_REAL(erfc)
397 DECLARE_PGMATH_MTH_VERSION_REAL(erfc_scaled)
398 DECLARE_PGMATH_ALL(exp)
399 DECLARE_PGMATH_MTH_VERSION_REAL(gamma)
400 extern "C" float __mth_i_hypot(float, float);
401 extern "C" double __mth_i_dhypot(double, double);
402 template <> constexpr auto Sym<L::F, I::hypot, float>{__mth_i_hypot};
403 template <> constexpr auto Sym<L::F, I::hypot, double>{__mth_i_dhypot};
404 template <> constexpr auto Sym<L::P, I::hypot, float>{__mth_i_hypot};
405 template <> constexpr auto Sym<L::P, I::hypot, double>{__mth_i_dhypot};
406 template <> constexpr auto Sym<L::R, I::hypot, float>{__mth_i_hypot};
407 template <> constexpr auto Sym<L::R, I::hypot, double>{__mth_i_dhypot};
408 DECLARE_PGMATH_ALL(log)
409 DECLARE_PGMATH_REAL(log10)
410 DECLARE_PGMATH_MTH_VERSION_REAL(log_gamma)
411 // no function for modulo in libpgmath
412 extern "C" float __fs_mod_1(float, float);
413 extern "C" double __fd_mod_1(double, double);
414 template <> constexpr auto Sym<L::F, I::mod, float>{__fs_mod_1};
415 template <> constexpr auto Sym<L::F, I::mod, double>{__fd_mod_1};
416 template <> constexpr auto Sym<L::P, I::mod, float>{__fs_mod_1};
417 template <> constexpr auto Sym<L::P, I::mod, double>{__fd_mod_1};
418 template <> constexpr auto Sym<L::R, I::mod, float>{__fs_mod_1};
419 template <> constexpr auto Sym<L::R, I::mod, double>{__fd_mod_1};
420 DECLARE_PGMATH_ALL2(pow)
421 DECLARE_PGMATH_ALL(sin)
422 DECLARE_PGMATH_ALL(sinh)
423 DECLARE_PGMATH_MTH_VERSION_REAL(sqrt)
424 DECLARE_PGMATH_COMPLEX(sqrt) // real versions are __mth_i...
425 DECLARE_PGMATH_ALL(tan)
426 DECLARE_PGMATH_ALL(tanh)
427 
428 // Fill the function map used for folding with libpgmath symbols
429 template <L Lib, typename HostT>
430 static void AddLibpgmathRealHostProcedures(
431     HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
432   static_assert(std::is_same_v<HostT, float> || std::is_same_v<HostT, double>);
433   HostRuntimeIntrinsicProcedure pgmathSymbols[]{
434       {"acos", Sym<Lib, I::acos, HostT>, true},
435       {"acosh", Sym<Lib, I::acosh, HostT>, true},
436       {"asin", Sym<Lib, I::asin, HostT>, true},
437       {"asinh", Sym<Lib, I::asinh, HostT>, true},
438       {"atan", Sym<Lib, I::atan, HostT>, true},
439       {"atan", Sym<Lib, I::atan2, HostT>,
440           true}, // atan is also the generic name for atan2
441       {"atanh", Sym<Lib, I::atanh, HostT>, true},
442       {"bessel_j0", Sym<Lib, I::bessel_j0, HostT>, true},
443       {"bessel_j1", Sym<Lib, I::bessel_j1, HostT>, true},
444       {"bessel_jn", Sym<Lib, I::bessel_jn, HostT>, true},
445       {"bessel_y0", Sym<Lib, I::bessel_y0, HostT>, true},
446       {"bessel_y1", Sym<Lib, I::bessel_y1, HostT>, true},
447       {"bessel_yn", Sym<Lib, I::bessel_yn, HostT>, true},
448       {"cos", Sym<Lib, I::cos, HostT>, true},
449       {"cosh", Sym<Lib, I::cosh, HostT>, true},
450       {"erf", Sym<Lib, I::erf, HostT>, true},
451       {"erfc", Sym<Lib, I::erfc, HostT>, true},
452       {"erfc_scaled", Sym<Lib, I::erfc_scaled, HostT>, true},
453       {"exp", Sym<Lib, I::exp, HostT>, true},
454       {"gamma", Sym<Lib, I::gamma, HostT>, true},
455       {"hypot", Sym<Lib, I::hypot, HostT>, true},
456       {"log", Sym<Lib, I::log, HostT>, true},
457       {"log10", Sym<Lib, I::log10, HostT>, true},
458       {"log_gamma", Sym<Lib, I::log_gamma, HostT>, true},
459       {"mod", Sym<Lib, I::mod, HostT>, true},
460       {"pow", Sym<Lib, I::pow, HostT>, true},
461       {"sin", Sym<Lib, I::sin, HostT>, true},
462       {"sinh", Sym<Lib, I::sinh, HostT>, true},
463       {"sqrt", Sym<Lib, I::sqrt, HostT>, true},
464       {"tan", Sym<Lib, I::tan, HostT>, true},
465       {"tanh", Sym<Lib, I::tanh, HostT>, true},
466   };
467 
468   for (auto sym : pgmathSymbols) {
469     hostIntrinsicLibrary.AddProcedure(std::move(sym));
470   }
471 }
472 
473 // Note: std::complex and _complex are layout compatible but are not guaranteed
474 // to be linkage compatible. For instance, on i386, float _Complex is returned
475 // by a pair of register but std::complex<float> is returned by structure
476 // address. To fix the issue, wrapper around C _Complex functions are defined
477 // below.
478 template <FuncPointer<float _Complex, float _Complex> func>
479 static std::complex<float> ComplexCFuncWrapper(std::complex<float> &arg) {
480   float _Complex res{func(*reinterpret_cast<float _Complex *>(&arg))};
481   return *reinterpret_cast<std::complex<float> *>(&res);
482 }
483 
484 template <FuncPointer<double _Complex, double _Complex> func>
485 static std::complex<double> ComplexCFuncWrapper(std::complex<double> &arg) {
486   double _Complex res{func(*reinterpret_cast<double _Complex *>(&arg))};
487   return *reinterpret_cast<std::complex<double> *>(&res);
488 }
489 
490 template <FuncPointer<float _Complex, float _Complex, float _Complex> func>
491 static std::complex<float> ComplexCFuncWrapper(
492     std::complex<float> &arg1, std::complex<float> &arg2) {
493   float _Complex res{func(*reinterpret_cast<float _Complex *>(&arg1),
494       *reinterpret_cast<float _Complex *>(&arg2))};
495   return *reinterpret_cast<std::complex<float> *>(&res);
496 }
497 
498 template <FuncPointer<double _Complex, double _Complex, double _Complex> func>
499 static std::complex<double> ComplexCFuncWrapper(
500     std::complex<double> &arg1, std::complex<double> &arg2) {
501   double _Complex res{func(*reinterpret_cast<double _Complex *>(&arg1),
502       *reinterpret_cast<double _Complex *>(&arg2))};
503   return *reinterpret_cast<std::complex<double> *>(&res);
504 }
505 
506 template <L Lib, typename HostT>
507 static void AddLibpgmathComplexHostProcedures(
508     HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
509   static_assert(std::is_same_v<HostT, float> || std::is_same_v<HostT, double>);
510   using CHostT = std::complex<HostT>;
511   // cmath is used to complement pgmath when symbols are not available
512   using CmathF = FuncPointer<CHostT, const CHostT &>;
513   HostRuntimeIntrinsicProcedure pgmathSymbols[]{
514       {"abs", FuncPointer<HostT, const CHostT &>{std::abs}, true},
515       {"acos", ComplexCFuncWrapper<Sym<Lib, I::acos, CHostT>>, true},
516       {"acosh", CmathF{std::acosh}, true},
517       {"asin", ComplexCFuncWrapper<Sym<Lib, I::asin, CHostT>>, true},
518       {"asinh", CmathF{std::asinh}, true},
519       {"atan", ComplexCFuncWrapper<Sym<Lib, I::atan, CHostT>>, true},
520       {"atanh", CmathF{std::atanh}, true},
521       {"cos", ComplexCFuncWrapper<Sym<Lib, I::cos, CHostT>>, true},
522       {"cosh", ComplexCFuncWrapper<Sym<Lib, I::cosh, CHostT>>, true},
523       {"exp", ComplexCFuncWrapper<Sym<Lib, I::exp, CHostT>>, true},
524       {"log", ComplexCFuncWrapper<Sym<Lib, I::log, CHostT>>, true},
525       {"pow", ComplexCFuncWrapper<Sym<Lib, I::pow, CHostT>>, true},
526       {"sin", ComplexCFuncWrapper<Sym<Lib, I::sin, CHostT>>, true},
527       {"sinh", ComplexCFuncWrapper<Sym<Lib, I::sinh, CHostT>>, true},
528       {"sqrt", ComplexCFuncWrapper<Sym<Lib, I::sqrt, CHostT>>, true},
529       {"tan", ComplexCFuncWrapper<Sym<Lib, I::tan, CHostT>>, true},
530       {"tanh", ComplexCFuncWrapper<Sym<Lib, I::tanh, CHostT>>, true},
531   };
532 
533   for (auto sym : pgmathSymbols) {
534     hostIntrinsicLibrary.AddProcedure(std::move(sym));
535   }
536 }
537 
538 template <L Lib>
539 static void InitHostIntrinsicLibraryWithLibpgmath(
540     HostIntrinsicProceduresLibrary &lib) {
541   if constexpr (host::FortranTypeExists<float>()) {
542     AddLibpgmathRealHostProcedures<Lib, float>(lib);
543   }
544   if constexpr (host::FortranTypeExists<double>()) {
545     AddLibpgmathRealHostProcedures<Lib, double>(lib);
546   }
547   if constexpr (host::FortranTypeExists<std::complex<float>>()) {
548     AddLibpgmathComplexHostProcedures<Lib, float>(lib);
549   }
550   if constexpr (host::FortranTypeExists<std::complex<double>>()) {
551     AddLibpgmathComplexHostProcedures<Lib, double>(lib);
552   }
553   // No long double functions in libpgmath
554   if constexpr (host::FortranTypeExists<long double>()) {
555     AddLibmRealHostProcedures<long double>(lib);
556   }
557   if constexpr (host::FortranTypeExists<std::complex<long double>>()) {
558     AddLibmComplexHostProcedures<long double>(lib);
559   }
560 }
561 } // namespace pgmath
562 #endif // LINK_WITH_LIBPGMATH
563 
564 // Define which host runtime functions will be used for folding
565 HostIntrinsicProceduresLibrary::HostIntrinsicProceduresLibrary() {
566   // TODO: When command line options regarding targeted numerical library is
567   // available, this needs to be revisited to take it into account. So far,
568   // default to libpgmath if F18 is built with it.
569 #if LINK_WITH_LIBPGMATH
570   // This looks and is stupid for now (until TODO above), but it is needed
571   // to silence clang warnings on unused symbols if all declared pgmath
572   // symbols are not used somewhere.
573   if (true) {
574     pgmath::InitHostIntrinsicLibraryWithLibpgmath<pgmath::L::P>(*this);
575   } else if (false) {
576     pgmath::InitHostIntrinsicLibraryWithLibpgmath<pgmath::L::F>(*this);
577   } else {
578     pgmath::InitHostIntrinsicLibraryWithLibpgmath<pgmath::L::R>(*this);
579   }
580 #else
581   InitHostIntrinsicLibraryWithLibm(*this);
582 #endif
583 }
584 } // namespace Fortran::evaluate
585