#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