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
68template <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
74template <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
83template <class _Tp>
84complex<_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
108template <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
114template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x) {
115  return log(__x) / log(_Tp(10));
116}
117
118// sqrt
119
120template <class _Tp>
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
137template <class _Tp>
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
157template <class _Tp>
158std::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
165template <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
173template <class _Tp>
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
201template <class _Tp>
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
234template <class _Tp>
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
264template <class _Tp>
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
278template <class _Tp>
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
294template <class _Tp>
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
316template <class _Tp>
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
324template <class _Tp>
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
358template <class _Tp>
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
366template <class _Tp>
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
374template <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
380template <class _Tp>
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