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
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>
167static 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
187complex<float> _STLP_CALL
188sqrt(const complex<float>& z) { return sqrtT(z); }
189
190complex<double>  _STLP_CALL
191sqrt(const complex<double>& z) { return sqrtT(z); }
192
193#if !defined (_STLP_NO_LONG_DOUBLE)
194complex<long double> _STLP_CALL
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>
202static 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
220template <class _Tp>
221static 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
251template <class _Tp>
252static 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
272template <class _Tp>
273static 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
281template <class _Tp>
282static 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
291template <class _Tp>
292static 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
301template <class _Tp>
302static 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
354