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 _STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z) 230 { 231 const float LN10_INVF = 1.f / ::log(10.f); 232 return log10T(z, LN10_INVF); 233 } 234 235 _STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z) 236 { 237 const double LN10_INV = 1. / ::log10(10.); 238 return log10T(z, LN10_INV); 239 } 240 241 #if !defined (_STLP_NO_LONG_DOUBLE) 242 _STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z) 243 { 244 const long double LN10_INVL = 1.l / ::log(10.l); 245 return log10T(z, LN10_INVL); 246 } 247 #endif 248 249 //---------------------------------------------------------------------- 250 // log 251 template <class _Tp> 252 static complex<_Tp> logT(const complex<_Tp>& z) { 253 complex<_Tp> r; 254 255 r._M_im = ::atan2(z._M_im, z._M_re); 256 r._M_re = ::log(::hypot(z._M_re, z._M_im)); 257 return r; 258 } 259 _STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z) 260 { return logT(z); } 261 262 _STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z) 263 { return logT(z); } 264 265 #ifndef _STLP_NO_LONG_DOUBLE 266 _STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z) 267 { return logT(z); } 268 # endif 269 270 //---------------------------------------------------------------------- 271 // pow 272 template <class _Tp> 273 static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) { 274 _Tp logr = ::log(a); 275 _Tp x = ::exp(logr * b._M_re); 276 _Tp y = logr * b._M_im; 277 278 return complex<_Tp>(x * ::cos(y), x * ::sin(y)); 279 } 280 281 template <class _Tp> 282 static complex<_Tp> powT(const complex<_Tp>& z_in, int n) { 283 complex<_Tp> z = z_in; 284 z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >()); 285 if (n < 0) 286 return _Tp(1.0) / z; 287 else 288 return z; 289 } 290 291 template <class _Tp> 292 static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) { 293 _Tp logr = ::log(::hypot(a._M_re,a._M_im)); 294 _Tp logi = ::atan2(a._M_im, a._M_re); 295 _Tp x = ::exp(logr * b); 296 _Tp y = logi * b; 297 298 return complex<_Tp>(x * ::cos(y), x * ::sin(y)); 299 } 300 301 template <class _Tp> 302 static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) { 303 _Tp logr = ::log(::hypot(a._M_re,a._M_im)); 304 _Tp logi = ::atan2(a._M_im, a._M_re); 305 _Tp x = ::exp(logr * b._M_re - logi * b._M_im); 306 _Tp y = logr * b._M_im + logi * b._M_re; 307 308 return complex<_Tp>(x * ::cos(y), x * ::sin(y)); 309 } 310 311 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b) 312 { return powT(a, b); } 313 314 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n) 315 { return powT(z_in, n); } 316 317 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b) 318 { return powT(a, b); } 319 320 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b) 321 { return powT(a, b); } 322 323 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b) 324 { return powT(a, b); } 325 326 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n) 327 { return powT(z_in, n); } 328 329 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b) 330 { return powT(a, b); } 331 332 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b) 333 { return powT(a, b); } 334 335 #if !defined (_STLP_NO_LONG_DOUBLE) 336 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a, 337 const complex<long double>& b) 338 { return powT(a, b); } 339 340 341 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n) 342 { return powT(z_in, n); } 343 344 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a, 345 const long double& b) 346 { return powT(a, b); } 347 348 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a, 349 const complex<long double>& b) 350 { return powT(a, b); } 351 #endif 352 353 _STLP_END_NAMESPACE 354
