Symaptic: os/ossrv/genericopenlibs/cppstdlib/stl/src/complex.cpp@260cb5ec6c19

     1 /*

     2  * Copyright (c) 1999

     3  * Silicon Graphics Computer Systems, Inc.

     4  *

     5  * Copyright (c) 1999

     6  * Boris Fomitchev

     7  *

     8  * This material is provided "as is", with absolutely no warranty expressed

     9  * or implied. Any use is at your own risk.

    10  *

    11  * Permission to use or copy this software for any purpose is hereby granted

    12  * without fee, provided the above notices are retained on all copies.

    13  * Permission to modify the code and to distribute modified code is granted,

    14  * provided the above notices are retained, and a notice that the code was

    15  * modified is included with the above copyright notice.

    16  *

    17  */

    19 #include "stlport_prefix.h"

    21 #include <numeric>

    22 #include <cmath>

    23 #include <complex>

    25 #if defined (_STLP_MSVC_LIB) && (_STLP_MSVC_LIB >= 1400)

    26 // hypot is deprecated.

    27 #  if defined (_STLP_MSVC)

    28 #    pragma warning (disable : 4996)

    29 #  elif defined (__ICL)

    30 #    pragma warning (disable : 1478)

    31 #  endif

    32 #endif

    34 _STLP_BEGIN_NAMESPACE

    36 // Complex division and square roots.

    38 // Absolute value

    39 _STLP_TEMPLATE_NULL

    40 _STLP_DECLSPEC float _STLP_CALL abs(const complex<float>& __z)

    41 { return ::hypot(__z._M_re, __z._M_im); }

    42 _STLP_TEMPLATE_NULL

    43 _STLP_DECLSPEC double _STLP_CALL abs(const complex<double>& __z)

    44 { return ::hypot(__z._M_re, __z._M_im); }

    46 #if !defined (_STLP_NO_LONG_DOUBLE)

    47 _STLP_TEMPLATE_NULL

    48 _STLP_DECLSPEC long double _STLP_CALL abs(const complex<long double>& __z)

    49 { return ::hypot(__z._M_re, __z._M_im); }

    50 #endif

    52 // Phase

    54 _STLP_TEMPLATE_NULL

    55 _STLP_DECLSPEC float _STLP_CALL arg(const complex<float>& __z)

    56 { return ::atan2(__z._M_im, __z._M_re); }

    58 _STLP_TEMPLATE_NULL

    59 _STLP_DECLSPEC double _STLP_CALL arg(const complex<double>& __z)

    60 { return ::atan2(__z._M_im, __z._M_re); }

    62 #if !defined (_STLP_NO_LONG_DOUBLE)

    63 _STLP_TEMPLATE_NULL

    64 _STLP_DECLSPEC long double _STLP_CALL arg(const complex<long double>& __z)

    65 { return ::atan2(__z._M_im, __z._M_re); }

    66 #endif

    68 // Construct a complex number from polar representation

    69 _STLP_TEMPLATE_NULL

    70 _STLP_DECLSPEC complex<float> _STLP_CALL polar(const float& __rho, const float& __phi)

    71 { return complex<float>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }

    72 _STLP_TEMPLATE_NULL

    73 _STLP_DECLSPEC complex<double> _STLP_CALL polar(const double& __rho, const double& __phi)

    74 { return complex<double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }

    76 #if !defined (_STLP_NO_LONG_DOUBLE)

    77 _STLP_TEMPLATE_NULL

    78 _STLP_DECLSPEC complex<long double> _STLP_CALL polar(const long double& __rho, const long double& __phi)

    79 { return complex<long double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }

    80 #endif

    82 // Division

    83 template <class _Tp>

    84 static void _divT(const _Tp& __z1_r, const _Tp& __z1_i,

    85                   const _Tp& __z2_r, const _Tp& __z2_i,

    86                   _Tp& __res_r, _Tp& __res_i) {

    87   _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;

    88   _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;

    90   if (__ar <= __ai) {

    91     _Tp __ratio = __z2_r / __z2_i;

    92     _Tp __denom = __z2_i * (1 + __ratio * __ratio);

    93     __res_r = (__z1_r * __ratio + __z1_i) / __denom;

    94     __res_i = (__z1_i * __ratio - __z1_r) / __denom;

    95   }

    96   else {

    97     _Tp __ratio = __z2_i / __z2_r;

    98     _Tp __denom = __z2_r * (1 + __ratio * __ratio);

    99     __res_r = (__z1_r + __z1_i * __ratio) / __denom;

   100     __res_i = (__z1_i - __z1_r * __ratio) / __denom;

   101   }

   102 }

   104 template <class _Tp>

   105 static void _divT(const _Tp& __z1_r,

   106                   const _Tp& __z2_r, const _Tp& __z2_i,

   107                   _Tp& __res_r, _Tp& __res_i) {

   108   _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;

   109   _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;

   111   if (__ar <= __ai) {

   112     _Tp __ratio = __z2_r / __z2_i;

   113     _Tp __denom = __z2_i * (1 + __ratio * __ratio);

   114     __res_r = (__z1_r * __ratio) / __denom;

   115     __res_i = - __z1_r / __denom;

   116   }

   117   else {

   118     _Tp __ratio = __z2_i / __z2_r;

   119     _Tp __denom = __z2_r * (1 + __ratio * __ratio);

   120     __res_r = __z1_r / __denom;

   121     __res_i = - (__z1_r * __ratio) / __denom;

   122   }

   123 }

   125 _STLP_DECLSPEC void _STLP_CALL

   126 complex<float>::_div(const float& __z1_r, const float& __z1_i,

   127                      const float& __z2_r, const float& __z2_i,

   128                      float& __res_r, float& __res_i)

   129 { _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }

   131 _STLP_DECLSPEC void _STLP_CALL

   132 complex<float>::_div(const float& __z1_r,

   133                      const float& __z2_r, const float& __z2_i,

   134                      float& __res_r, float& __res_i)

   135 { _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }

   138 _STLP_DECLSPEC void _STLP_CALL

   139 complex<double>::_div(const double& __z1_r, const double& __z1_i,

   140                       const double& __z2_r, const double& __z2_i,

   141                       double& __res_r, double& __res_i)

   142 { _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }

   144 _STLP_DECLSPEC void _STLP_CALL

   145 complex<double>::_div(const double& __z1_r,

   146                       const double& __z2_r, const double& __z2_i,

   147                       double& __res_r, double& __res_i)

   148 { _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }

   150 #if !defined (_STLP_NO_LONG_DOUBLE)

   151 _STLP_DECLSPEC void _STLP_CALL

   152 complex<long double>::_div(const long double& __z1_r, const long double& __z1_i,

   153                            const long double& __z2_r, const long double& __z2_i,

   154                            long double& __res_r, long double& __res_i)

   155 { _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }

   157 _STLP_DECLSPEC void _STLP_CALL

   158 complex<long double>::_div(const long double& __z1_r,

   159                            const long double& __z2_r, const long double& __z2_i,

   160                            long double& __res_r, long double& __res_i)

   161 { _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }

   162 #endif

   164 //----------------------------------------------------------------------

   165 // Square root

   166 template <class _Tp>

   167 static complex<_Tp> sqrtT(const complex<_Tp>& z) {

   168   _Tp re = z._M_re;

   169   _Tp im = z._M_im;

   170   _Tp mag = ::hypot(re, im);

   171   complex<_Tp> result;

   173   if (mag == 0.f) {

   174     result._M_re = result._M_im = 0.f;

   175   } else if (re > 0.f) {

   176     result._M_re = ::sqrt(0.5f * (mag + re));

   177     result._M_im = im/result._M_re/2.f;

   178   } else {

   179     result._M_im = ::sqrt(0.5f * (mag - re));

   180     if (im < 0.f)

   181       result._M_im = - result._M_im;

   182     result._M_re = im/result._M_im/2.f;

   183   }

   184   return result;

   185 }

   187 _STLP_DECLSPEC complex<float> _STLP_CALL

   188 sqrt(const complex<float>& z) { return sqrtT(z); }

   190 _STLP_DECLSPEC complex<double>  _STLP_CALL

   191 sqrt(const complex<double>& z) { return sqrtT(z); }

   193 #if !defined (_STLP_NO_LONG_DOUBLE)

   194 _STLP_DECLSPEC complex<long double> _STLP_CALL

   195 sqrt(const complex<long double>& z) { return sqrtT(z); }

   196 #endif

   198 // exp, log, pow for complex<float>, complex<double>, and complex<long double>

   199 //----------------------------------------------------------------------

   200 // exp

   201 template <class _Tp>

   202 static complex<_Tp> expT(const complex<_Tp>& z) {

   203   _Tp expx = ::exp(z._M_re);

   204   return complex<_Tp>(expx * ::cos(z._M_im),

   205                       expx * ::sin(z._M_im));

   206 }

   207 _STLP_DECLSPEC complex<float>  _STLP_CALL exp(const complex<float>& z)

   208 { return expT(z); }

   210 _STLP_DECLSPEC complex<double> _STLP_CALL exp(const complex<double>& z)

   211 { return expT(z); }

   213 #if !defined (_STLP_NO_LONG_DOUBLE)

   214 _STLP_DECLSPEC complex<long double> _STLP_CALL exp(const complex<long double>& z)

   215 { return expT(z); }

   216 #endif

   218 //----------------------------------------------------------------------

   219 // log10

   220 template <class _Tp>

   221 static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) {

   222   complex<_Tp> r;

   224   r._M_im = ::atan2(z._M_im, z._M_re) * ln10_inv;

   225   r._M_re = ::log10(::hypot(z._M_re, z._M_im));

   226   return r;

   227 }

   229 static const float LN10_INVF = 1.f / ::log(10.f);

   230 _STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z)

   231 { return log10T(z, LN10_INVF); }

   233 static const double LN10_INV = 1. / ::log10(10.);

   234 _STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z)

   235 { return log10T(z, LN10_INV); }

   237 #if !defined (_STLP_NO_LONG_DOUBLE)

   238 static const long double LN10_INVL = 1.l / ::log(10.l);

   239 _STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z)

   240 { return log10T(z, LN10_INVL); }

   241 #endif

   243 //----------------------------------------------------------------------

   244 // log

   245 template <class _Tp>

   246 static complex<_Tp> logT(const complex<_Tp>& z) {

   247   complex<_Tp> r;

   249   r._M_im = ::atan2(z._M_im, z._M_re);

   250   r._M_re = ::log(::hypot(z._M_re, z._M_im));

   251   return r;

   252 }

   253 _STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z)

   254 { return logT(z); }

   256 _STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)

   257 { return logT(z); }

   259 #ifndef _STLP_NO_LONG_DOUBLE

   260 _STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z)

   261 { return logT(z); }

   262 # endif

   264 //----------------------------------------------------------------------

   265 // pow

   266 template <class _Tp>

   267 static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) {

   268   _Tp logr = ::log(a);

   269   _Tp x = ::exp(logr * b._M_re);

   270   _Tp y = logr * b._M_im;

   272   return complex<_Tp>(x * ::cos(y), x * ::sin(y));

   273 }

   275 template <class _Tp>

   276 static complex<_Tp> powT(const complex<_Tp>& z_in, int n) {

   277   complex<_Tp> z = z_in;

   278   z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >());

   279   if (n < 0)

   280     return _Tp(1.0) / z;

   281   else

   282     return z;

   283 }

   285 template <class _Tp>

   286 static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) {

   287   _Tp logr = ::log(::hypot(a._M_re,a._M_im));

   288   _Tp logi = ::atan2(a._M_im, a._M_re);

   289   _Tp x = ::exp(logr * b);

   290   _Tp y = logi * b;

   292   return complex<_Tp>(x * ::cos(y), x * ::sin(y));

   293 }

   295 template <class _Tp>

   296 static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) {

   297   _Tp logr = ::log(::hypot(a._M_re,a._M_im));

   298   _Tp logi = ::atan2(a._M_im, a._M_re);

   299   _Tp x = ::exp(logr * b._M_re - logi * b._M_im);

   300   _Tp y = logr * b._M_im + logi * b._M_re;

   302   return complex<_Tp>(x * ::cos(y), x * ::sin(y));

   303 }

   305 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)

   306 { return powT(a, b); }

   308 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)

   309 { return powT(z_in, n); }

   311 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)

   312 { return powT(a, b); }

   314 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)

   315 { return powT(a, b); }

   317 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)

   318 { return powT(a, b); }

   320 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)

   321 { return powT(z_in, n); }

   323 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)

   324 { return powT(a, b); }

   326 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b)

   327 { return powT(a, b); }

   329 #if !defined (_STLP_NO_LONG_DOUBLE)

   330 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a,

   331                                                    const complex<long double>& b)

   332 { return powT(a, b); }

   335 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)

   336 { return powT(z_in, n); }

   338 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,

   339                                                    const long double& b)

   340 { return powT(a, b); }

   342 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,

   343                                                    const complex<long double>& b)

   344 { return powT(a, b); }

   345 #endif

   347 _STLP_END_NAMESPACE

author	sl
	Tue, 10 Jun 2014 14:32:02 +0200
changeset 1	260cb5ec6c19
permissions	-rw-r--r--