ReactOS 0.4.15-dev-8339-g4028de8
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
83template <class _Tp>
84static 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
104template <class _Tp>
105static 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
125void _STLP_CALL
126complex<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
131void _STLP_CALL
132complex<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
138void _STLP_CALL
139complex<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
144void _STLP_CALL
145complex<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)
151void _STLP_CALL
152complex<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
157void _STLP_CALL
158complex<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
166template <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
188sqrt(const complex<float>& z) { return sqrtT(z); }
189
191sqrt(const complex<double>& z) { return sqrtT(z); }
192
193#if !defined (_STLP_NO_LONG_DOUBLE)
195sqrt(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
201template <class _Tp>
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
220template <class _Tp>
221static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) {
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
229static const float LN10_INVF = 1.f / ::log(10.f);
231{ return log10T(z, LN10_INVF); }
232
233static const double LN10_INV = 1. / ::log10(10.);
235{ return log10T(z, LN10_INV); }
236
237#if !defined (_STLP_NO_LONG_DOUBLE)
238static const long double LN10_INVL = 1.l / ::log(10.l);
240{ return log10T(z, LN10_INVL); }
241#endif
242
243//----------------------------------------------------------------------
244// log
245template <class _Tp>
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
266template <class _Tp>
267static 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
275template <class _Tp>
276static 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
285template <class _Tp>
286static 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
295template <class _Tp>
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
#define _STLP_CALL
Definition: _bc.h:131
_STLP_DECLSPEC complex< float > _STLP_CALL cos(const complex< float > &)
_STLP_DECLSPEC complex< float > _STLP_CALL sin(const complex< float > &)
#define _STLP_PRIV
Definition: _dm.h:70
_Tp __power(_Tp __x, _Integer __n, _MonoidOperation __opr)
Definition: _numeric.c:74
valarray< _Tp > atan2(const valarray< _Tp > &__x, const valarray< _Tp > &__y)
Definition: _valarray.h:928
static const long double LN10_INVL
Definition: complex.cpp:238
static complex< _Tp > sqrtT(const complex< _Tp > &z)
Definition: complex.cpp:167
static const double LN10_INV
Definition: complex.cpp:233
complex< float > _STLP_CALL sqrt(const complex< float > &z)
Definition: complex.cpp:188
_STLP_DECLSPEC complex< float > _STLP_CALL log10(const complex< float > &z)
Definition: complex.cpp:230
static complex< _Tp > logT(const complex< _Tp > &z)
Definition: complex.cpp:246
_STLP_TEMPLATE_NULL _STLP_DECLSPEC complex< float > _STLP_CALL polar(const float &__rho, const float &__phi)
Definition: complex.cpp:70
static complex< _Tp > log10T(const complex< _Tp > &z, const _Tp &ln10_inv)
Definition: complex.cpp:221
static complex< _Tp > powT(const _Tp &a, const complex< _Tp > &b)
Definition: complex.cpp:267
static complex< _Tp > expT(const complex< _Tp > &z)
Definition: complex.cpp:202
static const float LN10_INVF
Definition: complex.cpp:229
_STLP_DECLSPEC complex< float > _STLP_CALL pow(const float &a, const complex< float > &b)
Definition: complex.cpp:305
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 abs(i)
Definition: fconv.c:206
#define _STLP_TEMPLATE_NULL
Definition: features.h:652
#define _STLP_BEGIN_NAMESPACE
Definition: features.h:501
#define _STLP_END_NAMESPACE
Definition: features.h:503
#define _STLP_DECLSPEC
Definition: features.h:982
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
GLdouble n
Definition: glext.h:7729
GLboolean GLboolean GLboolean b
Definition: glext.h:6204
GLboolean GLboolean GLboolean GLboolean a
Definition: glext.h:6204
GLuint64EXT * result
Definition: glext.h:11304
GLdouble GLdouble z
Definition: glext.h:5874
_Check_return_ __CRT_INLINE double hypot(_In_ double x, _In_ double y)
Definition: math.h:240
#define b
Definition: ke_i.h:79
DWORD exp
Definition: msg.c:16058
#define log(outFile, fmt,...)
Definition: util.h:15
value_type _M_im
Definition: _complex.h:171
value_type _M_re
Definition: _complex.h:170
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
void * arg
Definition: msvc.h:10