clang/lib/Headers/openmp_wrappers/complex_cmath.h

   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