lib/c++/stlport/complex.cpp

   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  */
  18
  19 #include "stlport_prefix.h"
  20
  21 #include <numeric>
  22 #include <cmath>
  23 #include <complex>
  24
  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
  33
  34 _STLP_BEGIN_NAMESPACE
  35
  36 // Complex division and square roots.
  37
  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); }
  45
  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
  51
  52 // Phase
  53
  54 _STLP_TEMPLATE_NULL
  55 _STLP_DECLSPEC float _STLP_CALL arg(const complex<float>& __z)
  56 { return ::atan2(__z._M_im, __z._M_re); }
  57
  58 _STLP_TEMPLATE_NULL
  59 _STLP_DECLSPEC double _STLP_CALL arg(const complex<double>& __z)
  60 { return ::atan2(__z._M_im, __z._M_re); }
  61
  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
  67
  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)); }
  75
  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
  81
  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;
  89
  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 }
 103
 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;
 110
 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 }
 124
 125 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); }
 130
 131 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); }
 136
 137
 138 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); }
 143
 144 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); }
 149
 150 #if !defined (_STLP_NO_LONG_DOUBLE)
 151 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); }
 156
 157 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
 163
 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;
 172
 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 }
 186
 187 complex<float> _STLP_CALL
 188 sqrt(const complex<float>& z) { return sqrtT(z); }
 189
 190 complex<double>  _STLP_CALL
 191 sqrt(const complex<double>& z) { return sqrtT(z); }
 192
 193 #if !defined (_STLP_NO_LONG_DOUBLE)
 194 complex<long double> _STLP_CALL
 195 sqrt(const complex<long double>& z) { return sqrtT(z); }
 196 #endif
 197
 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); }
 209
 210 _STLP_DECLSPEC complex<double> _STLP_CALL exp(const complex<double>& z)
 211 { return expT(z); }
 212
 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
 217
 218 //----------------------------------------------------------------------
 219 // log10
 220 template <class _Tp>
 221 static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) {
 222   complex<_Tp> r;
 223
 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 }
 228
 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); }
 232
 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); }
 236
 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
 242
 243 //----------------------------------------------------------------------
 244 // log
 245 template <class _Tp>
 246 static complex<_Tp> logT(const complex<_Tp>& z) {
 247   complex<_Tp> r;
 248
 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); }
 255
 256 _STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)
 257 { return logT(z); }
 258
 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
 263
 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;
 271
 272   return complex<_Tp>(x * ::cos(y), x * ::sin(y));
 273 }
 274
 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 }
 284
 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;
 291
 292   return complex<_Tp>(x * ::cos(y), x * ::sin(y));
 293 }
 294
 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;
 301
 302   return complex<_Tp>(x * ::cos(y), x * ::sin(y));
 303 }
 304
 305 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)
 306 { return powT(a, b); }
 307
 308 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)
 309 { return powT(z_in, n); }
 310
 311 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)
 312 { return powT(a, b); }
 313
 314 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)
 315 { return powT(a, b); }
 316
 317 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)
 318 { return powT(a, b); }
 319
 320 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)
 321 { return powT(z_in, n); }
 322
 323 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)
 324 { return powT(a, b); }
 325
 326 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b)
 327 { return powT(a, b); }
 328
 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); }
 333
 334
 335 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)
 336 { return powT(z_in, n); }
 337
 338 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
 339                                                    const long double& b)
 340 { return powT(a, b); }
 341
 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
 346
 347 _STLP_END_NAMESPACE