1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/genericopenlibs/cppstdlib/stl/src/complex.cpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,347 @@
1.4 +/*
1.5 + * Copyright (c) 1999
1.6 + * Silicon Graphics Computer Systems, Inc.
1.7 + *
1.8 + * Copyright (c) 1999
1.9 + * Boris Fomitchev
1.10 + *
1.11 + * This material is provided "as is", with absolutely no warranty expressed
1.12 + * or implied. Any use is at your own risk.
1.13 + *
1.14 + * Permission to use or copy this software for any purpose is hereby granted
1.15 + * without fee, provided the above notices are retained on all copies.
1.16 + * Permission to modify the code and to distribute modified code is granted,
1.17 + * provided the above notices are retained, and a notice that the code was
1.18 + * modified is included with the above copyright notice.
1.19 + *
1.20 + */
1.21 +
1.22 +#include "stlport_prefix.h"
1.23 +
1.24 +#include <numeric>
1.25 +#include <cmath>
1.26 +#include <complex>
1.27 +
1.28 +#if defined (_STLP_MSVC_LIB) && (_STLP_MSVC_LIB >= 1400)
1.29 +// hypot is deprecated.
1.30 +# if defined (_STLP_MSVC)
1.31 +# pragma warning (disable : 4996)
1.32 +# elif defined (__ICL)
1.33 +# pragma warning (disable : 1478)
1.34 +# endif
1.35 +#endif
1.36 +
1.37 +_STLP_BEGIN_NAMESPACE
1.38 +
1.39 +// Complex division and square roots.
1.40 +
1.41 +// Absolute value
1.42 +_STLP_TEMPLATE_NULL
1.43 +_STLP_DECLSPEC float _STLP_CALL abs(const complex<float>& __z)
1.44 +{ return ::hypot(__z._M_re, __z._M_im); }
1.45 +_STLP_TEMPLATE_NULL
1.46 +_STLP_DECLSPEC double _STLP_CALL abs(const complex<double>& __z)
1.47 +{ return ::hypot(__z._M_re, __z._M_im); }
1.48 +
1.49 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.50 +_STLP_TEMPLATE_NULL
1.51 +_STLP_DECLSPEC long double _STLP_CALL abs(const complex<long double>& __z)
1.52 +{ return ::hypot(__z._M_re, __z._M_im); }
1.53 +#endif
1.54 +
1.55 +// Phase
1.56 +
1.57 +_STLP_TEMPLATE_NULL
1.58 +_STLP_DECLSPEC float _STLP_CALL arg(const complex<float>& __z)
1.59 +{ return ::atan2(__z._M_im, __z._M_re); }
1.60 +
1.61 +_STLP_TEMPLATE_NULL
1.62 +_STLP_DECLSPEC double _STLP_CALL arg(const complex<double>& __z)
1.63 +{ return ::atan2(__z._M_im, __z._M_re); }
1.64 +
1.65 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.66 +_STLP_TEMPLATE_NULL
1.67 +_STLP_DECLSPEC long double _STLP_CALL arg(const complex<long double>& __z)
1.68 +{ return ::atan2(__z._M_im, __z._M_re); }
1.69 +#endif
1.70 +
1.71 +// Construct a complex number from polar representation
1.72 +_STLP_TEMPLATE_NULL
1.73 +_STLP_DECLSPEC complex<float> _STLP_CALL polar(const float& __rho, const float& __phi)
1.74 +{ return complex<float>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
1.75 +_STLP_TEMPLATE_NULL
1.76 +_STLP_DECLSPEC complex<double> _STLP_CALL polar(const double& __rho, const double& __phi)
1.77 +{ return complex<double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
1.78 +
1.79 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.80 +_STLP_TEMPLATE_NULL
1.81 +_STLP_DECLSPEC complex<long double> _STLP_CALL polar(const long double& __rho, const long double& __phi)
1.82 +{ return complex<long double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
1.83 +#endif
1.84 +
1.85 +// Division
1.86 +template <class _Tp>
1.87 +static void _divT(const _Tp& __z1_r, const _Tp& __z1_i,
1.88 + const _Tp& __z2_r, const _Tp& __z2_i,
1.89 + _Tp& __res_r, _Tp& __res_i) {
1.90 + _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
1.91 + _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
1.92 +
1.93 + if (__ar <= __ai) {
1.94 + _Tp __ratio = __z2_r / __z2_i;
1.95 + _Tp __denom = __z2_i * (1 + __ratio * __ratio);
1.96 + __res_r = (__z1_r * __ratio + __z1_i) / __denom;
1.97 + __res_i = (__z1_i * __ratio - __z1_r) / __denom;
1.98 + }
1.99 + else {
1.100 + _Tp __ratio = __z2_i / __z2_r;
1.101 + _Tp __denom = __z2_r * (1 + __ratio * __ratio);
1.102 + __res_r = (__z1_r + __z1_i * __ratio) / __denom;
1.103 + __res_i = (__z1_i - __z1_r * __ratio) / __denom;
1.104 + }
1.105 +}
1.106 +
1.107 +template <class _Tp>
1.108 +static void _divT(const _Tp& __z1_r,
1.109 + const _Tp& __z2_r, const _Tp& __z2_i,
1.110 + _Tp& __res_r, _Tp& __res_i) {
1.111 + _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
1.112 + _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
1.113 +
1.114 + if (__ar <= __ai) {
1.115 + _Tp __ratio = __z2_r / __z2_i;
1.116 + _Tp __denom = __z2_i * (1 + __ratio * __ratio);
1.117 + __res_r = (__z1_r * __ratio) / __denom;
1.118 + __res_i = - __z1_r / __denom;
1.119 + }
1.120 + else {
1.121 + _Tp __ratio = __z2_i / __z2_r;
1.122 + _Tp __denom = __z2_r * (1 + __ratio * __ratio);
1.123 + __res_r = __z1_r / __denom;
1.124 + __res_i = - (__z1_r * __ratio) / __denom;
1.125 + }
1.126 +}
1.127 +
1.128 +_STLP_DECLSPEC void _STLP_CALL
1.129 +complex<float>::_div(const float& __z1_r, const float& __z1_i,
1.130 + const float& __z2_r, const float& __z2_i,
1.131 + float& __res_r, float& __res_i)
1.132 +{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
1.133 +
1.134 +_STLP_DECLSPEC void _STLP_CALL
1.135 +complex<float>::_div(const float& __z1_r,
1.136 + const float& __z2_r, const float& __z2_i,
1.137 + float& __res_r, float& __res_i)
1.138 +{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
1.139 +
1.140 +
1.141 +_STLP_DECLSPEC void _STLP_CALL
1.142 +complex<double>::_div(const double& __z1_r, const double& __z1_i,
1.143 + const double& __z2_r, const double& __z2_i,
1.144 + double& __res_r, double& __res_i)
1.145 +{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
1.146 +
1.147 +_STLP_DECLSPEC void _STLP_CALL
1.148 +complex<double>::_div(const double& __z1_r,
1.149 + const double& __z2_r, const double& __z2_i,
1.150 + double& __res_r, double& __res_i)
1.151 +{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
1.152 +
1.153 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.154 +_STLP_DECLSPEC void _STLP_CALL
1.155 +complex<long double>::_div(const long double& __z1_r, const long double& __z1_i,
1.156 + const long double& __z2_r, const long double& __z2_i,
1.157 + long double& __res_r, long double& __res_i)
1.158 +{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
1.159 +
1.160 +_STLP_DECLSPEC void _STLP_CALL
1.161 +complex<long double>::_div(const long double& __z1_r,
1.162 + const long double& __z2_r, const long double& __z2_i,
1.163 + long double& __res_r, long double& __res_i)
1.164 +{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
1.165 +#endif
1.166 +
1.167 +//----------------------------------------------------------------------
1.168 +// Square root
1.169 +template <class _Tp>
1.170 +static complex<_Tp> sqrtT(const complex<_Tp>& z) {
1.171 + _Tp re = z._M_re;
1.172 + _Tp im = z._M_im;
1.173 + _Tp mag = ::hypot(re, im);
1.174 + complex<_Tp> result;
1.175 +
1.176 + if (mag == 0.f) {
1.177 + result._M_re = result._M_im = 0.f;
1.178 + } else if (re > 0.f) {
1.179 + result._M_re = ::sqrt(0.5f * (mag + re));
1.180 + result._M_im = im/result._M_re/2.f;
1.181 + } else {
1.182 + result._M_im = ::sqrt(0.5f * (mag - re));
1.183 + if (im < 0.f)
1.184 + result._M_im = - result._M_im;
1.185 + result._M_re = im/result._M_im/2.f;
1.186 + }
1.187 + return result;
1.188 +}
1.189 +
1.190 +_STLP_DECLSPEC complex<float> _STLP_CALL
1.191 +sqrt(const complex<float>& z) { return sqrtT(z); }
1.192 +
1.193 +_STLP_DECLSPEC complex<double> _STLP_CALL
1.194 +sqrt(const complex<double>& z) { return sqrtT(z); }
1.195 +
1.196 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.197 +_STLP_DECLSPEC complex<long double> _STLP_CALL
1.198 +sqrt(const complex<long double>& z) { return sqrtT(z); }
1.199 +#endif
1.200 +
1.201 +// exp, log, pow for complex<float>, complex<double>, and complex<long double>
1.202 +//----------------------------------------------------------------------
1.203 +// exp
1.204 +template <class _Tp>
1.205 +static complex<_Tp> expT(const complex<_Tp>& z) {
1.206 + _Tp expx = ::exp(z._M_re);
1.207 + return complex<_Tp>(expx * ::cos(z._M_im),
1.208 + expx * ::sin(z._M_im));
1.209 +}
1.210 +_STLP_DECLSPEC complex<float> _STLP_CALL exp(const complex<float>& z)
1.211 +{ return expT(z); }
1.212 +
1.213 +_STLP_DECLSPEC complex<double> _STLP_CALL exp(const complex<double>& z)
1.214 +{ return expT(z); }
1.215 +
1.216 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.217 +_STLP_DECLSPEC complex<long double> _STLP_CALL exp(const complex<long double>& z)
1.218 +{ return expT(z); }
1.219 +#endif
1.220 +
1.221 +//----------------------------------------------------------------------
1.222 +// log10
1.223 +template <class _Tp>
1.224 +static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) {
1.225 + complex<_Tp> r;
1.226 +
1.227 + r._M_im = ::atan2(z._M_im, z._M_re) * ln10_inv;
1.228 + r._M_re = ::log10(::hypot(z._M_re, z._M_im));
1.229 + return r;
1.230 +}
1.231 +
1.232 +static const float LN10_INVF = 1.f / ::log(10.f);
1.233 +_STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z)
1.234 +{ return log10T(z, LN10_INVF); }
1.235 +
1.236 +static const double LN10_INV = 1. / ::log10(10.);
1.237 +_STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z)
1.238 +{ return log10T(z, LN10_INV); }
1.239 +
1.240 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.241 +static const long double LN10_INVL = 1.l / ::log(10.l);
1.242 +_STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z)
1.243 +{ return log10T(z, LN10_INVL); }
1.244 +#endif
1.245 +
1.246 +//----------------------------------------------------------------------
1.247 +// log
1.248 +template <class _Tp>
1.249 +static complex<_Tp> logT(const complex<_Tp>& z) {
1.250 + complex<_Tp> r;
1.251 +
1.252 + r._M_im = ::atan2(z._M_im, z._M_re);
1.253 + r._M_re = ::log(::hypot(z._M_re, z._M_im));
1.254 + return r;
1.255 +}
1.256 +_STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z)
1.257 +{ return logT(z); }
1.258 +
1.259 +_STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)
1.260 +{ return logT(z); }
1.261 +
1.262 +#ifndef _STLP_NO_LONG_DOUBLE
1.263 +_STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z)
1.264 +{ return logT(z); }
1.265 +# endif
1.266 +
1.267 +//----------------------------------------------------------------------
1.268 +// pow
1.269 +template <class _Tp>
1.270 +static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) {
1.271 + _Tp logr = ::log(a);
1.272 + _Tp x = ::exp(logr * b._M_re);
1.273 + _Tp y = logr * b._M_im;
1.274 +
1.275 + return complex<_Tp>(x * ::cos(y), x * ::sin(y));
1.276 +}
1.277 +
1.278 +template <class _Tp>
1.279 +static complex<_Tp> powT(const complex<_Tp>& z_in, int n) {
1.280 + complex<_Tp> z = z_in;
1.281 + z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >());
1.282 + if (n < 0)
1.283 + return _Tp(1.0) / z;
1.284 + else
1.285 + return z;
1.286 +}
1.287 +
1.288 +template <class _Tp>
1.289 +static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) {
1.290 + _Tp logr = ::log(::hypot(a._M_re,a._M_im));
1.291 + _Tp logi = ::atan2(a._M_im, a._M_re);
1.292 + _Tp x = ::exp(logr * b);
1.293 + _Tp y = logi * b;
1.294 +
1.295 + return complex<_Tp>(x * ::cos(y), x * ::sin(y));
1.296 +}
1.297 +
1.298 +template <class _Tp>
1.299 +static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) {
1.300 + _Tp logr = ::log(::hypot(a._M_re,a._M_im));
1.301 + _Tp logi = ::atan2(a._M_im, a._M_re);
1.302 + _Tp x = ::exp(logr * b._M_re - logi * b._M_im);
1.303 + _Tp y = logr * b._M_im + logi * b._M_re;
1.304 +
1.305 + return complex<_Tp>(x * ::cos(y), x * ::sin(y));
1.306 +}
1.307 +
1.308 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)
1.309 +{ return powT(a, b); }
1.310 +
1.311 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)
1.312 +{ return powT(z_in, n); }
1.313 +
1.314 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)
1.315 +{ return powT(a, b); }
1.316 +
1.317 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)
1.318 +{ return powT(a, b); }
1.319 +
1.320 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)
1.321 +{ return powT(a, b); }
1.322 +
1.323 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)
1.324 +{ return powT(z_in, n); }
1.325 +
1.326 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)
1.327 +{ return powT(a, b); }
1.328 +
1.329 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b)
1.330 +{ return powT(a, b); }
1.331 +
1.332 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.333 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a,
1.334 + const complex<long double>& b)
1.335 +{ return powT(a, b); }
1.336 +
1.337 +
1.338 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)
1.339 +{ return powT(z_in, n); }
1.340 +
1.341 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
1.342 + const long double& b)
1.343 +{ return powT(a, b); }
1.344 +
1.345 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
1.346 + const complex<long double>& b)
1.347 +{ return powT(a, b); }
1.348 +#endif
1.349 +
1.350 +_STLP_END_NAMESPACE