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
abs(const complex<float> & __z)40 _STLP_DECLSPEC float _STLP_CALL abs(const complex<float>& __z)
41 { return ::hypot(__z._M_re, __z._M_im); }
42 _STLP_TEMPLATE_NULL
abs(const complex<double> & __z)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
abs(const complex<long double> & __z)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
arg(const complex<float> & __z)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
arg(const complex<double> & __z)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
arg(const complex<long double> & __z)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
polar(const float & __rho,const float & __phi)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
polar(const double & __rho,const double & __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)
77 _STLP_TEMPLATE_NULL
polar(const long double & __rho,const long double & __phi)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>
_divT(const _Tp & __z1_r,const _Tp & __z1_i,const _Tp & __z2_r,const _Tp & __z2_i,_Tp & __res_r,_Tp & __res_i)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>
_divT(const _Tp & __z1_r,const _Tp & __z2_r,const _Tp & __z2_i,_Tp & __res_r,_Tp & __res_i)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
_div(const float & __z1_r,const float & __z1_i,const float & __z2_r,const float & __z2_i,float & __res_r,float & __res_i)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
_div(const float & __z1_r,const float & __z2_r,const float & __z2_i,float & __res_r,float & __res_i)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
_div(const double & __z1_r,const double & __z1_i,const double & __z2_r,const double & __z2_i,double & __res_r,double & __res_i)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
_div(const double & __z1_r,const double & __z2_r,const double & __z2_i,double & __res_r,double & __res_i)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
_div(const long double & __z1_r,const long double & __z1_i,const long double & __z2_r,const long double & __z2_i,long double & __res_r,long double & __res_i)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
_div(const long double & __z1_r,const long double & __z2_r,const long double & __z2_i,long double & __res_r,long double & __res_i)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>
sqrtT(const complex<_Tp> & z)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
sqrt(const complex<float> & z)188 sqrt(const complex<float>& z) { return sqrtT(z); }
189
190 complex<double> _STLP_CALL
sqrt(const complex<double> & z)191 sqrt(const complex<double>& z) { return sqrtT(z); }
192
193 #if !defined (_STLP_NO_LONG_DOUBLE)
194 complex<long double> _STLP_CALL
sqrt(const complex<long double> & z)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>
expT(const complex<_Tp> & z)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 }
exp(const complex<float> & z)207 _STLP_DECLSPEC complex<float> _STLP_CALL exp(const complex<float>& z)
208 { return expT(z); }
209
exp(const complex<double> & z)210 _STLP_DECLSPEC complex<double> _STLP_CALL exp(const complex<double>& z)
211 { return expT(z); }
212
213 #if !defined (_STLP_NO_LONG_DOUBLE)
exp(const complex<long double> & z)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>
log10T(const complex<_Tp> & z,const _Tp & ln10_inv)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
log10(const complex<float> & z)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
log10(const complex<double> & z)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)
log10(const complex<long double> & z)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>
logT(const complex<_Tp> & z)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 }
log(const complex<float> & z)259 _STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z)
260 { return logT(z); }
261
log(const complex<double> & z)262 _STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)
263 { return logT(z); }
264
265 #ifndef _STLP_NO_LONG_DOUBLE
log(const complex<long double> & z)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>
powT(const _Tp & a,const complex<_Tp> & b)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>
powT(const complex<_Tp> & z_in,int n)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>
powT(const complex<_Tp> & a,const _Tp & b)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>
powT(const complex<_Tp> & a,const complex<_Tp> & b)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
pow(const float & a,const complex<float> & b)311 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)
312 { return powT(a, b); }
313
pow(const complex<float> & z_in,int n)314 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)
315 { return powT(z_in, n); }
316
pow(const complex<float> & a,const float & b)317 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)
318 { return powT(a, b); }
319
pow(const complex<float> & a,const complex<float> & b)320 _STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)
321 { return powT(a, b); }
322
pow(const double & a,const complex<double> & b)323 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)
324 { return powT(a, b); }
325
pow(const complex<double> & z_in,int n)326 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)
327 { return powT(z_in, n); }
328
pow(const complex<double> & a,const double & b)329 _STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)
330 { return powT(a, b); }
331
pow(const complex<double> & a,const complex<double> & b)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)
pow(const long double & a,const complex<long double> & b)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
pow(const complex<long double> & z_in,int n)341 _STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)
342 { return powT(z_in, n); }
343
pow(const complex<long double> & a,const long double & b)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
pow(const complex<long double> & a,const complex<long double> & b)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