|
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 _STLP_DECLSPEC 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 _STLP_DECLSPEC 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 _STLP_DECLSPEC 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 _STLP_DECLSPEC 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 _STLP_DECLSPEC 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 _STLP_DECLSPEC 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 _STLP_DECLSPEC complex<float> _STLP_CALL |
|
188 sqrt(const complex<float>& z) { return sqrtT(z); } |
|
189 |
|
190 _STLP_DECLSPEC complex<double> _STLP_CALL |
|
191 sqrt(const complex<double>& z) { return sqrtT(z); } |
|
192 |
|
193 #if !defined (_STLP_NO_LONG_DOUBLE) |
|
194 _STLP_DECLSPEC 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 static const float LN10_INVF = 1.f / ::log(10.f); |
|
230 _STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z) |
|
231 { return log10T(z, LN10_INVF); } |
|
232 |
|
233 static const double LN10_INV = 1. / ::log10(10.); |
|
234 _STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z) |
|
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); |
|
239 _STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z) |
|
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 } |
|
253 _STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z) |
|
254 { return logT(z); } |
|
255 |
|
256 _STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z) |
|
257 { return logT(z); } |
|
258 |
|
259 #ifndef _STLP_NO_LONG_DOUBLE |
|
260 _STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z) |
|
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 |
|
305 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b) |
|
306 { return powT(a, b); } |
|
307 |
|
308 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n) |
|
309 { return powT(z_in, n); } |
|
310 |
|
311 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b) |
|
312 { return powT(a, b); } |
|
313 |
|
314 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b) |
|
315 { return powT(a, b); } |
|
316 |
|
317 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b) |
|
318 { return powT(a, b); } |
|
319 |
|
320 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n) |
|
321 { return powT(z_in, n); } |
|
322 |
|
323 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b) |
|
324 { return powT(a, b); } |
|
325 |
|
326 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b) |
|
327 { return powT(a, b); } |
|
328 |
|
329 #if !defined (_STLP_NO_LONG_DOUBLE) |
|
330 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a, |
|
331 const complex<long double>& b) |
|
332 { return powT(a, b); } |
|
333 |
|
334 |
|
335 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n) |
|
336 { return powT(z_in, n); } |
|
337 |
|
338 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a, |
|
339 const long double& b) |
|
340 { return powT(a, b); } |
|
341 |
|
342 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a, |
|
343 const complex<long double>& b) |
|
344 { return powT(a, b); } |
|
345 #endif |
|
346 |
|
347 _STLP_END_NAMESPACE |