diff -r 000000000000 -r e4d67989cc36 genericopenlibs/cppstdlib/stl/src/complex.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/genericopenlibs/cppstdlib/stl/src/complex.cpp Tue Feb 02 02:01:42 2010 +0200 @@ -0,0 +1,347 @@ +/* + * Copyright (c) 1999 + * Silicon Graphics Computer Systems, Inc. + * + * Copyright (c) 1999 + * Boris Fomitchev + * + * This material is provided "as is", with absolutely no warranty expressed + * or implied. Any use is at your own risk. + * + * Permission to use or copy this software for any purpose is hereby granted + * without fee, provided the above notices are retained on all copies. + * Permission to modify the code and to distribute modified code is granted, + * provided the above notices are retained, and a notice that the code was + * modified is included with the above copyright notice. + * + */ + +#include "stlport_prefix.h" + +#include +#include +#include + +#if defined (_STLP_MSVC_LIB) && (_STLP_MSVC_LIB >= 1400) +// hypot is deprecated. +# if defined (_STLP_MSVC) +# pragma warning (disable : 4996) +# elif defined (__ICL) +# pragma warning (disable : 1478) +# endif +#endif + +_STLP_BEGIN_NAMESPACE + +// Complex division and square roots. + +// Absolute value +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC float _STLP_CALL abs(const complex& __z) +{ return ::hypot(__z._M_re, __z._M_im); } +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC double _STLP_CALL abs(const complex& __z) +{ return ::hypot(__z._M_re, __z._M_im); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC long double _STLP_CALL abs(const complex& __z) +{ return ::hypot(__z._M_re, __z._M_im); } +#endif + +// Phase + +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC float _STLP_CALL arg(const complex& __z) +{ return ::atan2(__z._M_im, __z._M_re); } + +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC double _STLP_CALL arg(const complex& __z) +{ return ::atan2(__z._M_im, __z._M_re); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC long double _STLP_CALL arg(const complex& __z) +{ return ::atan2(__z._M_im, __z._M_re); } +#endif + +// Construct a complex number from polar representation +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC complex _STLP_CALL polar(const float& __rho, const float& __phi) +{ return complex(__rho * ::cos(__phi), __rho * ::sin(__phi)); } +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC complex _STLP_CALL polar(const double& __rho, const double& __phi) +{ return complex(__rho * ::cos(__phi), __rho * ::sin(__phi)); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_TEMPLATE_NULL +_STLP_DECLSPEC complex _STLP_CALL polar(const long double& __rho, const long double& __phi) +{ return complex(__rho * ::cos(__phi), __rho * ::sin(__phi)); } +#endif + +// Division +template +static void _divT(const _Tp& __z1_r, const _Tp& __z1_i, + const _Tp& __z2_r, const _Tp& __z2_i, + _Tp& __res_r, _Tp& __res_i) { + _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r; + _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i; + + if (__ar <= __ai) { + _Tp __ratio = __z2_r / __z2_i; + _Tp __denom = __z2_i * (1 + __ratio * __ratio); + __res_r = (__z1_r * __ratio + __z1_i) / __denom; + __res_i = (__z1_i * __ratio - __z1_r) / __denom; + } + else { + _Tp __ratio = __z2_i / __z2_r; + _Tp __denom = __z2_r * (1 + __ratio * __ratio); + __res_r = (__z1_r + __z1_i * __ratio) / __denom; + __res_i = (__z1_i - __z1_r * __ratio) / __denom; + } +} + +template +static void _divT(const _Tp& __z1_r, + const _Tp& __z2_r, const _Tp& __z2_i, + _Tp& __res_r, _Tp& __res_i) { + _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r; + _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i; + + if (__ar <= __ai) { + _Tp __ratio = __z2_r / __z2_i; + _Tp __denom = __z2_i * (1 + __ratio * __ratio); + __res_r = (__z1_r * __ratio) / __denom; + __res_i = - __z1_r / __denom; + } + else { + _Tp __ratio = __z2_i / __z2_r; + _Tp __denom = __z2_r * (1 + __ratio * __ratio); + __res_r = __z1_r / __denom; + __res_i = - (__z1_r * __ratio) / __denom; + } +} + +_STLP_DECLSPEC void _STLP_CALL +complex::_div(const float& __z1_r, const float& __z1_i, + const float& __z2_r, const float& __z2_i, + float& __res_r, float& __res_i) +{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); } + +_STLP_DECLSPEC void _STLP_CALL +complex::_div(const float& __z1_r, + const float& __z2_r, const float& __z2_i, + float& __res_r, float& __res_i) +{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); } + + +_STLP_DECLSPEC void _STLP_CALL +complex::_div(const double& __z1_r, const double& __z1_i, + const double& __z2_r, const double& __z2_i, + double& __res_r, double& __res_i) +{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); } + +_STLP_DECLSPEC void _STLP_CALL +complex::_div(const double& __z1_r, + const double& __z2_r, const double& __z2_i, + double& __res_r, double& __res_i) +{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_DECLSPEC void _STLP_CALL +complex::_div(const long double& __z1_r, const long double& __z1_i, + const long double& __z2_r, const long double& __z2_i, + long double& __res_r, long double& __res_i) +{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); } + +_STLP_DECLSPEC void _STLP_CALL +complex::_div(const long double& __z1_r, + const long double& __z2_r, const long double& __z2_i, + long double& __res_r, long double& __res_i) +{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); } +#endif + +//---------------------------------------------------------------------- +// Square root +template +static complex<_Tp> sqrtT(const complex<_Tp>& z) { + _Tp re = z._M_re; + _Tp im = z._M_im; + _Tp mag = ::hypot(re, im); + complex<_Tp> result; + + if (mag == 0.f) { + result._M_re = result._M_im = 0.f; + } else if (re > 0.f) { + result._M_re = ::sqrt(0.5f * (mag + re)); + result._M_im = im/result._M_re/2.f; + } else { + result._M_im = ::sqrt(0.5f * (mag - re)); + if (im < 0.f) + result._M_im = - result._M_im; + result._M_re = im/result._M_im/2.f; + } + return result; +} + +_STLP_DECLSPEC complex _STLP_CALL +sqrt(const complex& z) { return sqrtT(z); } + +_STLP_DECLSPEC complex _STLP_CALL +sqrt(const complex& z) { return sqrtT(z); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_DECLSPEC complex _STLP_CALL +sqrt(const complex& z) { return sqrtT(z); } +#endif + +// exp, log, pow for complex, complex, and complex +//---------------------------------------------------------------------- +// exp +template +static complex<_Tp> expT(const complex<_Tp>& z) { + _Tp expx = ::exp(z._M_re); + return complex<_Tp>(expx * ::cos(z._M_im), + expx * ::sin(z._M_im)); +} +_STLP_DECLSPEC complex _STLP_CALL exp(const complex& z) +{ return expT(z); } + +_STLP_DECLSPEC complex _STLP_CALL exp(const complex& z) +{ return expT(z); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_DECLSPEC complex _STLP_CALL exp(const complex& z) +{ return expT(z); } +#endif + +//---------------------------------------------------------------------- +// log10 +template +static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) { + complex<_Tp> r; + + r._M_im = ::atan2(z._M_im, z._M_re) * ln10_inv; + r._M_re = ::log10(::hypot(z._M_re, z._M_im)); + return r; +} + +static const float LN10_INVF = 1.f / ::log(10.f); +_STLP_DECLSPEC complex _STLP_CALL log10(const complex& z) +{ return log10T(z, LN10_INVF); } + +static const double LN10_INV = 1. / ::log10(10.); +_STLP_DECLSPEC complex _STLP_CALL log10(const complex& z) +{ return log10T(z, LN10_INV); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +static const long double LN10_INVL = 1.l / ::log(10.l); +_STLP_DECLSPEC complex _STLP_CALL log10(const complex& z) +{ return log10T(z, LN10_INVL); } +#endif + +//---------------------------------------------------------------------- +// log +template +static complex<_Tp> logT(const complex<_Tp>& z) { + complex<_Tp> r; + + r._M_im = ::atan2(z._M_im, z._M_re); + r._M_re = ::log(::hypot(z._M_re, z._M_im)); + return r; +} +_STLP_DECLSPEC complex _STLP_CALL log(const complex& z) +{ return logT(z); } + +_STLP_DECLSPEC complex _STLP_CALL log(const complex& z) +{ return logT(z); } + +#ifndef _STLP_NO_LONG_DOUBLE +_STLP_DECLSPEC complex _STLP_CALL log(const complex& z) +{ return logT(z); } +# endif + +//---------------------------------------------------------------------- +// pow +template +static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) { + _Tp logr = ::log(a); + _Tp x = ::exp(logr * b._M_re); + _Tp y = logr * b._M_im; + + return complex<_Tp>(x * ::cos(y), x * ::sin(y)); +} + +template +static complex<_Tp> powT(const complex<_Tp>& z_in, int n) { + complex<_Tp> z = z_in; + z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >()); + if (n < 0) + return _Tp(1.0) / z; + else + return z; +} + +template +static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) { + _Tp logr = ::log(::hypot(a._M_re,a._M_im)); + _Tp logi = ::atan2(a._M_im, a._M_re); + _Tp x = ::exp(logr * b); + _Tp y = logi * b; + + return complex<_Tp>(x * ::cos(y), x * ::sin(y)); +} + +template +static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) { + _Tp logr = ::log(::hypot(a._M_re,a._M_im)); + _Tp logi = ::atan2(a._M_im, a._M_re); + _Tp x = ::exp(logr * b._M_re - logi * b._M_im); + _Tp y = logr * b._M_im + logi * b._M_re; + + return complex<_Tp>(x * ::cos(y), x * ::sin(y)); +} + +_STLP_DECLSPEC complex _STLP_CALL pow(const float& a, const complex& b) +{ return powT(a, b); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& z_in, int n) +{ return powT(z_in, n); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& a, const float& b) +{ return powT(a, b); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& a, const complex& b) +{ return powT(a, b); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const double& a, const complex& b) +{ return powT(a, b); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& z_in, int n) +{ return powT(z_in, n); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& a, const double& b) +{ return powT(a, b); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& a, const complex& b) +{ return powT(a, b); } + +#if !defined (_STLP_NO_LONG_DOUBLE) +_STLP_DECLSPEC complex _STLP_CALL pow(const long double& a, + const complex& b) +{ return powT(a, b); } + + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& z_in, int n) +{ return powT(z_in, n); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& a, + const long double& b) +{ return powT(a, b); } + +_STLP_DECLSPEC complex _STLP_CALL pow(const complex& a, + const complex& b) +{ return powT(a, b); } +#endif + +_STLP_END_NAMESPACE