394 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8//9// std::complex header copied from the libcxx source and simplified for use in10// OpenMP target offload regions.11//12//===----------------------------------------------------------------------===//13 14#ifndef _OPENMP15#error "This file is for OpenMP compilation only."16#endif17 18#ifndef __cplusplus19#error "This file is for C++ compilation only."20#endif21 22#ifndef _LIBCPP_COMPLEX23#define _LIBCPP_COMPLEX24 25#include <cmath>26#include <type_traits>27 28#define __DEVICE__ static constexpr __attribute__((nothrow))29 30namespace std {31 32// abs33 34template <class _Tp> __DEVICE__ _Tp abs(const std::complex<_Tp> &__c) {35 return hypot(__c.real(), __c.imag());36}37 38// arg39 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>::type47arg(_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// norm57 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// conj67#ifdef _GLIBCXX20_CONSTEXPR68#define CXX20_CONSTEXPR_DEVICE __DEVICE__69#else70#define CXX20_CONSTEXPR_DEVICE71#endif72template <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// proj78 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// polar87 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// log112 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// log10118 119template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x) {120 return log(__x) / log(_Tp(10));121}122 123// sqrt124 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// exp141 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// pow161 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// asinh177 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// acosh205 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 else217 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// atanh238 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// sinh268 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// cosh282 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// tanh298 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// asin320 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// acos328 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// atan362 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// sin370 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// cos378 379template <class _Tp> std::complex<_Tp> cos(const std::complex<_Tp> &__x) {380 return cosh(complex<_Tp>(-__x.imag(), __x.real()));381}382 383// tan384 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 std392 393#endif394