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