1 //===------------------------- __complex_cmath.h --------------------------===//
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 // std::complex header copied from the libcxx source and simplified for use in
10 // OpenMP target offload regions.
11 //
12 //===----------------------------------------------------------------------===//
13 
14 #ifndef _OPENMP
15 #error "This file is for OpenMP compilation only."
16 #endif
17 
18 #ifndef __cplusplus
19 #error "This file is for C++ compilation only."
20 #endif
21 
22 #ifndef _LIBCPP_COMPLEX
23 #define _LIBCPP_COMPLEX
24 
25 #include <cmath>
26 #include <type_traits>
27 
28 #define __DEVICE__ static constexpr __attribute__((nothrow))
29 
30 namespace std {
31 
32 // abs
33 
abs(const std::complex<_Tp> & __c)34 template <class _Tp> __DEVICE__ _Tp abs(const std::complex<_Tp> &__c) {
35   return hypot(__c.real(), __c.imag());
36 }
37 
38 // arg
39 
arg(const std::complex<_Tp> & __c)40 template <class _Tp> __DEVICE__ _Tp arg(const std::complex<_Tp> &__c) {
41   return atan2(__c.imag(), __c.real());
42 }
43 
44 template <class _Tp>
45 typename enable_if<is_integral<_Tp>::value || is_same<_Tp, double>::value,
46                    double>::type
arg(_Tp __re)47 arg(_Tp __re) {
48   return atan2(0., __re);
49 }
50 
51 template <class _Tp>
arg(_Tp __re)52 typename enable_if<is_same<_Tp, float>::value, float>::type arg(_Tp __re) {
53   return atan2f(0.F, __re);
54 }
55 
56 // norm
57 
norm(const std::complex<_Tp> & __c)58 template <class _Tp> __DEVICE__ _Tp norm(const std::complex<_Tp> &__c) {
59   if (std::isinf(__c.real()))
60     return abs(__c.real());
61   if (std::isinf(__c.imag()))
62     return abs(__c.imag());
63   return __c.real() * __c.real() + __c.imag() * __c.imag();
64 }
65 
66 // conj
67 
conj(const std::complex<_Tp> & __c)68 template <class _Tp> std::complex<_Tp> conj(const std::complex<_Tp> &__c) {
69   return std::complex<_Tp>(__c.real(), -__c.imag());
70 }
71 
72 // proj
73 
proj(const std::complex<_Tp> & __c)74 template <class _Tp> std::complex<_Tp> proj(const std::complex<_Tp> &__c) {
75   std::complex<_Tp> __r = __c;
76   if (std::isinf(__c.real()) || std::isinf(__c.imag()))
77     __r = std::complex<_Tp>(INFINITY, copysign(_Tp(0), __c.imag()));
78   return __r;
79 }
80 
81 // polar
82 
83 template <class _Tp>
84 complex<_Tp> polar(const _Tp &__rho, const _Tp &__theta = _Tp()) {
85   if (std::isnan(__rho) || signbit(__rho))
86     return std::complex<_Tp>(_Tp(NAN), _Tp(NAN));
87   if (std::isnan(__theta)) {
88     if (std::isinf(__rho))
89       return std::complex<_Tp>(__rho, __theta);
90     return std::complex<_Tp>(__theta, __theta);
91   }
92   if (std::isinf(__theta)) {
93     if (std::isinf(__rho))
94       return std::complex<_Tp>(__rho, _Tp(NAN));
95     return std::complex<_Tp>(_Tp(NAN), _Tp(NAN));
96   }
97   _Tp __x = __rho * cos(__theta);
98   if (std::isnan(__x))
99     __x = 0;
100   _Tp __y = __rho * sin(__theta);
101   if (std::isnan(__y))
102     __y = 0;
103   return std::complex<_Tp>(__x, __y);
104 }
105 
106 // log
107 
log(const std::complex<_Tp> & __x)108 template <class _Tp> std::complex<_Tp> log(const std::complex<_Tp> &__x) {
109   return std::complex<_Tp>(log(abs(__x)), arg(__x));
110 }
111 
112 // log10
113 
log10(const std::complex<_Tp> & __x)114 template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x) {
115   return log(__x) / log(_Tp(10));
116 }
117 
118 // sqrt
119 
120 template <class _Tp>
sqrt(const std::complex<_Tp> & __x)121 __DEVICE__ std::complex<_Tp> sqrt(const std::complex<_Tp> &__x) {
122   if (std::isinf(__x.imag()))
123     return std::complex<_Tp>(_Tp(INFINITY), __x.imag());
124   if (std::isinf(__x.real())) {
125     if (__x.real() > _Tp(0))
126       return std::complex<_Tp>(__x.real(), std::isnan(__x.imag())
127                                                ? __x.imag()
128                                                : copysign(_Tp(0), __x.imag()));
129     return std::complex<_Tp>(std::isnan(__x.imag()) ? __x.imag() : _Tp(0),
130                              copysign(__x.real(), __x.imag()));
131   }
132   return polar(sqrt(abs(__x)), arg(__x) / _Tp(2));
133 }
134 
135 // exp
136 
137 template <class _Tp>
exp(const std::complex<_Tp> & __x)138 __DEVICE__ std::complex<_Tp> exp(const std::complex<_Tp> &__x) {
139   _Tp __i = __x.imag();
140   if (std::isinf(__x.real())) {
141     if (__x.real() < _Tp(0)) {
142       if (!std::isfinite(__i))
143         __i = _Tp(1);
144     } else if (__i == 0 || !std::isfinite(__i)) {
145       if (std::isinf(__i))
146         __i = _Tp(NAN);
147       return std::complex<_Tp>(__x.real(), __i);
148     }
149   } else if (std::isnan(__x.real()) && __x.imag() == 0)
150     return __x;
151   _Tp __e = exp(__x.real());
152   return std::complex<_Tp>(__e * cos(__i), __e * sin(__i));
153 }
154 
155 // pow
156 
157 template <class _Tp>
pow(const std::complex<_Tp> & __x,const std::complex<_Tp> & __y)158 std::complex<_Tp> pow(const std::complex<_Tp> &__x,
159                       const std::complex<_Tp> &__y) {
160   return exp(__y * log(__x));
161 }
162 
163 // __sqr, computes pow(x, 2)
164 
__sqr(const std::complex<_Tp> & __x)165 template <class _Tp> std::complex<_Tp> __sqr(const std::complex<_Tp> &__x) {
166   return std::complex<_Tp>((__x.real() - __x.imag()) *
167                                (__x.real() + __x.imag()),
168                            _Tp(2) * __x.real() * __x.imag());
169 }
170 
171 // asinh
172 
173 template <class _Tp>
asinh(const std::complex<_Tp> & __x)174 __DEVICE__ std::complex<_Tp> asinh(const std::complex<_Tp> &__x) {
175   const _Tp __pi(atan2(+0., -0.));
176   if (std::isinf(__x.real())) {
177     if (std::isnan(__x.imag()))
178       return __x;
179     if (std::isinf(__x.imag()))
180       return std::complex<_Tp>(__x.real(),
181                                copysign(__pi * _Tp(0.25), __x.imag()));
182     return std::complex<_Tp>(__x.real(), copysign(_Tp(0), __x.imag()));
183   }
184   if (std::isnan(__x.real())) {
185     if (std::isinf(__x.imag()))
186       return std::complex<_Tp>(__x.imag(), __x.real());
187     if (__x.imag() == 0)
188       return __x;
189     return std::complex<_Tp>(__x.real(), __x.real());
190   }
191   if (std::isinf(__x.imag()))
192     return std::complex<_Tp>(copysign(__x.imag(), __x.real()),
193                              copysign(__pi / _Tp(2), __x.imag()));
194   std::complex<_Tp> __z = log(__x + sqrt(__sqr(__x) + _Tp(1)));
195   return std::complex<_Tp>(copysign(__z.real(), __x.real()),
196                            copysign(__z.imag(), __x.imag()));
197 }
198 
199 // acosh
200 
201 template <class _Tp>
acosh(const std::complex<_Tp> & __x)202 __DEVICE__ std::complex<_Tp> acosh(const std::complex<_Tp> &__x) {
203   const _Tp __pi(atan2(+0., -0.));
204   if (std::isinf(__x.real())) {
205     if (std::isnan(__x.imag()))
206       return std::complex<_Tp>(abs(__x.real()), __x.imag());
207     if (std::isinf(__x.imag())) {
208       if (__x.real() > 0)
209         return std::complex<_Tp>(__x.real(),
210                                  copysign(__pi * _Tp(0.25), __x.imag()));
211       else
212         return std::complex<_Tp>(-__x.real(),
213                                  copysign(__pi * _Tp(0.75), __x.imag()));
214     }
215     if (__x.real() < 0)
216       return std::complex<_Tp>(-__x.real(), copysign(__pi, __x.imag()));
217     return std::complex<_Tp>(__x.real(), copysign(_Tp(0), __x.imag()));
218   }
219   if (std::isnan(__x.real())) {
220     if (std::isinf(__x.imag()))
221       return std::complex<_Tp>(abs(__x.imag()), __x.real());
222     return std::complex<_Tp>(__x.real(), __x.real());
223   }
224   if (std::isinf(__x.imag()))
225     return std::complex<_Tp>(abs(__x.imag()),
226                              copysign(__pi / _Tp(2), __x.imag()));
227   std::complex<_Tp> __z = log(__x + sqrt(__sqr(__x) - _Tp(1)));
228   return std::complex<_Tp>(copysign(__z.real(), _Tp(0)),
229                            copysign(__z.imag(), __x.imag()));
230 }
231 
232 // atanh
233 
234 template <class _Tp>
atanh(const std::complex<_Tp> & __x)235 __DEVICE__ std::complex<_Tp> atanh(const std::complex<_Tp> &__x) {
236   const _Tp __pi(atan2(+0., -0.));
237   if (std::isinf(__x.imag())) {
238     return std::complex<_Tp>(copysign(_Tp(0), __x.real()),
239                              copysign(__pi / _Tp(2), __x.imag()));
240   }
241   if (std::isnan(__x.imag())) {
242     if (std::isinf(__x.real()) || __x.real() == 0)
243       return std::complex<_Tp>(copysign(_Tp(0), __x.real()), __x.imag());
244     return std::complex<_Tp>(__x.imag(), __x.imag());
245   }
246   if (std::isnan(__x.real())) {
247     return std::complex<_Tp>(__x.real(), __x.real());
248   }
249   if (std::isinf(__x.real())) {
250     return std::complex<_Tp>(copysign(_Tp(0), __x.real()),
251                              copysign(__pi / _Tp(2), __x.imag()));
252   }
253   if (abs(__x.real()) == _Tp(1) && __x.imag() == _Tp(0)) {
254     return std::complex<_Tp>(copysign(_Tp(INFINITY), __x.real()),
255                              copysign(_Tp(0), __x.imag()));
256   }
257   std::complex<_Tp> __z = log((_Tp(1) + __x) / (_Tp(1) - __x)) / _Tp(2);
258   return std::complex<_Tp>(copysign(__z.real(), __x.real()),
259                            copysign(__z.imag(), __x.imag()));
260 }
261 
262 // sinh
263 
264 template <class _Tp>
sinh(const std::complex<_Tp> & __x)265 __DEVICE__ std::complex<_Tp> sinh(const std::complex<_Tp> &__x) {
266   if (std::isinf(__x.real()) && !std::isfinite(__x.imag()))
267     return std::complex<_Tp>(__x.real(), _Tp(NAN));
268   if (__x.real() == 0 && !std::isfinite(__x.imag()))
269     return std::complex<_Tp>(__x.real(), _Tp(NAN));
270   if (__x.imag() == 0 && !std::isfinite(__x.real()))
271     return __x;
272   return std::complex<_Tp>(sinh(__x.real()) * cos(__x.imag()),
273                            cosh(__x.real()) * sin(__x.imag()));
274 }
275 
276 // cosh
277 
278 template <class _Tp>
cosh(const std::complex<_Tp> & __x)279 __DEVICE__ std::complex<_Tp> cosh(const std::complex<_Tp> &__x) {
280   if (std::isinf(__x.real()) && !std::isfinite(__x.imag()))
281     return std::complex<_Tp>(abs(__x.real()), _Tp(NAN));
282   if (__x.real() == 0 && !std::isfinite(__x.imag()))
283     return std::complex<_Tp>(_Tp(NAN), __x.real());
284   if (__x.real() == 0 && __x.imag() == 0)
285     return std::complex<_Tp>(_Tp(1), __x.imag());
286   if (__x.imag() == 0 && !std::isfinite(__x.real()))
287     return std::complex<_Tp>(abs(__x.real()), __x.imag());
288   return std::complex<_Tp>(cosh(__x.real()) * cos(__x.imag()),
289                            sinh(__x.real()) * sin(__x.imag()));
290 }
291 
292 // tanh
293 
294 template <class _Tp>
tanh(const std::complex<_Tp> & __x)295 __DEVICE__ std::complex<_Tp> tanh(const std::complex<_Tp> &__x) {
296   if (std::isinf(__x.real())) {
297     if (!std::isfinite(__x.imag()))
298       return std::complex<_Tp>(_Tp(1), _Tp(0));
299     return std::complex<_Tp>(_Tp(1),
300                              copysign(_Tp(0), sin(_Tp(2) * __x.imag())));
301   }
302   if (std::isnan(__x.real()) && __x.imag() == 0)
303     return __x;
304   _Tp __2r(_Tp(2) * __x.real());
305   _Tp __2i(_Tp(2) * __x.imag());
306   _Tp __d(cosh(__2r) + cos(__2i));
307   _Tp __2rsh(sinh(__2r));
308   if (std::isinf(__2rsh) && std::isinf(__d))
309     return std::complex<_Tp>(__2rsh > _Tp(0) ? _Tp(1) : _Tp(-1),
310                              __2i > _Tp(0) ? _Tp(0) : _Tp(-0.));
311   return std::complex<_Tp>(__2rsh / __d, sin(__2i) / __d);
312 }
313 
314 // asin
315 
316 template <class _Tp>
asin(const std::complex<_Tp> & __x)317 __DEVICE__ std::complex<_Tp> asin(const std::complex<_Tp> &__x) {
318   std::complex<_Tp> __z = asinh(complex<_Tp>(-__x.imag(), __x.real()));
319   return std::complex<_Tp>(__z.imag(), -__z.real());
320 }
321 
322 // acos
323 
324 template <class _Tp>
acos(const std::complex<_Tp> & __x)325 __DEVICE__ std::complex<_Tp> acos(const std::complex<_Tp> &__x) {
326   const _Tp __pi(atan2(+0., -0.));
327   if (std::isinf(__x.real())) {
328     if (std::isnan(__x.imag()))
329       return std::complex<_Tp>(__x.imag(), __x.real());
330     if (std::isinf(__x.imag())) {
331       if (__x.real() < _Tp(0))
332         return std::complex<_Tp>(_Tp(0.75) * __pi, -__x.imag());
333       return std::complex<_Tp>(_Tp(0.25) * __pi, -__x.imag());
334     }
335     if (__x.real() < _Tp(0))
336       return std::complex<_Tp>(__pi,
337                                signbit(__x.imag()) ? -__x.real() : __x.real());
338     return std::complex<_Tp>(_Tp(0),
339                              signbit(__x.imag()) ? __x.real() : -__x.real());
340   }
341   if (std::isnan(__x.real())) {
342     if (std::isinf(__x.imag()))
343       return std::complex<_Tp>(__x.real(), -__x.imag());
344     return std::complex<_Tp>(__x.real(), __x.real());
345   }
346   if (std::isinf(__x.imag()))
347     return std::complex<_Tp>(__pi / _Tp(2), -__x.imag());
348   if (__x.real() == 0 && (__x.imag() == 0 || isnan(__x.imag())))
349     return std::complex<_Tp>(__pi / _Tp(2), -__x.imag());
350   std::complex<_Tp> __z = log(__x + sqrt(__sqr(__x) - _Tp(1)));
351   if (signbit(__x.imag()))
352     return std::complex<_Tp>(abs(__z.imag()), abs(__z.real()));
353   return std::complex<_Tp>(abs(__z.imag()), -abs(__z.real()));
354 }
355 
356 // atan
357 
358 template <class _Tp>
atan(const std::complex<_Tp> & __x)359 __DEVICE__ std::complex<_Tp> atan(const std::complex<_Tp> &__x) {
360   std::complex<_Tp> __z = atanh(complex<_Tp>(-__x.imag(), __x.real()));
361   return std::complex<_Tp>(__z.imag(), -__z.real());
362 }
363 
364 // sin
365 
366 template <class _Tp>
sin(const std::complex<_Tp> & __x)367 __DEVICE__ std::complex<_Tp> sin(const std::complex<_Tp> &__x) {
368   std::complex<_Tp> __z = sinh(complex<_Tp>(-__x.imag(), __x.real()));
369   return std::complex<_Tp>(__z.imag(), -__z.real());
370 }
371 
372 // cos
373 
cos(const std::complex<_Tp> & __x)374 template <class _Tp> std::complex<_Tp> cos(const std::complex<_Tp> &__x) {
375   return cosh(complex<_Tp>(-__x.imag(), __x.real()));
376 }
377 
378 // tan
379 
380 template <class _Tp>
tan(const std::complex<_Tp> & __x)381 __DEVICE__ std::complex<_Tp> tan(const std::complex<_Tp> &__x) {
382   std::complex<_Tp> __z = tanh(complex<_Tp>(-__x.imag(), __x.real()));
383   return std::complex<_Tp>(__z.imag(), -__z.real());
384 }
385 
386 } // namespace std
387 
388 #endif
389