/home/ntakagi/work/STLport-5.1.5/src/complex.cppGo to the documentation of this file.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 Mon Mar 10 15:32:16 2008 by ![]() |