1#pragma once 2 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 bool 268operator==(const complex<_Tp> &__x, const complex<_Tp> &__y) { 269 return __x.real() == __y.real() && __x.imag() == __y.imag(); 270} 271 272template <class _Tp> 273inline bool 274operator==(const complex<_Tp> &__x, const _Tp &__y) { 275 return __x.real() == __y && __x.imag() == 0; 276} 277 278template <class _Tp> 279inline bool 280operator==(const _Tp &__x, const complex<_Tp> &__y) { 281 return __x == __y.real() && 0 == __y.imag(); 282} 283 284template <class _Tp> 285inline bool 286operator!=(const complex<_Tp> &__x, const complex<_Tp> &__y) { 287 return !(__x == __y); 288} 289 290template <class _Tp> 291inline bool 292operator!=(const complex<_Tp> &__x, const _Tp &__y) { 293 return !(__x == __y); 294} 295 296template <class _Tp> 297inline bool 298operator!=(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// arg 305 306template <class _Tp> _Tp arg(const std::complex<_Tp> &__c); 307 308// norm 309 310template <class _Tp> _Tp norm(const std::complex<_Tp> &__c); 311 312// conj 313 314template <class _Tp> std::complex<_Tp> conj(const std::complex<_Tp> &__c); 315 316// proj 317 318template <class _Tp> std::complex<_Tp> proj(const std::complex<_Tp> &__c); 319 320// polar 321 322template <class _Tp> 323complex<_Tp> polar(const _Tp &__rho, const _Tp &__theta = _Tp()); 324 325// log 326 327template <class _Tp> std::complex<_Tp> log(const std::complex<_Tp> &__x); 328 329// log10 330 331template <class _Tp> std::complex<_Tp> log10(const std::complex<_Tp> &__x); 332 333// sqrt 334 335template <class _Tp> 336std::complex<_Tp> sqrt(const std::complex<_Tp> &__x); 337 338// exp 339 340template <class _Tp> 341std::complex<_Tp> exp(const std::complex<_Tp> &__x); 342 343// pow 344 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// asinh 354 355template <class _Tp> 356std::complex<_Tp> asinh(const std::complex<_Tp> &__x); 357 358// acosh 359 360template <class _Tp> 361std::complex<_Tp> acosh(const std::complex<_Tp> &__x); 362 363// atanh 364 365template <class _Tp> 366std::complex<_Tp> atanh(const std::complex<_Tp> &__x); 367 368// sinh 369 370template <class _Tp> 371std::complex<_Tp> sinh(const std::complex<_Tp> &__x); 372 373// cosh 374 375template <class _Tp> 376std::complex<_Tp> cosh(const std::complex<_Tp> &__x); 377 378// tanh 379 380template <class _Tp> 381std::complex<_Tp> tanh(const std::complex<_Tp> &__x); 382 383// asin 384 385template <class _Tp> 386std::complex<_Tp> asin(const std::complex<_Tp> &__x); 387 388// acos 389 390template <class _Tp> 391std::complex<_Tp> acos(const std::complex<_Tp> &__x); 392 393// atan 394 395template <class _Tp> 396std::complex<_Tp> atan(const std::complex<_Tp> &__x); 397 398// sin 399 400template <class _Tp> 401std::complex<_Tp> sin(const std::complex<_Tp> &__x); 402 403// cos 404 405template <class _Tp> std::complex<_Tp> cos(const std::complex<_Tp> &__x); 406 407// tan 408 409template <class _Tp> 410std::complex<_Tp> tan(const std::complex<_Tp> &__x); 411 412} // namespace std 413