1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/build/stlport/src/complex.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,353 @@
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 +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 +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 +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 +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 +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 +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 +complex<float> _STLP_CALL
1.191 +sqrt(const complex<float>& z) { return sqrtT(z); }
1.192 +
1.193 +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 +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 +_STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z)
1.233 +{
1.234 + const float LN10_INVF = 1.f / ::log(10.f);
1.235 + return log10T(z, LN10_INVF);
1.236 +}
1.237 +
1.238 +_STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z)
1.239 +{
1.240 + const double LN10_INV = 1. / ::log10(10.);
1.241 + return log10T(z, LN10_INV);
1.242 +}
1.243 +
1.244 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.245 +_STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z)
1.246 +{
1.247 + const long double LN10_INVL = 1.l / ::log(10.l);
1.248 + return log10T(z, LN10_INVL);
1.249 +}
1.250 +#endif
1.251 +
1.252 +//----------------------------------------------------------------------
1.253 +// log
1.254 +template <class _Tp>
1.255 +static complex<_Tp> logT(const complex<_Tp>& z) {
1.256 + complex<_Tp> r;
1.257 +
1.258 + r._M_im = ::atan2(z._M_im, z._M_re);
1.259 + r._M_re = ::log(::hypot(z._M_re, z._M_im));
1.260 + return r;
1.261 +}
1.262 +_STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z)
1.263 +{ return logT(z); }
1.264 +
1.265 +_STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)
1.266 +{ return logT(z); }
1.267 +
1.268 +#ifndef _STLP_NO_LONG_DOUBLE
1.269 +_STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z)
1.270 +{ return logT(z); }
1.271 +# endif
1.272 +
1.273 +//----------------------------------------------------------------------
1.274 +// pow
1.275 +template <class _Tp>
1.276 +static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) {
1.277 + _Tp logr = ::log(a);
1.278 + _Tp x = ::exp(logr * b._M_re);
1.279 + _Tp y = logr * b._M_im;
1.280 +
1.281 + return complex<_Tp>(x * ::cos(y), x * ::sin(y));
1.282 +}
1.283 +
1.284 +template <class _Tp>
1.285 +static complex<_Tp> powT(const complex<_Tp>& z_in, int n) {
1.286 + complex<_Tp> z = z_in;
1.287 + z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >());
1.288 + if (n < 0)
1.289 + return _Tp(1.0) / z;
1.290 + else
1.291 + return z;
1.292 +}
1.293 +
1.294 +template <class _Tp>
1.295 +static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) {
1.296 + _Tp logr = ::log(::hypot(a._M_re,a._M_im));
1.297 + _Tp logi = ::atan2(a._M_im, a._M_re);
1.298 + _Tp x = ::exp(logr * b);
1.299 + _Tp y = logi * b;
1.300 +
1.301 + return complex<_Tp>(x * ::cos(y), x * ::sin(y));
1.302 +}
1.303 +
1.304 +template <class _Tp>
1.305 +static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) {
1.306 + _Tp logr = ::log(::hypot(a._M_re,a._M_im));
1.307 + _Tp logi = ::atan2(a._M_im, a._M_re);
1.308 + _Tp x = ::exp(logr * b._M_re - logi * b._M_im);
1.309 + _Tp y = logr * b._M_im + logi * b._M_re;
1.310 +
1.311 + return complex<_Tp>(x * ::cos(y), x * ::sin(y));
1.312 +}
1.313 +
1.314 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)
1.315 +{ return powT(a, b); }
1.316 +
1.317 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)
1.318 +{ return powT(z_in, n); }
1.319 +
1.320 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)
1.321 +{ return powT(a, b); }
1.322 +
1.323 +_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)
1.324 +{ return powT(a, b); }
1.325 +
1.326 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)
1.327 +{ return powT(a, b); }
1.328 +
1.329 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)
1.330 +{ return powT(z_in, n); }
1.331 +
1.332 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)
1.333 +{ return powT(a, b); }
1.334 +
1.335 +_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b)
1.336 +{ return powT(a, b); }
1.337 +
1.338 +#if !defined (_STLP_NO_LONG_DOUBLE)
1.339 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a,
1.340 + const complex<long double>& b)
1.341 +{ return powT(a, b); }
1.342 +
1.343 +
1.344 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)
1.345 +{ return powT(z_in, n); }
1.346 +
1.347 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
1.348 + const long double& b)
1.349 +{ return powT(a, b); }
1.350 +
1.351 +_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
1.352 + const complex<long double>& b)
1.353 +{ return powT(a, b); }
1.354 +#endif
1.355 +
1.356 +_STLP_END_NAMESPACE
