ReactOS  0.4.14-dev-114-gc8cbd56
complex.cpp
Go to the documentation of this file.
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 
35 
36 // Complex division and square roots.
37 
38 // Absolute value
41 { return ::hypot(__z._M_re, __z._M_im); }
44 { return ::hypot(__z._M_re, __z._M_im); }
45 
46 #if !defined (_STLP_NO_LONG_DOUBLE)
49 { return ::hypot(__z._M_re, __z._M_im); }
50 #endif
51 
52 // Phase
53 
56 { return ::atan2(__z._M_im, __z._M_re); }
57 
60 { return ::atan2(__z._M_im, __z._M_re); }
61 
62 #if !defined (_STLP_NO_LONG_DOUBLE)
65 { return ::atan2(__z._M_im, __z._M_re); }
66 #endif
67 
68 // Construct a complex number from polar representation
70 _STLP_DECLSPEC complex<float> _STLP_CALL polar(const float& __rho, const float& __phi)
71 { return complex<float>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
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)
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>
168  _Tp re = z._M_re;
169  _Tp im = z._M_im;
170  _Tp mag = ::hypot(re, im);
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 
188 sqrt(const complex<float>& z) { return sqrtT(z); }
189 
191 sqrt(const complex<double>& z) { return sqrtT(z); }
192 
193 #if !defined (_STLP_NO_LONG_DOUBLE)
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 }
208 { return expT(z); }
209 
211 { return expT(z); }
212 
213 #if !defined (_STLP_NO_LONG_DOUBLE)
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);
231 { return log10T(z, LN10_INVF); }
232 
233 static const double LN10_INV = 1. / ::log10(10.);
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);
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 }
254 { return logT(z); }
255 
257 { return logT(z); }
258 
259 #ifndef _STLP_NO_LONG_DOUBLE
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 
306 { return powT(a, b); }
307 
309 { return powT(z_in, n); }
310 
312 { return powT(a, b); }
313 
315 { return powT(a, b); }
316 
318 { return powT(a, b); }
319 
321 { return powT(z_in, n); }
322 
324 { return powT(a, b); }
325 
327 { return powT(a, b); }
328 
329 #if !defined (_STLP_NO_LONG_DOUBLE)
331  const complex<long double>& b)
332 { return powT(a, b); }
333 
334 
336 { return powT(z_in, n); }
337 
339  const long double& b)
340 { return powT(a, b); }
341 
343  const complex<long double>& b)
344 { return powT(a, b); }
345 #endif
346 
_STLP_TEMPLATE_NULL _STLP_DECLSPEC complex< float > _STLP_CALL polar(const float &__rho, const float &__phi)
Definition: complex.cpp:70
static void _STLP_CALL _div(const value_type &__z1_r, const value_type &__z1_i, const value_type &__z2_r, const value_type &__z2_i, value_type &__res_r, value_type &__res_i)
Definition: _complex.c:44
valarray< _Tp > atan2(const valarray< _Tp > &__x, const valarray< _Tp > &__y)
Definition: _valarray.h:928
_STLP_TEMPLATE_NULL _STLP_DECLSPEC float _STLP_CALL arg(const complex< float > &__z)
Definition: complex.cpp:55
_Tp __power(_Tp __x, _Integer __n, _MonoidOperation __opr)
Definition: _numeric.c:74
static const float LN10_INVF
Definition: complex.cpp:229
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
static complex< _Tp > expT(const complex< _Tp > &z)
Definition: complex.cpp:202
GLdouble n
Definition: glext.h:7729
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
static complex< _Tp > sqrtT(const complex< _Tp > &z)
Definition: complex.cpp:167
static complex< _Tp > powT(const _Tp &a, const complex< _Tp > &b)
Definition: complex.cpp:267
_STLP_DECLSPEC complex< float > _STLP_CALL pow(const float &a, const complex< float > &b)
Definition: complex.cpp:305
complex< float > _STLP_CALL sqrt(const complex< float > &z)
Definition: complex.cpp:188
GLdouble GLdouble z
Definition: glext.h:5874
_STLP_DECLSPEC complex< float > _STLP_CALL cos(const complex< float > &)
_STLP_DECLSPEC complex< float > _STLP_CALL log10(const complex< float > &z)
Definition: complex.cpp:230
_Check_return_ __CRT_INLINE double hypot(_In_ double x, _In_ double y)
Definition: math.h:223
value_type _M_im
Definition: _complex.h:314
#define b
Definition: ke_i.h:79
r l[0]
Definition: byte_order.h:167
GLfloat f
Definition: glext.h:7540
GLboolean GLboolean GLboolean b
Definition: glext.h:6204
static complex< _Tp > log10T(const complex< _Tp > &z, const _Tp &ln10_inv)
Definition: complex.cpp:221
value_type _M_re
Definition: _complex.h:313
_STLP_BEGIN_NAMESPACE _STLP_TEMPLATE_NULL _STLP_DECLSPEC float _STLP_CALL abs(const complex< float > &__z)
Definition: complex.cpp:40
static void _divT(const _Tp &__z1_r, const _Tp &__z1_i, const _Tp &__z2_r, const _Tp &__z2_i, _Tp &__res_r, _Tp &__res_i)
Definition: complex.cpp:84
#define _STLP_PRIV
Definition: _dm.h:70
#define _STLP_DECLSPEC
Definition: features.h:982
_STLP_DECLSPEC complex< float > _STLP_CALL exp(const complex< float > &z)
Definition: complex.cpp:207
static complex< _Tp > logT(const complex< _Tp > &z)
Definition: complex.cpp:246
#define _STLP_TEMPLATE_NULL
Definition: features.h:652
value_type _M_re
Definition: _complex.h:451
#define _STLP_END_NAMESPACE
Definition: features.h:503
value_type _M_im
Definition: _complex.h:452
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
_STLP_DECLSPEC complex< float > _STLP_CALL log(const complex< float > &z)
Definition: complex.cpp:253
static const double LN10_INV
Definition: complex.cpp:233
static const long double LN10_INVL
Definition: complex.cpp:238
#define _STLP_BEGIN_NAMESPACE
Definition: features.h:501
_STLP_DECLSPEC complex< float > _STLP_CALL sin(const complex< float > &)
GLboolean GLboolean GLboolean GLboolean a
Definition: glext.h:6204
GLuint64EXT * result
Definition: glext.h:11304
#define _STLP_CALL
Definition: _bc.h:131