Warning: This file is not a C or C++ file. It does not have highlighting.
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 | |
34 | template <class _Tp> __DEVICE__ _Tp abs(const std::complex<_Tp> &__c) { |
35 | return hypot(__c.real(), __c.imag()); |
36 | } |
37 | |
38 | // arg |
39 | |
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 |
47 | arg(_Tp __re) { |
48 | return atan2(0., __re); |
49 | } |
50 | |
51 | template <class _Tp> |
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 | |
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 | |
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 | |
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 | |
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 | |
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> |
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> |
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> |
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 | |
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> |
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> |
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> |
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> |
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> |
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> |
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> |
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> |
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> |
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> |
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 | |
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> |
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 |
Warning: This file is not a C or C++ file. It does not have highlighting.