// xcomplex internal header
// NOTE: no include guard

		// TEMPLATE FUNCTION operator+
_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator+(const _CMPLX(_Ty)& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// add complex to complex
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp += _Right);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator+(const _CMPLX(_Ty)& _Left,
		const _Ty& _Right)
	{	// add real to complex
	_CMPLX(_Ty) _Tmp(_Left);
	_Tmp.real(_Tmp.real() + _Right);
	return (_Tmp);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator+(const _Ty& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// add complex to real
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp += _Right);
	}

		// TEMPLATE FUNCTION operator-
_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator-(const _CMPLX(_Ty)& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// subtract complex from complex
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp -= _Right);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator-(const _CMPLX(_Ty)& _Left,
		const _Ty& _Right)
	{	// subtract real from complex
	_CMPLX(_Ty) _Tmp(_Left);
	_Tmp.real(_Tmp.real() - _Right);
	return (_Tmp);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator-(const _Ty& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// subtract complex from real
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp -= _Right);
	}

		// TEMPLATE FUNCTION operator*
_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator*(const _CMPLX(_Ty)& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// multiply complex by complex
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp *= _Right);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator*(const _CMPLX(_Ty)& _Left,
		const _Ty& _Right)
	{	// multiply complex by real
	_CMPLX(_Ty) _Tmp(_Left);
	_Tmp.real(_Tmp.real() * _Right);
	_Tmp.imag(_Tmp.imag() * _Right);
	return (_Tmp);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator*(const _Ty& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// multiply real by complex
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp *= _Right);
	}

		// TEMPLATE FUNCTION operator/
_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator/(const _CMPLX(_Ty)& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// divide complex by complex
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp /= _Right);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator/(const _CMPLX(_Ty)& _Left,
		const _Ty& _Right)
	{	// divide complex by real
	_CMPLX(_Ty) _Tmp(_Left);
	_Tmp.real(_Tmp.real() / _Right);
	_Tmp.imag(_Tmp.imag() / _Right);
	return (_Tmp);
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator/(const _Ty& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// divide real by complex
	_CMPLX(_Ty) _Tmp(_Left);
	return (_Tmp /= _Right);
	}

		// TEMPLATE FUNCTION UNARY operator+
_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator+(const _CMPLX(_Ty)& _Left)
	{	// return +complex
	return (_CMPLX(_Ty)(_Left));
	}

		// TEMPLATE FUNCTION UNARY operator-
_TMPLT(_Ty) inline
	_CMPLX(_Ty) operator-(const _CMPLX(_Ty)& _Left)
	{	// return -complex
	return (_CMPLX(_Ty)(-_Left.real(), -_Left.imag()));
	}

		// TEMPLATE FUNCTION operator==
_TMPLT(_Ty) inline
	bool operator==(const _CMPLX(_Ty)& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// test complex equal to complex
	return (_Left.real() == _Right.real()
		&& _Left.imag() == _Right.imag());
	}

_TMPLT(_Ty) inline
	bool operator==(const _CMPLX(_Ty)& _Left,
		const _Ty& _Right)
	{	// test real equal to complex
	return (_Left.real() == _Right
		&& _Left.imag() == 0);
	}

_TMPLT(_Ty) inline
	bool operator==(const _Ty& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// test complex equal to real
	return (_Left == _Right.real()
		&& 0 == _Right.imag());
	}

		// TEMPLATE FUNCTION operator!=
_TMPLT(_Ty) inline
	bool operator!=(const _CMPLX(_Ty)& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// test complex not equal to complex
	return (!(_Left == _Right));
	}

_TMPLT(_Ty) inline
	bool operator!=(const _CMPLX(_Ty)& _Left,
		const _Ty& _Right)
	{	// test real not equal to complex
	return (!(_Left == _Right));
	}

_TMPLT(_Ty) inline
	bool operator!=(const _Ty& _Left,
		const _CMPLX(_Ty)& _Right)
	{	// test complex not equal to real
	return (!(_Left == _Right));
	}

		// TEMPLATE FUNCTION imag
_TMPLT(_Ty) inline
	_Ty imag(const _CMPLX(_Ty)& _Left)
	{	// return imaginary component
	return (_Left.imag());
	}

		// TEMPLATE FUNCTION real
_TMPLT(_Ty) inline
	_Ty real(const _CMPLX(_Ty)& _Left)
	{	// return real component
	return (_Left.real());
	}

		// TEMPLATE FUNCTION _Fabs
_TMPLT(_Ty) inline
	_Ty _Fabs(const _CMPLX(_Ty)& _Left, int *_Pexp)
	{	// return magnitude and scale factor
	*_Pexp = 0;
	_Ty _Av = real(_Left);
	_Ty _Bv = imag(_Left);

	if (_CTR(_Ty)::_Isinf(_Av) || _CTR(_Ty)::_Isinf(_Bv))
		return (_CTR(_Ty)::_Infv(_Bv));	// at least one component is INF
	else if (_CTR(_Ty)::_Isnan(_Av))
		return (_Av);	//	real component is NaN
	else if (_CTR(_Ty)::_Isnan(_Bv))
		return (_Bv);	// imaginary component is NaN
	else
		{	// neither component is NaN or INF
		if (_Av < 0)
			_Av = -_Av;
		if (_Bv < 0)
			_Bv = -_Bv;
		if (_Av < _Bv)
			{	// ensure that |_Bv| <= |_Av|
			_Ty _Tmp = _Av;
			_Av = _Bv;
			_Bv = _Tmp;
			}

		if (_Av == 0)
			return (_Av);	// |0| == 0
		if (1 <= _Av)
			*_Pexp = 2, _Av = _Av * (_Ty)0.25, _Bv = _Bv * (_Ty)0.25;
		else
			*_Pexp = -2, _Av = _Av * 4, _Bv = _Bv * 4;

		_Ty _Tmp = _Av - _Bv;
		if (_Tmp == _Av)
			return (_Av);	// _Bv unimportant
		else if (_Bv < _Tmp)
			{	// use simple approximation
			const _Ty _Qv = _Av / _Bv;
			return (_Av + _Bv / (_Qv + _CTR(_Ty)::sqrt(_Qv * _Qv + 1)));
			}
		else
			{	// use 1 1/2 precision to preserve bits
			static const _Ty _Root2 =
				(_Ty)1.4142135623730950488016887242096981L;
			static const _Ty _Oneplusroot2high =
				(_Ty)(10125945.0 / 4194304.0);	// exact if prec >= 24 bits
			static const _Ty _Oneplusroot2low =
				(_Ty)1.4341252375973918872420969807856967e-7L;

			const _Ty _Qv = _Tmp / _Bv;
			const _Ty _Rv = (_Qv + 2) * _Qv;
			const _Ty _Sv = _Rv / (_Root2 + _CTR(_Ty)::sqrt(_Rv + 2))
				+ _Oneplusroot2low + _Qv + _Oneplusroot2high;
			return (_Av + _Bv / _Sv);
			}
		}
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) sqrt(const _CMPLX(_Ty)& _Left);

		// TEMPLATE FUNCTION abs
_TMPLT(_Ty) inline
	_Ty abs(const _CMPLX(_Ty)& _Left)
	{	// return |complex| as real
	int _Leftexp;
	_Ty _Rho = _Fabs(_Left, &_Leftexp);	// get magnitude and scale factor

	if (_Leftexp == 0)
		return (_Rho);	// no scale factor
	else
		return (_CTR(_Ty)::ldexp(_Rho, _Leftexp));	// scale result
	}

		// TEMPLATE FUNCTION acos
_TMPLT(_Ty) inline
	_CMPLX(_Ty) acos(const _CMPLX(_Ty)& _Left)
	{	// return acos
	static const _Ty _Arcbig = (_Ty)0.25L
		* _CTR(_Ty)::sqrt(_CTR(_Ty)::_Flt_max());
	static const _Ty _Pi = (_Ty)3.1415926535897932384626433832795029L;

	const _Ty _Re = real(_Left);
	const _Ty _Im = imag(_Left);
	_Ty _Ux;
	_Ty _Vx;

	if (_CTR(_Ty)::_Isnan(_Re) || _CTR(_Ty)::_Isnan(_Im))
		{	// at least one NaN
		_Ux = _CTR(_Ty)::_Nanv(_Re);
		_Vx = _Ux;
		}
	else if (_CTR(_Ty)::_Isinf(_Re))
		{	// (+/-Inf, not NaN)
		if (_CTR(_Ty)::_Isinf(_Im))
			if (_Re < 0)
				_Ux = (_Ty)0.75 * _Pi;	// (-Inf, +/-Inf)
			else
				_Ux = (_Ty)0.25 * _Pi;	// (+Inf, +/-Inf)
		else
			if (_Re < 0)
				_Ux = _Pi;	// (-Inf, finite)
			else
				_Ux = 0;	// (+Inf, finite)
		_Vx = -_CTR(_Ty)::_Infv(_Re);
		if (_Im < 0)
			_Vx = -_Vx;
		}
	else if (_CTR(_Ty)::_Isinf(_Im))
		{	// (finite, finite)
		_Ux = (_Ty)0.50 * _Pi;	// (finite, +/-Inf)
		_Vx = -_Im;
		}
	else
		{	// (finite, finite)
		const _CMPLX(_Ty) _Wx = sqrt(_CMPLX(_Ty)(1 + _Re, -_Im));
		const _CMPLX(_Ty) _Zx = sqrt(_CMPLX(_Ty)(1 - _Re, -_Im));
		const _Ty _Wr = real(_Wx);
		const _Ty _Wi = imag(_Wx);
		const _Ty _Zr = real(_Zx);
		const _Ty _Zi = imag(_Zx);
		_Ty _Alfa, _Beta;

		_Ux = 2 * _CTR(_Ty)::atan2(_Zr, _Wr);

		if (_Arcbig < _Wr)
			{	// real parts large
			_Alfa = _Wr;
			_Beta = _Zi + _Wi * (_Zr / _Alfa);
			}
		else if (_Arcbig < _Wi)
			{	// imag parts large
			_Alfa = _Wi;
			_Beta = _Wr * (_Zi / _Alfa) + _Zr;
			}
		else if (_Wi < -_Arcbig)
			{	// imag part of w large negative
			_Alfa = -_Wi;
			_Beta = _Wr * (_Zi / _Alfa) + _Zr;
			}
		else
			{	// shouldn't overflow
			_Alfa = 0;
			_Beta = _Wr * _Zi + _Wi * _Zr;	// Im(w * z)
			}

		_Vx = _CTR(_Ty)::asinh(_Beta);
		if (_Alfa == 0)
			;
		else if (0 <= _Ux)
			_Vx += _CTR(_Ty)::log(_Alfa);
		else
			_Vx -= _CTR(_Ty)::log(_Alfa);	// asinh(a*b) = asinh(a)+log(b)
		}
	return (_CMPLX(_Ty)(_Ux, _Vx));
	}

		// TEMPLATE FUNCTION acosh
_TMPLT(_Ty) inline
	_CMPLX(_Ty) acosh(const _CMPLX(_Ty)& _Left)
	{	// return acosh
	static const _Ty _Arcbig = (_Ty)0.25L
		* _CTR(_Ty)::sqrt(_CTR(_Ty)::_Flt_max());
	static const _Ty _Pi = (_Ty)3.1415926535897932384626433832795029L;

	const _Ty _Re = real(_Left);
	_Ty _Im = imag(_Left);
	_Ty _Ux;
	_Ty _Vx;

	if (_CTR(_Ty)::_Isnan(_Re) || _CTR(_Ty)::_Isnan(_Im))
		{	// at least one NaN
		_Ux = _CTR(_Ty)::_Nanv(_Re);
		_Vx = _Ux;
		}
	else if (_CTR(_Ty)::_Isinf(_Re))
		{	// (+/-Inf, not NaN)
		_Ux = _CTR(_Ty)::_Infv(_Re);

		if (_CTR(_Ty)::_Isinf(_Im))
			if (_Re < 0)
				_Vx = (_Ty)0.75 * _Pi;	// (-Inf, +/-Inf)
			else
				_Vx = (_Ty)0.25 * _Pi;	// (+Inf, +/-Inf)
		else
			if (_Re < 0)
				_Vx = _Pi;	// (-Inf, finite)
			else
				_Vx = 0;	// (+Inf, finite)
		if (_Im < 0)
			_Vx = -_Vx;
		}
	else if (_CTR(_Ty)::_Isinf(_Im))
		{	// (finite, +/-Inf)
		_Ux = _CTR(_Ty)::_Infv(_Re);
		_Vx = _Im < 0 ? -(_Ty)0.50 * _Pi : (_Ty)0.50 * _Pi;
		}
	else
		{	// (finite, finite)
		const _CMPLX(_Ty) _Wx = sqrt(_CMPLX(_Ty)(_Re - 1, -_Im));
		const _CMPLX(_Ty) _Zx = sqrt(_CMPLX(_Ty)(_Re + 1, _Im));
		const _Ty _Wr = real(_Wx);
		const _Ty _Wi = imag(_Wx);
		const _Ty _Zr = real(_Zx);
		const _Ty _Zi = imag(_Zx);
		_Ty _Alfa, _Beta;

		if (_Arcbig < _Wr)
			{	// real parts large
			_Alfa = _Wr;
			_Beta = _Zr - _Wi * (_Zi / _Alfa);
			}
		else if (_Arcbig < _Wi)
			{	// imag parts large
			_Alfa = _Wi;
			_Beta = _Wr * (_Zr / _Alfa) - _Zi;
			}
		else if (_Wi < -_Arcbig)
			{	// imag part of w large negative
			_Alfa = -_Wi;
			_Beta = _Wr * (_Zr / _Alfa) - _Zi;
			}
		else
			{	// shouldn't overflow
			_Alfa = 0;
			_Beta = _Wr * _Zr - _Wi * _Zi;	// Re(w * z)
			}

		_Ux = _CTR(_Ty)::asinh(_Beta);
		if (_Alfa == 0)
			;
		else if (0 <= _Ux)
			_Ux += _CTR(_Ty)::log(_Alfa);
		else
			_Ux -= _CTR(_Ty)::log(_Alfa);	// asinh(a*b) = asinh(a)+log(b)

		bool _Neg = true;
		if (_Im < 0)
			_Im = -_Im;
		else
			_Neg = false;
		_Vx = 2 * _CTR(_Ty)::atan2(
			imag(sqrt(_CMPLX(_Ty)(_Re - 1, _Im))), _Zr);
		if (_Neg)
			_Vx = -_Vx;
		}
	return (_CMPLX(_Ty)(_Ux, _Vx));
	}

		// TEMPLATE FUNCTION asinh
_TMPLT(_Ty) inline
	_CMPLX(_Ty) asinh(const _CMPLX(_Ty)& _Left)
	{	// return asinh
	static const _Ty _Arcbig = (_Ty)0.25L
		* _CTR(_Ty)::sqrt(_CTR(_Ty)::_Flt_max());
	static const _Ty _Pi = (_Ty)3.1415926535897932384626433832795029L;

	const _Ty _Re = real(_Left);
	_Ty _Im = imag(_Left);
	_Ty _Ux;
	_Ty _Vx;

	if (_CTR(_Ty)::_Isnan(_Re) || _CTR(_Ty)::_Isnan(_Im))
		{	// at least one NaN/Inf
		_Ux = _CTR(_Ty)::_Nanv(_Re);
		_Vx = _Ux;
		}
	else if (_CTR(_Ty)::_Isinf(_Re))
		{	// (+/-Inf, not NaN)
		_Ux = _CTR(_Ty)::_Infv(_Re);

		if (_CTR(_Ty)::_Isinf(_Im))
			{	// (+/-Inf, +/-Inf)
			_Ux = _Re;
			_Vx = _Im < 0 ? -(_Ty)0.25 * _Pi : (_Ty)0.25 * _Pi;
			}
		else
			{	// (+/-Inf, finite)
			_Ux = _Re;
			_Vx = _Im < 0 ? -(_Ty)0 : (_Ty)0;
			}
		}
	else if (_CTR(_Ty)::_Isinf(_Im))
		{	// (finite, +/-Inf)
		_Ux = _Im;
		_Vx = _Im < 0 ? -(_Ty)0.50 * _Pi : (_Ty)0.50 * _Pi;
		}
	else
		{	// (finite, finite)
		const _CMPLX(_Ty) _Wx = sqrt(_CMPLX(_Ty)(1 - _Im, _Re));
		const _CMPLX(_Ty) _Zx = sqrt(_CMPLX(_Ty)(1 + _Im, -_Re));
		const _Ty _Wr = real(_Wx);
		const _Ty _Wi = imag(_Wx);
		const _Ty _Zr = real(_Zx);
		const _Ty _Zi = imag(_Zx);
		_Ty _Alfa, _Beta;

		if (_Arcbig < _Wr)
			{	// real parts large
			_Alfa = _Wr;
			_Beta = _Wi * (_Zr / _Alfa) - _Zi;
			}
		else if (_Arcbig < _Wi)
			{	// imag parts large
			_Alfa = _Wi;
			_Beta = _Zr - _Wr * (_Zi / _Alfa);
			}
		else if (_Wi < -_Arcbig)
			{	// imag part of w large negative
			_Alfa = -_Wi;
			_Beta = _Zr - _Wr * (_Zi / _Alfa);
			}
		else
			{	// shouldn't overflow
			_Alfa = 0;
			_Beta = _Wi * _Zr - _Wr * _Zi;	// Im(w * conj(z))
			}

		_Ux = _CTR(_Ty)::asinh(_Beta);
		if (_Alfa == 0)
			;
		else if (0 <= _Ux)
			_Ux += _CTR(_Ty)::log(_Alfa);
		else
			_Ux -= _CTR(_Ty)::log(_Alfa);	// asinh(a*b) = asinh(a)+log(b)

		_Vx = _CTR(_Ty)::atan2(_Im, real(_Wx * _Zx));
		}

	return (_CMPLX(_Ty)(_Ux, _Vx));
	}

		// TEMPLATE FUNCTION asin
_TMPLT(_Ty) inline
	_CMPLX(_Ty) asin(const _CMPLX(_Ty)& _Left)
	{	// return asin
	_CMPLX(_Ty) _Asinh = _STD asinh(_CMPLX(_Ty)(-imag(_Left), real(_Left)));

	return (_CMPLX(_Ty)(imag(_Asinh), -real(_Asinh)));
	}

		// TEMPLATE FUNCTION atanh
_TMPLT(_Ty) inline
	_CMPLX(_Ty) atanh(const _CMPLX(_Ty)& _Left)
	{	// return atanh
	static const _Ty _Arcbig = (_Ty)0.25L
		* _CTR(_Ty)::sqrt(_CTR(_Ty)::_Flt_max());
	static const _Ty _Piby2 = (_Ty)1.5707963267948966192313216916397514L;

	_Ty _Re = real(_Left);
	_Ty _Im = imag(_Left);
	_Ty _Ux;
	_Ty _Vx;

	if (_CTR(_Ty)::_Isnan(_Re) || _CTR(_Ty)::_Isnan(_Im))
		{	// at least one NaN
		_Ux = _CTR(_Ty)::_Nanv(_Re);
		_Vx = _Ux;
		}
	else if (_CTR(_Ty)::_Isinf(_Re))
		{	// (+/-Inf, not NaN)
		_Ux = _Re < 0 ? -(_Ty)0 : (_Ty)0;
		_Vx = _Im < 0 ? -_Piby2 : _Piby2;
		}
	else
		{	// (finite, not _NaN)
		_Ty _Magim = _Im < 0 ? -_Im : _Im;
		bool _Neg = _Re < 0;

		if (_Neg)
			_Re = -_Re;
		else
			_Im = -_Im;

		if (_Arcbig < _Re)
			{	// |re| is large
			_Ty _Fx = _Im / _Re;

			_Ux = 1 / _Re / (1 + _Fx * _Fx);
			_Vx = _Im < 0 ?  -_Piby2 : _Piby2;
			}
		else if (_Arcbig < _Magim)
			{	// |im| is large
			_Ty _Fx = _Re / _Im;

			_Ux = _Fx / _Im / (1 + _Fx * _Fx);
			_Vx = _Im < 0 ?  -_Piby2 : _Piby2;
			}
		else if (_Re != 1)
			{	// |re| is small
			_Ty _Refrom1 = 1 - _Re;
			_Ty _Imeps2 = _Magim * _Magim;

			_Ux = (_Ty)0.25 * _CTR(_Ty)::log1p(4 * _Re
				/ (_Refrom1 * _Refrom1 + _Imeps2));
			_Vx = (_Ty)0.50 * _CTR(_Ty)::atan2(2 * _Im,
				_Refrom1 * (1 + _Re) - _Imeps2);
			}
		else if (_Im == 0)
			{	// {+/-1, 0)
			_Ux = _CTR(_Ty)::_Infv(_Re);
			_Vx = _Im;
			}
		else
			{	// {+/-1, nonzero)
			_Ux = _CTR(_Ty)::log(
				_CTR(_Ty)::sqrt(_CTR(_Ty)::sqrt(4 + _Im * _Im))
					/ _CTR(_Ty)::sqrt(_Magim));
			_Vx = (_Ty)0.50 * (_Piby2 + _CTR(_Ty)::atan2(_Magim, 2));
			if (_Im < 0)
				_Vx = -_Vx;
			}

		if (_Neg)
			_Ux = -_Ux;
		else
			_Vx = -_Vx;
		}
	return (_CMPLX(_Ty)(_Ux, _Vx));
	}

		// TEMPLATE FUNCTION atan
_TMPLT(_Ty) inline
	_CMPLX(_Ty) atan(const _CMPLX(_Ty)& _Left)
	{	// return atan
	_CMPLX(_Ty) _Atanh = _STD atanh(_CMPLX(_Ty)(-imag(_Left), real(_Left)));

	return (_CMPLX(_Ty)(imag(_Atanh), -real(_Atanh)));
	}

		// TEMPLATE FUNCTION cosh
_TMPLT(_Ty) inline
	_CMPLX(_Ty) cosh(const _CMPLX(_Ty)& _Left)
	{	// return cosh(complex)
	return (_CMPLX(_Ty)(
		_CTR(_Ty)::_Cosh(real(_Left), _CTR(_Ty)::cos(imag(_Left))),
		_CTR(_Ty)::_Sinh(real(_Left), _CTR(_Ty)::sin(imag(_Left)))));
	}

		// TEMPLATE FUNCTION exp
_TMPLT(_Ty) inline
	_CMPLX(_Ty) exp(const _CMPLX(_Ty)& _Left)
	{	// return exp(complex)
	_Ty _Real(real(_Left)), _Imag(real(_Left));
	_CTR(_Ty)::_Exp(&_Real, _CTR(_Ty)::cos(imag(_Left)), 0);
	_CTR(_Ty)::_Exp(&_Imag, _CTR(_Ty)::sin(imag(_Left)), 0);
	return (_CMPLX(_Ty)(_Real, _Imag));
	}

		// TEMPLATE FUNCTION log
_TMPLT(_Ty) inline
	_CMPLX(_Ty) log(const _CMPLX(_Ty)& _Left)
	{	// return log(complex)
	_Ty _Theta = _CTR(_Ty)::atan2(imag(_Left), real(_Left));	// get phase

	if (_CTR(_Ty)::_Isnan(_Theta))
		return (_CMPLX(_Ty)(_Theta, _Theta));	// real or imag is NaN
	else
		{	// use 1 1/2 precision to preserve bits
		static const _Ty _Cm = (_Ty)(22713.0L / 32768.0L);
		static const _Ty _Cl =
			(_Ty)1.4286068203094172321214581765680755e-6L;
		int _Leftexp;
		_Ty _Rho = _Fabs(_Left, &_Leftexp);	// get magnitude and scale factor

		_Ty _Leftn = (_Ty)_Leftexp;
		_CMPLX(_Ty) _Tmp(
			_Rho == 0 ? -_CTR(_Ty)::_Infv(_Rho)	// log(0) == -INF
			: _CTR(_Ty)::_Isinf(_Rho) ? _Rho	// log(INF) == INF
			: _CTR(_Ty)::log(_Rho) + _Leftn * _Cl + _Leftn * _Cm,
				_Theta);
		return (_Tmp);
		}
	}

		// TEMPLATE FUNCTION pow

_TMPLT(_Ty) inline
	_CMPLX(_Ty) pow(const _CMPLX(_Ty)& _Left, int _Right)
	{	// return complex ^ integer
	_CMPLX(_Ty) _Tmp = _Left;
	unsigned int _Count = _Right;

	if (_Right < 0)
		_Count = 0 - _Count;	// safe negation as unsigned

	for (_CMPLX(_Ty) _Zv = _CMPLX(_Ty)(1); ; _Tmp *= _Tmp)
		{	// fold in _Left ^ (2 ^ _Count) as needed
		if ((_Count & 1) != 0)
			_Zv *= _Tmp;
		if ((_Count >>= 1) == 0)
			return (_Right < 0 ? _CMPLX(_Ty)(1) / _Zv : _Zv);
		}
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) _Pow(const _Ty& _Left, const _Ty& _Right)
	{	// return real ^ real
	if (0 <= _Left)
		return (_CTR(_Ty)::pow(_Left, _Right));
	else
		return (exp(_Right * log(_CMPLX(_Ty)(_Left))));
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) pow(const _CMPLX(_Ty)& _Left, const _Ty& _Right)
	{	// return complex ^ real
	if (imag(_Left) == 0)
		return (_Pow(real(_Left), _Right));
	else
		return (exp(_Right * log(_Left)));
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) pow(const _Ty& _Left, const _CMPLX(_Ty)& _Right)
	{	// return real ^ complex
	if (imag(_Right) == 0)
		return (_Pow(_Left, real(_Right)));
	else if (0 < _Left)
		return (exp(_Right * _CTR(_Ty)::log(_Left)));
	else
		return (exp(_Right * log(_CMPLX(_Ty)(_Left))));
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) pow(const _CMPLX(_Ty)& _Left, const _CMPLX(_Ty)& _Right)
	{	// return complex ^ complex
	if (imag(_Left) == 0)
		return (pow(real(_Left), _Right));
	else
		return (exp(_Right * log(_Left)));
	}

		// TEMPLATE FUNCTION sinh
_TMPLT(_Ty) inline
	_CMPLX(_Ty) sinh(const _CMPLX(_Ty)& _Left)
	{	// return sinh(complex)
	return (_CMPLX(_Ty)(
		_CTR(_Ty)::_Sinh(real(_Left), _CTR(_Ty)::cos(imag(_Left))),
		_CTR(_Ty)::_Cosh(real(_Left), _CTR(_Ty)::sin(imag(_Left)))));
	}

		// TEMPLATE FUNCTION sqrt
_TMPLT(_Ty) inline
	_CMPLX(_Ty) sqrt(const _CMPLX(_Ty)& _Left)
	{	// return sqrt(complex)
	int _Leftexp;
	_Ty _Rho = _Fabs(_Left, &_Leftexp);	// get magnitude and scale factor

	if (_Leftexp == 0)
		return (_CMPLX(_Ty)(_Rho, _Rho));	// argument is zero, INF, or NaN
	else
		{	// compute in safest quadrant
		_Ty _Realmag = _CTR(_Ty)::ldexp(real(_Left) < 0
			? - real(_Left) : real(_Left), -_Leftexp);
		_Rho = _CTR(_Ty)::ldexp(_CTR(_Ty)::sqrt(
			2 * (_Realmag + _Rho)), _Leftexp / 2 - 1);

		if (0 <= real(_Left))
			return (_CMPLX(_Ty)(_Rho, imag(_Left) / (2 * _Rho)));
		else if (imag(_Left) < 0)
			return (_CMPLX(_Ty)(-imag(_Left) / (2 * _Rho), -_Rho));
		else
			return (_CMPLX(_Ty)(imag(_Left) / (2 * _Rho), _Rho));
		}
	}

		// TEMPLATE FUNCTION tanh
_TMPLT(_Ty) inline
	_CMPLX(_Ty) tanh(const _CMPLX(_Ty)& _Left)
	{	// return tanh(complex)
	_Ty _Tv = _CTR(_Ty)::tan(imag(_Left));
	_Ty _Sv = _CTR(_Ty)::_Sinh(real(_Left), (_Ty)(1));
	_Ty _Bv = _Sv *((_Ty)(1) + _Tv * _Tv);
	_Ty _Dv = (_Ty)(1) + _Bv * _Sv;

	if (_CTR(_Ty)::_Isinf(_Dv))
		return (_CMPLX(_Ty)(_Sv < (_Ty)0 ? (_Ty)(-1) : (_Ty)(1),
			_Tv * (_Ty)(0)));
	else
		return (_CMPLX(_Ty)((_CTR(_Ty)::sqrt((_Ty)(1) + _Sv * _Sv))
			* _Bv / _Dv, _Tv / _Dv));
	}

		// TEMPLATE FUNCTION arg
_TMPLT(_Ty) inline
	_Ty arg(const _CMPLX(_Ty)& _Left)
	{	// return phase angle of complex as real
	return (_CTR(_Ty)::atan2(imag(_Left), real(_Left)));
	}

		// TEMPLATE FUNCTION conj
_TMPLT(_Ty) inline
	_CMPLX(_Ty) conj(const _CMPLX(_Ty)& _Left)
	{	// return complex conjugate
	return (_CMPLX(_Ty)(real(_Left), -imag(_Left)));
	}

		// TEMPLATE FUNCTION proj
_TMPLT(_Ty) inline
	_CMPLX(_Ty) proj(const _CMPLX(_Ty)& _Left)
	{	// return complex projection
	return (_CMPLX(_Ty)(
		_CTR(_Ty)::_Isinf(real(_Left)) || _CTR(_Ty)::_Isinf(real(_Left))
			? _CTR(_Ty)::_Infv(real(_Left)) : real(_Left),
		imag(_Left) < 0 ? -(_Ty)0 : (_Ty)0));
	}

		// TEMPLATE FUNCTION cos
_TMPLT(_Ty) inline
	_CMPLX(_Ty) cos(const _CMPLX(_Ty)& _Left)
	{	// return cos(complex)
	return (_CMPLX(_Ty)(
		_CTR(_Ty)::_Cosh(imag(_Left), _CTR(_Ty)::cos(real(_Left))),
		-_CTR(_Ty)::_Sinh(imag(_Left),
			_CTR(_Ty)::sin(real(_Left)))));
	}

		// TEMPLATE FUNCTION log10
_TMPLT(_Ty) inline
	_CMPLX(_Ty) log10(const _CMPLX(_Ty)& _Left)
	{	// return log10(complex)
	return (log(_Left) * (_Ty)0.43429448190325182765112891891660508L);
	}

		// TEMPLATE FUNCTION norm
_TMPLT(_Ty) inline
	_Ty norm(const _CMPLX(_Ty)& _Left)
	{	// return squared magnitude
	return (real(_Left) * real(_Left) + imag(_Left) * imag(_Left));
	}

		// TEMPLATE FUNCTION polar
_TMPLT(_Ty) inline
	_CMPLX(_Ty) polar(const _Ty& _Rho, const _Ty& _Theta)
	{	// return _Rho * exp(i * _Theta) as complex
	return (_CMPLX(_Ty)(_Rho * _CTR(_Ty)::cos(_Theta),
		_Rho * _CTR(_Ty)::sin(_Theta)));
	}

_TMPLT(_Ty) inline
	_CMPLX(_Ty) polar(const _Ty& _Rho)
	{	// return _Rho * exp(i * 0) as complex
	return (_CMPLX(_Ty)(_Rho, (_Ty)0));
	}

		// TEMPLATE FUNCTION sin
_TMPLT(_Ty) inline
	_CMPLX(_Ty) sin(const _CMPLX(_Ty)& _Left)
	{	// return sin(complex)
	return (_CMPLX(_Ty)(
		_CTR(_Ty)::_Cosh(imag(_Left), _CTR(_Ty)::sin(real(_Left))),
		_CTR(_Ty)::_Sinh(imag(_Left), _CTR(_Ty)::cos(real(_Left)))));
	}

		// TEMPLATE FUNCTION tan
_TMPLT(_Ty) inline
	_CMPLX(_Ty) tan(const _CMPLX(_Ty)& _Left)
	{	// return tan(complex)
	_CMPLX(_Ty) _Zv(tanh(_CMPLX(_Ty)(-imag(_Left), real(_Left))));
	return (_CMPLX(_Ty)(imag(_Zv), -real(_Zv)));
	}

		// ADDITIONAL OVERLOADS
#define _GENERIC_MATHC0X(FUN, VAL) \
template<class _Ty> inline \
	typename enable_if<_Is_numeric<_Ty>::value, \
		typename _Promote_to_float<_Ty>::type>::type \
	FUN(_Ty) \
	{ \
	typedef typename _Promote_to_float<_Ty>::type type; \
	return ((type)VAL); \
	}

#define _GENERIC_MATHC1X(FUN, VAL) \
template<class _Ty> inline \
	typename enable_if<_Is_numeric<_Ty>::value, \
		typename _Promote_to_float<_Ty>::type>::type \
	FUN(_Ty _Left) \
	{ \
	typedef typename _Promote_to_float<_Ty>::type type; \
	return ((type)VAL); \
	}

_GENERIC_MATHC0X(arg, 0)
_GENERIC_MATHC0X(imag, 0)
_GENERIC_MATHC1X(real, _Left)

_GENERIC_MATHC1X(norm, (_Left * _Left))

_GENERIC_MATHC1X(conj, _Left)
_GENERIC_MATHC1X(proj, _Left)

		// TEMPLATE FUNCTION polar
template<class _Ty1,
	class _Ty2> inline
	complex<typename _Common_float_type<_Ty1, _Ty2>::type>
	polar(const _Ty1& _Left, const _Ty2& _Right)
	{	// bring mixed types to a common type
	typedef complex<typename _Common_float_type<_Ty1, _Ty2>::type> type;
	return (_STD polar((type)_Left, (type)_Right));
	}

		// TEMPLATE FUNCTION pow
template<class _Ty1,
	class _Ty2> inline
	complex<typename _Common_float_type<_Ty1, _Ty2>::type>
	pow(const complex<_Ty1>& _Left, const complex<_Ty2>& _Right)
	{	// bring mixed types to a common type
	typedef complex<typename _Common_float_type<_Ty1, _Ty2>::type> type;
	return (_STD pow(type(_Left), type(_Right)));
	}

template<class _Ty1,
	class _Ty2> inline
	complex<typename _Common_float_type<_Ty1, _Ty2>::type>
	pow(const complex<_Ty1>& _Left, const _Ty2& _Right)
	{	// bring mixed types to a common type
	typedef complex<typename _Common_float_type<_Ty1, _Ty2>::type> type;
	return (_STD pow(type(_Left), type(_Right)));
	}

template<class _Ty1,
	class _Ty2> inline
	complex<typename _Common_float_type<_Ty1, _Ty2>::type>
	pow(const _Ty1& _Left, const complex<_Ty2>& _Right)
	{	// bring mixed types to a common type
	typedef complex<typename _Common_float_type<_Ty1, _Ty2>::type> type;
	return (_STD pow(type(_Left), type(_Right)));
	}

template<class _Ty1,
	class _Ty2> inline
	typename enable_if<is_integral<_Ty1>::value && is_integral<_Ty2>::value,
		complex<_Ty1> >::type
	pow(const complex<_Ty1>& _Left, _Ty2& _Right)
	{	// raise Gaussian integer to an integer power
	typedef complex<_Ty1> type;
	type _Ans = type(1, 0);

	if (_Right < 0)
		_Ans = type(0, 0);	// ignore 1/type(0, 0) error
	else if (0 < _Right)
		{	// raise to a positive power
		for (type _Factor = _Left; ; _Factor *= _Factor)
			{	// fold in _Left^(2^N))
			if ((_Right & 1) != 0)
				_Ans *= _Factor;
			if ((_Right >>= 1) == 0)
				break;
			}
		}
	return (_Ans);
	}

/*
 * Copyright (c) 1992-2012 by P.J. Plauger.  ALL RIGHTS RESERVED.
 * Consult your license regarding permissions and restrictions.
V6.00:0009 */
