413 lines · plain
1#pragma once2 3#include <cmath>4 5#define INFINITY (__builtin_inff())6#define NAN (__builtin_nanf (""))7 8namespace std {9 10// Taken from libc++11template <class _Tp>12class complex {13public:14 typedef _Tp value_type;15 16private:17 value_type __re_;18 value_type __im_;19 20public:21 complex(const value_type &__re = value_type(), const value_type &__im = value_type())22 : __re_(__re), __im_(__im) {}23 template <class _Xp>24 complex(const complex<_Xp> &__c)25 : __re_(__c.real()), __im_(__c.imag()) {}26 27 value_type real() const { return __re_; }28 value_type imag() const { return __im_; }29 30 void real(value_type __re) { __re_ = __re; }31 void imag(value_type __im) { __im_ = __im; }32 33 complex &operator=(const value_type &__re) {34 __re_ = __re;35 __im_ = value_type();36 return *this;37 }38 complex &operator+=(const value_type &__re) {39 __re_ += __re;40 return *this;41 }42 complex &operator-=(const value_type &__re) {43 __re_ -= __re;44 return *this;45 }46 complex &operator*=(const value_type &__re) {47 __re_ *= __re;48 __im_ *= __re;49 return *this;50 }51 complex &operator/=(const value_type &__re) {52 __re_ /= __re;53 __im_ /= __re;54 return *this;55 }56 57 template <class _Xp>58 complex &operator=(const complex<_Xp> &__c) {59 __re_ = __c.real();60 __im_ = __c.imag();61 return *this;62 }63 template <class _Xp>64 complex &operator+=(const complex<_Xp> &__c) {65 __re_ += __c.real();66 __im_ += __c.imag();67 return *this;68 }69 template <class _Xp>70 complex &operator-=(const complex<_Xp> &__c) {71 __re_ -= __c.real();72 __im_ -= __c.imag();73 return *this;74 }75 template <class _Xp>76 complex &operator*=(const complex<_Xp> &__c) {77 *this = *this * complex(__c.real(), __c.imag());78 return *this;79 }80 template <class _Xp>81 complex &operator/=(const complex<_Xp> &__c) {82 *this = *this / complex(__c.real(), __c.imag());83 return *this;84 }85};86 87template <class _Tp>88inline complex<_Tp>89operator+(const complex<_Tp> &__x, const complex<_Tp> &__y) {90 complex<_Tp> __t(__x);91 __t += __y;92 return __t;93}94 95template <class _Tp>96inline complex<_Tp>97operator+(const complex<_Tp> &__x, const _Tp &__y) {98 complex<_Tp> __t(__x);99 __t += __y;100 return __t;101}102 103template <class _Tp>104inline complex<_Tp>105operator+(const _Tp &__x, const complex<_Tp> &__y) {106 complex<_Tp> __t(__y);107 __t += __x;108 return __t;109}110 111template <class _Tp>112inline complex<_Tp>113operator-(const complex<_Tp> &__x, const complex<_Tp> &__y) {114 complex<_Tp> __t(__x);115 __t -= __y;116 return __t;117}118 119template <class _Tp>120inline complex<_Tp>121operator-(const complex<_Tp> &__x, const _Tp &__y) {122 complex<_Tp> __t(__x);123 __t -= __y;124 return __t;125}126 127template <class _Tp>128inline complex<_Tp>129operator-(const _Tp &__x, const complex<_Tp> &__y) {130 complex<_Tp> __t(-__y);131 __t += __x;132 return __t;133}134 135template <class _Tp>136complex<_Tp>137operator*(const complex<_Tp> &__z, const complex<_Tp> &__w) {138 _Tp __a = __z.real();139 _Tp __b = __z.imag();140 _Tp __c = __w.real();141 _Tp __d = __w.imag();142 _Tp __ac = __a * __c;143 _Tp __bd = __b * __d;144 _Tp __ad = __a * __d;145 _Tp __bc = __b * __c;146 _Tp __x = __ac - __bd;147 _Tp __y = __ad + __bc;148 if (std::isnan(__x) && std::isnan(__y)) {149 bool __recalc = false;150 if (std::isinf(__a) || std::isinf(__b)) {151 __a = copysign(std::isinf(__a) ? _Tp(1) : _Tp(0), __a);152 __b = copysign(std::isinf(__b) ? _Tp(1) : _Tp(0), __b);153 if (std::isnan(__c))154 __c = copysign(_Tp(0), __c);155 if (std::isnan(__d))156 __d = copysign(_Tp(0), __d);157 __recalc = true;158 }159 if (std::isinf(__c) || std::isinf(__d)) {160 __c = copysign(std::isinf(__c) ? _Tp(1) : _Tp(0), __c);161 __d = copysign(std::isinf(__d) ? _Tp(1) : _Tp(0), __d);162 if (std::isnan(__a))163 __a = copysign(_Tp(0), __a);164 if (std::isnan(__b))165 __b = copysign(_Tp(0), __b);166 __recalc = true;167 }168 if (!__recalc && (std::isinf(__ac) || std::isinf(__bd) ||169 std::isinf(__ad) || std::isinf(__bc))) {170 if (std::isnan(__a))171 __a = copysign(_Tp(0), __a);172 if (std::isnan(__b))173 __b = copysign(_Tp(0), __b);174 if (std::isnan(__c))175 __c = copysign(_Tp(0), __c);176 if (std::isnan(__d))177 __d = copysign(_Tp(0), __d);178 __recalc = true;179 }180 if (__recalc) {181 __x = _Tp(INFINITY) * (__a * __c - __b * __d);182 __y = _Tp(INFINITY) * (__a * __d + __b * __c);183 }184 }185 return complex<_Tp>(__x, __y);186}187 188template <class _Tp>189inline complex<_Tp>190operator*(const complex<_Tp> &__x, const _Tp &__y) {191 complex<_Tp> __t(__x);192 __t *= __y;193 return __t;194}195 196template <class _Tp>197inline complex<_Tp>198operator*(const _Tp &__x, const complex<_Tp> &__y) {199 complex<_Tp> __t(__y);200 __t *= __x;201 return __t;202}203 204template <class _Tp>205complex<_Tp>206operator/(const complex<_Tp> &__z, const complex<_Tp> &__w) {207 int __ilogbw = 0;208 _Tp __a = __z.real();209 _Tp __b = __z.imag();210 _Tp __c = __w.real();211 _Tp __d = __w.imag();212 _Tp __logbw = logb(fmax(fabs(__c), fabs(__d)));213 if (std::isfinite(__logbw)) {214 __ilogbw = static_cast<int>(__logbw);215 __c = scalbn(__c, -__ilogbw);216 __d = scalbn(__d, -__ilogbw);217 }218 _Tp __denom = __c * __c + __d * __d;219 _Tp __x = scalbn((__a * __c + __b * __d) / __denom, -__ilogbw);220 _Tp __y = scalbn((__b * __c - __a * __d) / __denom, -__ilogbw);221 if (std::isnan(__x) && std::isnan(__y)) {222 if ((__denom == _Tp(0)) && (!std::isnan(__a) || !std::isnan(__b))) {223 __x = copysign(_Tp(INFINITY), __c) * __a;224 __y = copysign(_Tp(INFINITY), __c) * __b;225 } else if ((std::isinf(__a) || std::isinf(__b)) && std::isfinite(__c) && std::isfinite(__d)) {226 __a = copysign(std::isinf(__a) ? _Tp(1) : _Tp(0), __a);227 __b = copysign(std::isinf(__b) ? _Tp(1) : _Tp(0), __b);228 __x = _Tp(INFINITY) * (__a * __c + __b * __d);229 __y = _Tp(INFINITY) * (__b * __c - __a * __d);230 } else if (std::isinf(__logbw) && __logbw > _Tp(0) && std::isfinite(__a) && std::isfinite(__b)) {231 __c = copysign(std::isinf(__c) ? _Tp(1) : _Tp(0), __c);232 __d = copysign(std::isinf(__d) ? _Tp(1) : _Tp(0), __d);233 __x = _Tp(0) * (__a * __c + __b * __d);234 __y = _Tp(0) * (__b * __c - __a * __d);235 }236 }237 return complex<_Tp>(__x, __y);238}239 240template <class _Tp>241inline complex<_Tp>242operator/(const complex<_Tp> &__x, const _Tp &__y) {243 return complex<_Tp>(__x.real() / __y, __x.imag() / __y);244}245 246template <class _Tp>247inline complex<_Tp>248operator/(const _Tp &__x, const complex<_Tp> &__y) {249 complex<_Tp> __t(__x);250 __t /= __y;251 return __t;252}253 254template <class _Tp>255inline complex<_Tp>256operator+(const complex<_Tp> &__x) {257 return __x;258}259 260template <class _Tp>261inline complex<_Tp>262operator-(const complex<_Tp> &__x) {263 return complex<_Tp>(-__x.real(), -__x.imag());264}265 266template <class _Tp>267inline bool268operator==(const complex<_Tp> &__x, const complex<_Tp> &__y) {269 return __x.real() == __y.real() && __x.imag() == __y.imag();270}271 272template <class _Tp>273inline bool274operator==(const complex<_Tp> &__x, const _Tp &__y) {275 return __x.real() == __y && __x.imag() == 0;276}277 278template <class _Tp>279inline bool280operator==(const _Tp &__x, const complex<_Tp> &__y) {281 return __x == __y.real() && 0 == __y.imag();282}283 284template <class _Tp>285inline bool286operator!=(const complex<_Tp> &__x, const complex<_Tp> &__y) {287 return !(__x == __y);288}289 290template <class _Tp>291inline bool292operator!=(const complex<_Tp> &__x, const _Tp &__y) {293 return !(__x == __y);294}295 296template <class _Tp>297inline bool298operator!=(const _Tp &__x, const complex<_Tp> &__y) {299 return !(__x == __y);300}301 302template <class _Tp> _Tp abs(const std::complex<_Tp> &__c);303 304// arg305 306template <class _Tp> _Tp arg(const std::complex<_Tp> &__c);307 308// norm309 310template <class _Tp> _Tp norm(const std::complex<_Tp> &__c);311 312// conj313 314template <class _Tp> std::complex<_Tp> conj(const std::complex<_Tp> &__c);315 316// proj317 318template <class _Tp> std::complex<_Tp> proj(const std::complex<_Tp> &__c);319 320// polar321 322template <class _Tp>323complex<_Tp> polar(const _Tp &__rho, const _Tp &__theta = _Tp());324 325// log326 327template <class _Tp> std::complex<_Tp> log(const std::complex<_Tp> &__x);328 329// log10330 331template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x);332 333// sqrt334 335template <class _Tp>336std::complex<_Tp> sqrt(const std::complex<_Tp> &__x);337 338// exp339 340template <class _Tp>341std::complex<_Tp> exp(const std::complex<_Tp> &__x);342 343// pow344 345template <class _Tp>346std::complex<_Tp> pow(const std::complex<_Tp> &__x,347 const std::complex<_Tp> &__y);348 349// __sqr, computes pow(x, 2)350 351template <class _Tp> std::complex<_Tp> __sqr(const std::complex<_Tp> &__x);352 353// asinh354 355template <class _Tp>356std::complex<_Tp> asinh(const std::complex<_Tp> &__x);357 358// acosh359 360template <class _Tp>361std::complex<_Tp> acosh(const std::complex<_Tp> &__x);362 363// atanh364 365template <class _Tp>366std::complex<_Tp> atanh(const std::complex<_Tp> &__x);367 368// sinh369 370template <class _Tp>371std::complex<_Tp> sinh(const std::complex<_Tp> &__x);372 373// cosh374 375template <class _Tp>376std::complex<_Tp> cosh(const std::complex<_Tp> &__x);377 378// tanh379 380template <class _Tp>381std::complex<_Tp> tanh(const std::complex<_Tp> &__x);382 383// asin384 385template <class _Tp>386std::complex<_Tp> asin(const std::complex<_Tp> &__x);387 388// acos389 390template <class _Tp>391std::complex<_Tp> acos(const std::complex<_Tp> &__x);392 393// atan394 395template <class _Tp>396std::complex<_Tp> atan(const std::complex<_Tp> &__x);397 398// sin399 400template <class _Tp>401std::complex<_Tp> sin(const std::complex<_Tp> &__x);402 403// cos404 405template <class _Tp> std::complex<_Tp> cos(const std::complex<_Tp> &__x);406 407// tan408 409template <class _Tp>410std::complex<_Tp> tan(const std::complex<_Tp> &__x);411 412} // namespace std413