| #pragma once |
| |
| #include <cmath> |
| |
| #define INFINITY (__builtin_inff()) |
| #define NAN (__builtin_nanf ("")) |
| |
| namespace std { |
| |
| // Taken from libc++ |
| template <class _Tp> |
| class complex { |
| public: |
| typedef _Tp value_type; |
| |
| private: |
| value_type __re_; |
| value_type __im_; |
| |
| public: |
| complex(const value_type &__re = value_type(), const value_type &__im = value_type()) |
| : __re_(__re), __im_(__im) {} |
| template <class _Xp> |
| complex(const complex<_Xp> &__c) |
| : __re_(__c.real()), __im_(__c.imag()) {} |
| |
| value_type real() const { return __re_; } |
| value_type imag() const { return __im_; } |
| |
| void real(value_type __re) { __re_ = __re; } |
| void imag(value_type __im) { __im_ = __im; } |
| |
| complex &operator=(const value_type &__re) { |
| __re_ = __re; |
| __im_ = value_type(); |
| return *this; |
| } |
| complex &operator+=(const value_type &__re) { |
| __re_ += __re; |
| return *this; |
| } |
| complex &operator-=(const value_type &__re) { |
| __re_ -= __re; |
| return *this; |
| } |
| complex &operator*=(const value_type &__re) { |
| __re_ *= __re; |
| __im_ *= __re; |
| return *this; |
| } |
| complex &operator/=(const value_type &__re) { |
| __re_ /= __re; |
| __im_ /= __re; |
| return *this; |
| } |
| |
| template <class _Xp> |
| complex &operator=(const complex<_Xp> &__c) { |
| __re_ = __c.real(); |
| __im_ = __c.imag(); |
| return *this; |
| } |
| template <class _Xp> |
| complex &operator+=(const complex<_Xp> &__c) { |
| __re_ += __c.real(); |
| __im_ += __c.imag(); |
| return *this; |
| } |
| template <class _Xp> |
| complex &operator-=(const complex<_Xp> &__c) { |
| __re_ -= __c.real(); |
| __im_ -= __c.imag(); |
| return *this; |
| } |
| template <class _Xp> |
| complex &operator*=(const complex<_Xp> &__c) { |
| *this = *this * complex(__c.real(), __c.imag()); |
| return *this; |
| } |
| template <class _Xp> |
| complex &operator/=(const complex<_Xp> &__c) { |
| *this = *this / complex(__c.real(), __c.imag()); |
| return *this; |
| } |
| }; |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator+(const complex<_Tp> &__x, const complex<_Tp> &__y) { |
| complex<_Tp> __t(__x); |
| __t += __y; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator+(const complex<_Tp> &__x, const _Tp &__y) { |
| complex<_Tp> __t(__x); |
| __t += __y; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator+(const _Tp &__x, const complex<_Tp> &__y) { |
| complex<_Tp> __t(__y); |
| __t += __x; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator-(const complex<_Tp> &__x, const complex<_Tp> &__y) { |
| complex<_Tp> __t(__x); |
| __t -= __y; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator-(const complex<_Tp> &__x, const _Tp &__y) { |
| complex<_Tp> __t(__x); |
| __t -= __y; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator-(const _Tp &__x, const complex<_Tp> &__y) { |
| complex<_Tp> __t(-__y); |
| __t += __x; |
| return __t; |
| } |
| |
| template <class _Tp> |
| complex<_Tp> |
| operator*(const complex<_Tp> &__z, const complex<_Tp> &__w) { |
| _Tp __a = __z.real(); |
| _Tp __b = __z.imag(); |
| _Tp __c = __w.real(); |
| _Tp __d = __w.imag(); |
| _Tp __ac = __a * __c; |
| _Tp __bd = __b * __d; |
| _Tp __ad = __a * __d; |
| _Tp __bc = __b * __c; |
| _Tp __x = __ac - __bd; |
| _Tp __y = __ad + __bc; |
| if (std::isnan(__x) && std::isnan(__y)) { |
| bool __recalc = false; |
| if (std::isinf(__a) || std::isinf(__b)) { |
| __a = copysign(std::isinf(__a) ? _Tp(1) : _Tp(0), __a); |
| __b = copysign(std::isinf(__b) ? _Tp(1) : _Tp(0), __b); |
| if (std::isnan(__c)) |
| __c = copysign(_Tp(0), __c); |
| if (std::isnan(__d)) |
| __d = copysign(_Tp(0), __d); |
| __recalc = true; |
| } |
| if (std::isinf(__c) || std::isinf(__d)) { |
| __c = copysign(std::isinf(__c) ? _Tp(1) : _Tp(0), __c); |
| __d = copysign(std::isinf(__d) ? _Tp(1) : _Tp(0), __d); |
| if (std::isnan(__a)) |
| __a = copysign(_Tp(0), __a); |
| if (std::isnan(__b)) |
| __b = copysign(_Tp(0), __b); |
| __recalc = true; |
| } |
| if (!__recalc && (std::isinf(__ac) || std::isinf(__bd) || |
| std::isinf(__ad) || std::isinf(__bc))) { |
| if (std::isnan(__a)) |
| __a = copysign(_Tp(0), __a); |
| if (std::isnan(__b)) |
| __b = copysign(_Tp(0), __b); |
| if (std::isnan(__c)) |
| __c = copysign(_Tp(0), __c); |
| if (std::isnan(__d)) |
| __d = copysign(_Tp(0), __d); |
| __recalc = true; |
| } |
| if (__recalc) { |
| __x = _Tp(INFINITY) * (__a * __c - __b * __d); |
| __y = _Tp(INFINITY) * (__a * __d + __b * __c); |
| } |
| } |
| return complex<_Tp>(__x, __y); |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator*(const complex<_Tp> &__x, const _Tp &__y) { |
| complex<_Tp> __t(__x); |
| __t *= __y; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator*(const _Tp &__x, const complex<_Tp> &__y) { |
| complex<_Tp> __t(__y); |
| __t *= __x; |
| return __t; |
| } |
| |
| template <class _Tp> |
| complex<_Tp> |
| operator/(const complex<_Tp> &__z, const complex<_Tp> &__w) { |
| int __ilogbw = 0; |
| _Tp __a = __z.real(); |
| _Tp __b = __z.imag(); |
| _Tp __c = __w.real(); |
| _Tp __d = __w.imag(); |
| _Tp __logbw = logb(fmax(fabs(__c), fabs(__d))); |
| if (std::isfinite(__logbw)) { |
| __ilogbw = static_cast<int>(__logbw); |
| __c = scalbn(__c, -__ilogbw); |
| __d = scalbn(__d, -__ilogbw); |
| } |
| _Tp __denom = __c * __c + __d * __d; |
| _Tp __x = scalbn((__a * __c + __b * __d) / __denom, -__ilogbw); |
| _Tp __y = scalbn((__b * __c - __a * __d) / __denom, -__ilogbw); |
| if (std::isnan(__x) && std::isnan(__y)) { |
| if ((__denom == _Tp(0)) && (!std::isnan(__a) || !std::isnan(__b))) { |
| __x = copysign(_Tp(INFINITY), __c) * __a; |
| __y = copysign(_Tp(INFINITY), __c) * __b; |
| } else if ((std::isinf(__a) || std::isinf(__b)) && std::isfinite(__c) && std::isfinite(__d)) { |
| __a = copysign(std::isinf(__a) ? _Tp(1) : _Tp(0), __a); |
| __b = copysign(std::isinf(__b) ? _Tp(1) : _Tp(0), __b); |
| __x = _Tp(INFINITY) * (__a * __c + __b * __d); |
| __y = _Tp(INFINITY) * (__b * __c - __a * __d); |
| } else if (std::isinf(__logbw) && __logbw > _Tp(0) && std::isfinite(__a) && std::isfinite(__b)) { |
| __c = copysign(std::isinf(__c) ? _Tp(1) : _Tp(0), __c); |
| __d = copysign(std::isinf(__d) ? _Tp(1) : _Tp(0), __d); |
| __x = _Tp(0) * (__a * __c + __b * __d); |
| __y = _Tp(0) * (__b * __c - __a * __d); |
| } |
| } |
| return complex<_Tp>(__x, __y); |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator/(const complex<_Tp> &__x, const _Tp &__y) { |
| return complex<_Tp>(__x.real() / __y, __x.imag() / __y); |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator/(const _Tp &__x, const complex<_Tp> &__y) { |
| complex<_Tp> __t(__x); |
| __t /= __y; |
| return __t; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator+(const complex<_Tp> &__x) { |
| return __x; |
| } |
| |
| template <class _Tp> |
| inline complex<_Tp> |
| operator-(const complex<_Tp> &__x) { |
| return complex<_Tp>(-__x.real(), -__x.imag()); |
| } |
| |
| template <class _Tp> |
| inline bool |
| operator==(const complex<_Tp> &__x, const complex<_Tp> &__y) { |
| return __x.real() == __y.real() && __x.imag() == __y.imag(); |
| } |
| |
| template <class _Tp> |
| inline bool |
| operator==(const complex<_Tp> &__x, const _Tp &__y) { |
| return __x.real() == __y && __x.imag() == 0; |
| } |
| |
| template <class _Tp> |
| inline bool |
| operator==(const _Tp &__x, const complex<_Tp> &__y) { |
| return __x == __y.real() && 0 == __y.imag(); |
| } |
| |
| template <class _Tp> |
| inline bool |
| operator!=(const complex<_Tp> &__x, const complex<_Tp> &__y) { |
| return !(__x == __y); |
| } |
| |
| template <class _Tp> |
| inline bool |
| operator!=(const complex<_Tp> &__x, const _Tp &__y) { |
| return !(__x == __y); |
| } |
| |
| template <class _Tp> |
| inline bool |
| operator!=(const _Tp &__x, const complex<_Tp> &__y) { |
| return !(__x == __y); |
| } |
| |
| template <class _Tp> _Tp abs(const std::complex<_Tp> &__c); |
| |
| // arg |
| |
| template <class _Tp> _Tp arg(const std::complex<_Tp> &__c); |
| |
| // norm |
| |
| template <class _Tp> _Tp norm(const std::complex<_Tp> &__c); |
| |
| // conj |
| |
| template <class _Tp> std::complex<_Tp> conj(const std::complex<_Tp> &__c); |
| |
| // proj |
| |
| template <class _Tp> std::complex<_Tp> proj(const std::complex<_Tp> &__c); |
| |
| // polar |
| |
| template <class _Tp> |
| complex<_Tp> polar(const _Tp &__rho, const _Tp &__theta = _Tp()); |
| |
| // log |
| |
| template <class _Tp> std::complex<_Tp> log(const std::complex<_Tp> &__x); |
| |
| // log10 |
| |
| template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x); |
| |
| // sqrt |
| |
| template <class _Tp> |
| std::complex<_Tp> sqrt(const std::complex<_Tp> &__x); |
| |
| // exp |
| |
| template <class _Tp> |
| std::complex<_Tp> exp(const std::complex<_Tp> &__x); |
| |
| // pow |
| |
| template <class _Tp> |
| std::complex<_Tp> pow(const std::complex<_Tp> &__x, |
| const std::complex<_Tp> &__y); |
| |
| // __sqr, computes pow(x, 2) |
| |
| template <class _Tp> std::complex<_Tp> __sqr(const std::complex<_Tp> &__x); |
| |
| // asinh |
| |
| template <class _Tp> |
| std::complex<_Tp> asinh(const std::complex<_Tp> &__x); |
| |
| // acosh |
| |
| template <class _Tp> |
| std::complex<_Tp> acosh(const std::complex<_Tp> &__x); |
| |
| // atanh |
| |
| template <class _Tp> |
| std::complex<_Tp> atanh(const std::complex<_Tp> &__x); |
| |
| // sinh |
| |
| template <class _Tp> |
| std::complex<_Tp> sinh(const std::complex<_Tp> &__x); |
| |
| // cosh |
| |
| template <class _Tp> |
| std::complex<_Tp> cosh(const std::complex<_Tp> &__x); |
| |
| // tanh |
| |
| template <class _Tp> |
| std::complex<_Tp> tanh(const std::complex<_Tp> &__x); |
| |
| // asin |
| |
| template <class _Tp> |
| std::complex<_Tp> asin(const std::complex<_Tp> &__x); |
| |
| // acos |
| |
| template <class _Tp> |
| std::complex<_Tp> acos(const std::complex<_Tp> &__x); |
| |
| // atan |
| |
| template <class _Tp> |
| std::complex<_Tp> atan(const std::complex<_Tp> &__x); |
| |
| // sin |
| |
| template <class _Tp> |
| std::complex<_Tp> sin(const std::complex<_Tp> &__x); |
| |
| // cos |
| |
| template <class _Tp> std::complex<_Tp> cos(const std::complex<_Tp> &__x); |
| |
| // tan |
| |
| template <class _Tp> |
| std::complex<_Tp> tan(const std::complex<_Tp> &__x); |
| |
| } // namespace std |