ReactOS Development > Doxygen

00001 /*
00002  * Copyright (c) 1999
00003  * Silicon Graphics Computer Systems, Inc.
00004  *
00005  * Copyright (c) 1999
00006  * Boris Fomitchev
00007  *
00008  * This material is provided "as is", with absolutely no warranty expressed
00009  * or implied. Any use is at your own risk.
00010  *
00011  * Permission to use or copy this software for any purpose is hereby granted
00012  * without fee, provided the above notices are retained on all copies.
00013  * Permission to modify the code and to distribute modified code is granted,
00014  * provided the above notices are retained, and a notice that the code was
00015  * modified is included with the above copyright notice.
00016  *
00017  */
00018 
00019 #include "stlport_prefix.h"
00020 
00021 #include <numeric>
00022 #include <cmath>
00023 #include <complex>
00024 
00025 #if defined (_STLP_MSVC_LIB) && (_STLP_MSVC_LIB >= 1400)
00026 // hypot is deprecated.
00027 #  if defined (_STLP_MSVC)
00028 #    pragma warning (disable : 4996)
00029 #  elif defined (__ICL)
00030 #    pragma warning (disable : 1478)
00031 #  endif
00032 #endif
00033 
00034 _STLP_BEGIN_NAMESPACE
00035 
00036 // Complex division and square roots.
00037 
00038 // Absolute value
00039 _STLP_TEMPLATE_NULL
00040 _STLP_DECLSPEC float _STLP_CALL abs(const complex<float>& __z)
00041 { return ::hypot(__z._M_re, __z._M_im); }
00042 _STLP_TEMPLATE_NULL
00043 _STLP_DECLSPEC double _STLP_CALL abs(const complex<double>& __z)
00044 { return ::hypot(__z._M_re, __z._M_im); }
00045 
00046 #if !defined (_STLP_NO_LONG_DOUBLE)
00047 _STLP_TEMPLATE_NULL
00048 _STLP_DECLSPEC long double _STLP_CALL abs(const complex<long double>& __z)
00049 { return ::hypot(__z._M_re, __z._M_im); }
00050 #endif
00051 
00052 // Phase
00053 
00054 _STLP_TEMPLATE_NULL
00055 _STLP_DECLSPEC float _STLP_CALL arg(const complex<float>& __z)
00056 { return ::atan2(__z._M_im, __z._M_re); }
00057 
00058 _STLP_TEMPLATE_NULL
00059 _STLP_DECLSPEC double _STLP_CALL arg(const complex<double>& __z)
00060 { return ::atan2(__z._M_im, __z._M_re); }
00061 
00062 #if !defined (_STLP_NO_LONG_DOUBLE)
00063 _STLP_TEMPLATE_NULL
00064 _STLP_DECLSPEC long double _STLP_CALL arg(const complex<long double>& __z)
00065 { return ::atan2(__z._M_im, __z._M_re); }
00066 #endif
00067 
00068 // Construct a complex number from polar representation
00069 _STLP_TEMPLATE_NULL
00070 _STLP_DECLSPEC complex<float> _STLP_CALL polar(const float& __rho, const float& __phi)
00071 { return complex<float>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
00072 _STLP_TEMPLATE_NULL
00073 _STLP_DECLSPEC complex<double> _STLP_CALL polar(const double& __rho, const double& __phi)
00074 { return complex<double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
00075 
00076 #if !defined (_STLP_NO_LONG_DOUBLE)
00077 _STLP_TEMPLATE_NULL
00078 _STLP_DECLSPEC complex<long double> _STLP_CALL polar(const long double& __rho, const long double& __phi)
00079 { return complex<long double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
00080 #endif
00081 
00082 // Division
00083 template <class _Tp>
00084 static void _divT(const _Tp& __z1_r, const _Tp& __z1_i,
00085                   const _Tp& __z2_r, const _Tp& __z2_i,
00086                   _Tp& __res_r, _Tp& __res_i) {
00087   _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
00088   _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
00089 
00090   if (__ar <= __ai) {
00091     _Tp __ratio = __z2_r / __z2_i;
00092     _Tp __denom = __z2_i * (1 + __ratio * __ratio);
00093     __res_r = (__z1_r * __ratio + __z1_i) / __denom;
00094     __res_i = (__z1_i * __ratio - __z1_r) / __denom;
00095   }
00096   else {
00097     _Tp __ratio = __z2_i / __z2_r;
00098     _Tp __denom = __z2_r * (1 + __ratio * __ratio);
00099     __res_r = (__z1_r + __z1_i * __ratio) / __denom;
00100     __res_i = (__z1_i - __z1_r * __ratio) / __denom;
00101   }
00102 }
00103 
00104 template <class _Tp>
00105 static void _divT(const _Tp& __z1_r,
00106                   const _Tp& __z2_r, const _Tp& __z2_i,
00107                   _Tp& __res_r, _Tp& __res_i) {
00108   _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
00109   _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
00110 
00111   if (__ar <= __ai) {
00112     _Tp __ratio = __z2_r / __z2_i;
00113     _Tp __denom = __z2_i * (1 + __ratio * __ratio);
00114     __res_r = (__z1_r * __ratio) / __denom;
00115     __res_i = - __z1_r / __denom;
00116   }
00117   else {
00118     _Tp __ratio = __z2_i / __z2_r;
00119     _Tp __denom = __z2_r * (1 + __ratio * __ratio);
00120     __res_r = __z1_r / __denom;
00121     __res_i = - (__z1_r * __ratio) / __denom;
00122   }
00123 }
00124 
00125 void _STLP_CALL
00126 complex<float>::_div(const float& __z1_r, const float& __z1_i,
00127                      const float& __z2_r, const float& __z2_i,
00128                      float& __res_r, float& __res_i)
00129 { _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
00130 
00131 void _STLP_CALL
00132 complex<float>::_div(const float& __z1_r,
00133                      const float& __z2_r, const float& __z2_i,
00134                      float& __res_r, float& __res_i)
00135 { _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
00136 
00137 
00138 void  _STLP_CALL
00139 complex<double>::_div(const double& __z1_r, const double& __z1_i,
00140                       const double& __z2_r, const double& __z2_i,
00141                       double& __res_r, double& __res_i)
00142 { _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
00143 
00144 void _STLP_CALL
00145 complex<double>::_div(const double& __z1_r,
00146                       const double& __z2_r, const double& __z2_i,
00147                       double& __res_r, double& __res_i)
00148 { _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
00149 
00150 #if !defined (_STLP_NO_LONG_DOUBLE)
00151 void  _STLP_CALL
00152 complex<long double>::_div(const long double& __z1_r, const long double& __z1_i,
00153                            const long double& __z2_r, const long double& __z2_i,
00154                            long double& __res_r, long double& __res_i)
00155 { _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
00156 
00157 void _STLP_CALL
00158 complex<long double>::_div(const long double& __z1_r,
00159                            const long double& __z2_r, const long double& __z2_i,
00160                            long double& __res_r, long double& __res_i)
00161 { _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
00162 #endif
00163 
00164 //----------------------------------------------------------------------
00165 // Square root
00166 template <class _Tp>
00167 static complex<_Tp> sqrtT(const complex<_Tp>& z) {
00168   _Tp re = z._M_re;
00169   _Tp im = z._M_im;
00170   _Tp mag = ::hypot(re, im);
00171   complex<_Tp> result;
00172 
00173   if (mag == 0.f) {
00174     result._M_re = result._M_im = 0.f;
00175   } else if (re > 0.f) {
00176     result._M_re = ::sqrt(0.5f * (mag + re));
00177     result._M_im = im/result._M_re/2.f;
00178   } else {
00179     result._M_im = ::sqrt(0.5f * (mag - re));
00180     if (im < 0.f)
00181       result._M_im = - result._M_im;
00182     result._M_re = im/result._M_im/2.f;
00183   }
00184   return result;
00185 }
00186 
00187 complex<float> _STLP_CALL
00188 sqrt(const complex<float>& z) { return sqrtT(z); }
00189 
00190 complex<double>  _STLP_CALL
00191 sqrt(const complex<double>& z) { return sqrtT(z); }
00192 
00193 #if !defined (_STLP_NO_LONG_DOUBLE)
00194 complex<long double> _STLP_CALL
00195 sqrt(const complex<long double>& z) { return sqrtT(z); }
00196 #endif
00197 
00198 // exp, log, pow for complex<float>, complex<double>, and complex<long double>
00199 //----------------------------------------------------------------------
00200 // exp
00201 template <class _Tp>
00202 static complex<_Tp> expT(const complex<_Tp>& z) {
00203   _Tp expx = ::exp(z._M_re);
00204   return complex<_Tp>(expx * ::cos(z._M_im),
00205                       expx * ::sin(z._M_im));
00206 }
00207 _STLP_DECLSPEC complex<float>  _STLP_CALL exp(const complex<float>& z)
00208 { return expT(z); }
00209 
00210 _STLP_DECLSPEC complex<double> _STLP_CALL exp(const complex<double>& z)
00211 { return expT(z); }
00212 
00213 #if !defined (_STLP_NO_LONG_DOUBLE)
00214 _STLP_DECLSPEC complex<long double> _STLP_CALL exp(const complex<long double>& z)
00215 { return expT(z); }
00216 #endif
00217 
00218 //----------------------------------------------------------------------
00219 // log10
00220 template <class _Tp>
00221 static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) {
00222   complex<_Tp> r;
00223 
00224   r._M_im = ::atan2(z._M_im, z._M_re) * ln10_inv;
00225   r._M_re = ::log10(::hypot(z._M_re, z._M_im));
00226   return r;
00227 }
00228 
00229 static const float LN10_INVF = 1.f / ::log(10.f);
00230 _STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z)
00231 { return log10T(z, LN10_INVF); }
00232 
00233 static const double LN10_INV = 1. / ::log10(10.);
00234 _STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z)
00235 { return log10T(z, LN10_INV); }
00236 
00237 #if !defined (_STLP_NO_LONG_DOUBLE)
00238 static const long double LN10_INVL = 1.l / ::log(10.l);
00239 _STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z)
00240 { return log10T(z, LN10_INVL); }
00241 #endif
00242 
00243 //----------------------------------------------------------------------
00244 // log
00245 template <class _Tp>
00246 static complex<_Tp> logT(const complex<_Tp>& z) {
00247   complex<_Tp> r;
00248 
00249   r._M_im = ::atan2(z._M_im, z._M_re);
00250   r._M_re = ::log(::hypot(z._M_re, z._M_im));
00251   return r;
00252 }
00253 _STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z)
00254 { return logT(z); }
00255 
00256 _STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)
00257 { return logT(z); }
00258 
00259 #ifndef _STLP_NO_LONG_DOUBLE
00260 _STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z)
00261 { return logT(z); }
00262 # endif
00263 
00264 //----------------------------------------------------------------------
00265 // pow
00266 template <class _Tp>
00267 static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) {
00268   _Tp logr = ::log(a);
00269   _Tp x = ::exp(logr * b._M_re);
00270   _Tp y = logr * b._M_im;
00271 
00272   return complex<_Tp>(x * ::cos(y), x * ::sin(y));
00273 }
00274 
00275 template <class _Tp>
00276 static complex<_Tp> powT(const complex<_Tp>& z_in, int n) {
00277   complex<_Tp> z = z_in;
00278   z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >());
00279   if (n < 0)
00280     return _Tp(1.0) / z;
00281   else
00282     return z;
00283 }
00284 
00285 template <class _Tp>
00286 static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) {
00287   _Tp logr = ::log(::hypot(a._M_re,a._M_im));
00288   _Tp logi = ::atan2(a._M_im, a._M_re);
00289   _Tp x = ::exp(logr * b);
00290   _Tp y = logi * b;
00291 
00292   return complex<_Tp>(x * ::cos(y), x * ::sin(y));
00293 }
00294 
00295 template <class _Tp>
00296 static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) {
00297   _Tp logr = ::log(::hypot(a._M_re,a._M_im));
00298   _Tp logi = ::atan2(a._M_im, a._M_re);
00299   _Tp x = ::exp(logr * b._M_re - logi * b._M_im);
00300   _Tp y = logr * b._M_im + logi * b._M_re;
00301 
00302   return complex<_Tp>(x * ::cos(y), x * ::sin(y));
00303 }
00304 
00305 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)
00306 { return powT(a, b); }
00307 
00308 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)
00309 { return powT(z_in, n); }
00310 
00311 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)
00312 { return powT(a, b); }
00313 
00314 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)
00315 { return powT(a, b); }
00316 
00317 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)
00318 { return powT(a, b); }
00319 
00320 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)
00321 { return powT(z_in, n); }
00322 
00323 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)
00324 { return powT(a, b); }
00325 
00326 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b)
00327 { return powT(a, b); }
00328 
00329 #if !defined (_STLP_NO_LONG_DOUBLE)
00330 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a,
00331                                                    const complex<long double>& b)
00332 { return powT(a, b); }
00333 
00334 
00335 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)
00336 { return powT(z_in, n); }
00337 
00338 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
00339                                                    const long double& b)
00340 { return powT(a, b); }
00341 
00342 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
00343                                                    const complex<long double>& b)
00344 { return powT(a, b); }
00345 #endif
00346 
00347 _STLP_END_NAMESPACE
Generated on Fri May 25 2012 04:33:39 for ReactOS by
1.7.6.1