complex_cmath.h (12559B) - Raw
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 #ifdef _GLIBCXX20_CONSTEXPR 68 #define CXX20_CONSTEXPR_DEVICE __DEVICE__ 69 #else 70 #define CXX20_CONSTEXPR_DEVICE 71 #endif 72 template <class _Tp> 73 CXX20_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 79 template <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 88 template <class _Tp> 89 complex<_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 113 template <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 119 template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x) { 120 return log(__x) / log(_Tp(10)); 121 } 122 123 // sqrt 124 125 template <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 142 template <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 162 template <class _Tp> 163 std::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 170 template <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 178 template <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 206 template <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 239 template <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 269 template <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 283 template <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 299 template <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 321 template <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 329 template <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 363 template <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 371 template <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 379 template <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 385 template <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