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