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