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