The Tor Browser / file comparison
comparison: build/stlport/src/complex.cpp
build/stlport/src/complex.cpp
- changeset 0
- 6474c204b198
equal
deleted
inserted
replaced
| -1:000000000000 | 0:9fcfcaec45bd |
|---|---|
| 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 |
