os/ossrv/genericopenlibs/cppstdlib/stl/src/complex.cpp
author sl@SLION-WIN7.fritz.box
Fri, 15 Jun 2012 03:10:57 +0200
changeset 0 bde4ae8d615e
permissions -rw-r--r--
First public contribution.
sl@0
     1
/*
sl@0
     2
 * Copyright (c) 1999
sl@0
     3
 * Silicon Graphics Computer Systems, Inc.
sl@0
     4
 *
sl@0
     5
 * Copyright (c) 1999
sl@0
     6
 * Boris Fomitchev
sl@0
     7
 *
sl@0
     8
 * This material is provided "as is", with absolutely no warranty expressed
sl@0
     9
 * or implied. Any use is at your own risk.
sl@0
    10
 *
sl@0
    11
 * Permission to use or copy this software for any purpose is hereby granted
sl@0
    12
 * without fee, provided the above notices are retained on all copies.
sl@0
    13
 * Permission to modify the code and to distribute modified code is granted,
sl@0
    14
 * provided the above notices are retained, and a notice that the code was
sl@0
    15
 * modified is included with the above copyright notice.
sl@0
    16
 *
sl@0
    17
 */
sl@0
    18
sl@0
    19
#include "stlport_prefix.h"
sl@0
    20
sl@0
    21
#include <numeric>
sl@0
    22
#include <cmath>
sl@0
    23
#include <complex>
sl@0
    24
sl@0
    25
#if defined (_STLP_MSVC_LIB) && (_STLP_MSVC_LIB >= 1400)
sl@0
    26
// hypot is deprecated.
sl@0
    27
#  if defined (_STLP_MSVC)
sl@0
    28
#    pragma warning (disable : 4996)
sl@0
    29
#  elif defined (__ICL)
sl@0
    30
#    pragma warning (disable : 1478)
sl@0
    31
#  endif
sl@0
    32
#endif
sl@0
    33
sl@0
    34
_STLP_BEGIN_NAMESPACE
sl@0
    35
sl@0
    36
// Complex division and square roots.
sl@0
    37
sl@0
    38
// Absolute value
sl@0
    39
_STLP_TEMPLATE_NULL
sl@0
    40
_STLP_DECLSPEC float _STLP_CALL abs(const complex<float>& __z)
sl@0
    41
{ return ::hypot(__z._M_re, __z._M_im); }
sl@0
    42
_STLP_TEMPLATE_NULL
sl@0
    43
_STLP_DECLSPEC double _STLP_CALL abs(const complex<double>& __z)
sl@0
    44
{ return ::hypot(__z._M_re, __z._M_im); }
sl@0
    45
sl@0
    46
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
    47
_STLP_TEMPLATE_NULL
sl@0
    48
_STLP_DECLSPEC long double _STLP_CALL abs(const complex<long double>& __z)
sl@0
    49
{ return ::hypot(__z._M_re, __z._M_im); }
sl@0
    50
#endif
sl@0
    51
sl@0
    52
// Phase
sl@0
    53
sl@0
    54
_STLP_TEMPLATE_NULL
sl@0
    55
_STLP_DECLSPEC float _STLP_CALL arg(const complex<float>& __z)
sl@0
    56
{ return ::atan2(__z._M_im, __z._M_re); }
sl@0
    57
sl@0
    58
_STLP_TEMPLATE_NULL
sl@0
    59
_STLP_DECLSPEC double _STLP_CALL arg(const complex<double>& __z)
sl@0
    60
{ return ::atan2(__z._M_im, __z._M_re); }
sl@0
    61
sl@0
    62
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
    63
_STLP_TEMPLATE_NULL
sl@0
    64
_STLP_DECLSPEC long double _STLP_CALL arg(const complex<long double>& __z)
sl@0
    65
{ return ::atan2(__z._M_im, __z._M_re); }
sl@0
    66
#endif
sl@0
    67
sl@0
    68
// Construct a complex number from polar representation
sl@0
    69
_STLP_TEMPLATE_NULL
sl@0
    70
_STLP_DECLSPEC complex<float> _STLP_CALL polar(const float& __rho, const float& __phi)
sl@0
    71
{ return complex<float>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
sl@0
    72
_STLP_TEMPLATE_NULL
sl@0
    73
_STLP_DECLSPEC complex<double> _STLP_CALL polar(const double& __rho, const double& __phi)
sl@0
    74
{ return complex<double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
sl@0
    75
sl@0
    76
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
    77
_STLP_TEMPLATE_NULL
sl@0
    78
_STLP_DECLSPEC complex<long double> _STLP_CALL polar(const long double& __rho, const long double& __phi)
sl@0
    79
{ return complex<long double>(__rho * ::cos(__phi), __rho * ::sin(__phi)); }
sl@0
    80
#endif
sl@0
    81
sl@0
    82
// Division
sl@0
    83
template <class _Tp>
sl@0
    84
static void _divT(const _Tp& __z1_r, const _Tp& __z1_i,
sl@0
    85
                  const _Tp& __z2_r, const _Tp& __z2_i,
sl@0
    86
                  _Tp& __res_r, _Tp& __res_i) {
sl@0
    87
  _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
sl@0
    88
  _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
sl@0
    89
sl@0
    90
  if (__ar <= __ai) {
sl@0
    91
    _Tp __ratio = __z2_r / __z2_i;
sl@0
    92
    _Tp __denom = __z2_i * (1 + __ratio * __ratio);
sl@0
    93
    __res_r = (__z1_r * __ratio + __z1_i) / __denom;
sl@0
    94
    __res_i = (__z1_i * __ratio - __z1_r) / __denom;
sl@0
    95
  }
sl@0
    96
  else {
sl@0
    97
    _Tp __ratio = __z2_i / __z2_r;
sl@0
    98
    _Tp __denom = __z2_r * (1 + __ratio * __ratio);
sl@0
    99
    __res_r = (__z1_r + __z1_i * __ratio) / __denom;
sl@0
   100
    __res_i = (__z1_i - __z1_r * __ratio) / __denom;
sl@0
   101
  }
sl@0
   102
}
sl@0
   103
sl@0
   104
template <class _Tp>
sl@0
   105
static void _divT(const _Tp& __z1_r,
sl@0
   106
                  const _Tp& __z2_r, const _Tp& __z2_i,
sl@0
   107
                  _Tp& __res_r, _Tp& __res_i) {
sl@0
   108
  _Tp __ar = __z2_r >= 0 ? __z2_r : -__z2_r;
sl@0
   109
  _Tp __ai = __z2_i >= 0 ? __z2_i : -__z2_i;
sl@0
   110
sl@0
   111
  if (__ar <= __ai) {
sl@0
   112
    _Tp __ratio = __z2_r / __z2_i;
sl@0
   113
    _Tp __denom = __z2_i * (1 + __ratio * __ratio);
sl@0
   114
    __res_r = (__z1_r * __ratio) / __denom;
sl@0
   115
    __res_i = - __z1_r / __denom;
sl@0
   116
  }
sl@0
   117
  else {
sl@0
   118
    _Tp __ratio = __z2_i / __z2_r;
sl@0
   119
    _Tp __denom = __z2_r * (1 + __ratio * __ratio);
sl@0
   120
    __res_r = __z1_r / __denom;
sl@0
   121
    __res_i = - (__z1_r * __ratio) / __denom;
sl@0
   122
  }
sl@0
   123
}
sl@0
   124
sl@0
   125
_STLP_DECLSPEC void _STLP_CALL
sl@0
   126
complex<float>::_div(const float& __z1_r, const float& __z1_i,
sl@0
   127
                     const float& __z2_r, const float& __z2_i,
sl@0
   128
                     float& __res_r, float& __res_i)
sl@0
   129
{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
sl@0
   130
sl@0
   131
_STLP_DECLSPEC void _STLP_CALL
sl@0
   132
complex<float>::_div(const float& __z1_r,
sl@0
   133
                     const float& __z2_r, const float& __z2_i,
sl@0
   134
                     float& __res_r, float& __res_i)
sl@0
   135
{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
sl@0
   136
sl@0
   137
sl@0
   138
_STLP_DECLSPEC void _STLP_CALL
sl@0
   139
complex<double>::_div(const double& __z1_r, const double& __z1_i,
sl@0
   140
                      const double& __z2_r, const double& __z2_i,
sl@0
   141
                      double& __res_r, double& __res_i)
sl@0
   142
{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
sl@0
   143
sl@0
   144
_STLP_DECLSPEC void _STLP_CALL
sl@0
   145
complex<double>::_div(const double& __z1_r,
sl@0
   146
                      const double& __z2_r, const double& __z2_i,
sl@0
   147
                      double& __res_r, double& __res_i)
sl@0
   148
{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
sl@0
   149
sl@0
   150
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
   151
_STLP_DECLSPEC void _STLP_CALL
sl@0
   152
complex<long double>::_div(const long double& __z1_r, const long double& __z1_i,
sl@0
   153
                           const long double& __z2_r, const long double& __z2_i,
sl@0
   154
                           long double& __res_r, long double& __res_i)
sl@0
   155
{ _divT(__z1_r, __z1_i, __z2_r, __z2_i, __res_r, __res_i); }
sl@0
   156
sl@0
   157
_STLP_DECLSPEC void _STLP_CALL
sl@0
   158
complex<long double>::_div(const long double& __z1_r,
sl@0
   159
                           const long double& __z2_r, const long double& __z2_i,
sl@0
   160
                           long double& __res_r, long double& __res_i)
sl@0
   161
{ _divT(__z1_r, __z2_r, __z2_i, __res_r, __res_i); }
sl@0
   162
#endif
sl@0
   163
sl@0
   164
//----------------------------------------------------------------------
sl@0
   165
// Square root
sl@0
   166
template <class _Tp>
sl@0
   167
static complex<_Tp> sqrtT(const complex<_Tp>& z) {
sl@0
   168
  _Tp re = z._M_re;
sl@0
   169
  _Tp im = z._M_im;
sl@0
   170
  _Tp mag = ::hypot(re, im);
sl@0
   171
  complex<_Tp> result;
sl@0
   172
sl@0
   173
  if (mag == 0.f) {
sl@0
   174
    result._M_re = result._M_im = 0.f;
sl@0
   175
  } else if (re > 0.f) {
sl@0
   176
    result._M_re = ::sqrt(0.5f * (mag + re));
sl@0
   177
    result._M_im = im/result._M_re/2.f;
sl@0
   178
  } else {
sl@0
   179
    result._M_im = ::sqrt(0.5f * (mag - re));
sl@0
   180
    if (im < 0.f)
sl@0
   181
      result._M_im = - result._M_im;
sl@0
   182
    result._M_re = im/result._M_im/2.f;
sl@0
   183
  }
sl@0
   184
  return result;
sl@0
   185
}
sl@0
   186
sl@0
   187
_STLP_DECLSPEC complex<float> _STLP_CALL
sl@0
   188
sqrt(const complex<float>& z) { return sqrtT(z); }
sl@0
   189
sl@0
   190
_STLP_DECLSPEC complex<double>  _STLP_CALL
sl@0
   191
sqrt(const complex<double>& z) { return sqrtT(z); }
sl@0
   192
sl@0
   193
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
   194
_STLP_DECLSPEC complex<long double> _STLP_CALL
sl@0
   195
sqrt(const complex<long double>& z) { return sqrtT(z); }
sl@0
   196
#endif
sl@0
   197
sl@0
   198
// exp, log, pow for complex<float>, complex<double>, and complex<long double>
sl@0
   199
//----------------------------------------------------------------------
sl@0
   200
// exp
sl@0
   201
template <class _Tp>
sl@0
   202
static complex<_Tp> expT(const complex<_Tp>& z) {
sl@0
   203
  _Tp expx = ::exp(z._M_re);
sl@0
   204
  return complex<_Tp>(expx * ::cos(z._M_im),
sl@0
   205
                      expx * ::sin(z._M_im));
sl@0
   206
}
sl@0
   207
_STLP_DECLSPEC complex<float>  _STLP_CALL exp(const complex<float>& z)
sl@0
   208
{ return expT(z); }
sl@0
   209
sl@0
   210
_STLP_DECLSPEC complex<double> _STLP_CALL exp(const complex<double>& z)
sl@0
   211
{ return expT(z); }
sl@0
   212
sl@0
   213
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
   214
_STLP_DECLSPEC complex<long double> _STLP_CALL exp(const complex<long double>& z)
sl@0
   215
{ return expT(z); }
sl@0
   216
#endif
sl@0
   217
sl@0
   218
//----------------------------------------------------------------------
sl@0
   219
// log10
sl@0
   220
template <class _Tp>
sl@0
   221
static complex<_Tp> log10T(const complex<_Tp>& z, const _Tp& ln10_inv) {
sl@0
   222
  complex<_Tp> r;
sl@0
   223
sl@0
   224
  r._M_im = ::atan2(z._M_im, z._M_re) * ln10_inv;
sl@0
   225
  r._M_re = ::log10(::hypot(z._M_re, z._M_im));
sl@0
   226
  return r;
sl@0
   227
}
sl@0
   228
sl@0
   229
static const float LN10_INVF = 1.f / ::log(10.f);
sl@0
   230
_STLP_DECLSPEC complex<float> _STLP_CALL log10(const complex<float>& z)
sl@0
   231
{ return log10T(z, LN10_INVF); }
sl@0
   232
sl@0
   233
static const double LN10_INV = 1. / ::log10(10.);
sl@0
   234
_STLP_DECLSPEC complex<double> _STLP_CALL log10(const complex<double>& z)
sl@0
   235
{ return log10T(z, LN10_INV); }
sl@0
   236
sl@0
   237
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
   238
static const long double LN10_INVL = 1.l / ::log(10.l);
sl@0
   239
_STLP_DECLSPEC complex<long double> _STLP_CALL log10(const complex<long double>& z)
sl@0
   240
{ return log10T(z, LN10_INVL); }
sl@0
   241
#endif
sl@0
   242
sl@0
   243
//----------------------------------------------------------------------
sl@0
   244
// log
sl@0
   245
template <class _Tp>
sl@0
   246
static complex<_Tp> logT(const complex<_Tp>& z) {
sl@0
   247
  complex<_Tp> r;
sl@0
   248
sl@0
   249
  r._M_im = ::atan2(z._M_im, z._M_re);
sl@0
   250
  r._M_re = ::log(::hypot(z._M_re, z._M_im));
sl@0
   251
  return r;
sl@0
   252
}
sl@0
   253
_STLP_DECLSPEC complex<float> _STLP_CALL log(const complex<float>& z)
sl@0
   254
{ return logT(z); }
sl@0
   255
sl@0
   256
_STLP_DECLSPEC complex<double> _STLP_CALL log(const complex<double>& z)
sl@0
   257
{ return logT(z); }
sl@0
   258
sl@0
   259
#ifndef _STLP_NO_LONG_DOUBLE
sl@0
   260
_STLP_DECLSPEC complex<long double> _STLP_CALL log(const complex<long double>& z)
sl@0
   261
{ return logT(z); }
sl@0
   262
# endif
sl@0
   263
sl@0
   264
//----------------------------------------------------------------------
sl@0
   265
// pow
sl@0
   266
template <class _Tp>
sl@0
   267
static complex<_Tp> powT(const _Tp& a, const complex<_Tp>& b) {
sl@0
   268
  _Tp logr = ::log(a);
sl@0
   269
  _Tp x = ::exp(logr * b._M_re);
sl@0
   270
  _Tp y = logr * b._M_im;
sl@0
   271
sl@0
   272
  return complex<_Tp>(x * ::cos(y), x * ::sin(y));
sl@0
   273
}
sl@0
   274
sl@0
   275
template <class _Tp>
sl@0
   276
static complex<_Tp> powT(const complex<_Tp>& z_in, int n) {
sl@0
   277
  complex<_Tp> z = z_in;
sl@0
   278
  z = _STLP_PRIV __power(z, (n < 0 ? -n : n), multiplies< complex<_Tp> >());
sl@0
   279
  if (n < 0)
sl@0
   280
    return _Tp(1.0) / z;
sl@0
   281
  else
sl@0
   282
    return z;
sl@0
   283
}
sl@0
   284
sl@0
   285
template <class _Tp>
sl@0
   286
static complex<_Tp> powT(const complex<_Tp>& a, const _Tp& b) {
sl@0
   287
  _Tp logr = ::log(::hypot(a._M_re,a._M_im));
sl@0
   288
  _Tp logi = ::atan2(a._M_im, a._M_re);
sl@0
   289
  _Tp x = ::exp(logr * b);
sl@0
   290
  _Tp y = logi * b;
sl@0
   291
sl@0
   292
  return complex<_Tp>(x * ::cos(y), x * ::sin(y));
sl@0
   293
}
sl@0
   294
sl@0
   295
template <class _Tp>
sl@0
   296
static complex<_Tp> powT(const complex<_Tp>& a, const complex<_Tp>& b) {
sl@0
   297
  _Tp logr = ::log(::hypot(a._M_re,a._M_im));
sl@0
   298
  _Tp logi = ::atan2(a._M_im, a._M_re);
sl@0
   299
  _Tp x = ::exp(logr * b._M_re - logi * b._M_im);
sl@0
   300
  _Tp y = logr * b._M_im + logi * b._M_re;
sl@0
   301
sl@0
   302
  return complex<_Tp>(x * ::cos(y), x * ::sin(y));
sl@0
   303
}
sl@0
   304
sl@0
   305
_STLP_DECLSPEC complex<float> _STLP_CALL pow(const float& a, const complex<float>& b)
sl@0
   306
{ return powT(a, b); }
sl@0
   307
sl@0
   308
_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& z_in, int n)
sl@0
   309
{ return powT(z_in, n); }
sl@0
   310
sl@0
   311
_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const float& b)
sl@0
   312
{ return powT(a, b); }
sl@0
   313
sl@0
   314
_STLP_DECLSPEC complex<float> _STLP_CALL pow(const complex<float>& a, const complex<float>& b)
sl@0
   315
{ return powT(a, b); }
sl@0
   316
sl@0
   317
_STLP_DECLSPEC complex<double> _STLP_CALL pow(const double& a, const complex<double>& b)
sl@0
   318
{ return powT(a, b); }
sl@0
   319
sl@0
   320
_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& z_in, int n)
sl@0
   321
{ return powT(z_in, n); }
sl@0
   322
sl@0
   323
_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const double& b)
sl@0
   324
{ return powT(a, b); }
sl@0
   325
sl@0
   326
_STLP_DECLSPEC complex<double> _STLP_CALL pow(const complex<double>& a, const complex<double>& b)
sl@0
   327
{ return powT(a, b); }
sl@0
   328
sl@0
   329
#if !defined (_STLP_NO_LONG_DOUBLE)
sl@0
   330
_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const long double& a,
sl@0
   331
                                                   const complex<long double>& b)
sl@0
   332
{ return powT(a, b); }
sl@0
   333
sl@0
   334
sl@0
   335
_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& z_in, int n)
sl@0
   336
{ return powT(z_in, n); }
sl@0
   337
sl@0
   338
_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
sl@0
   339
                                                   const long double& b)
sl@0
   340
{ return powT(a, b); }
sl@0
   341
sl@0
   342
_STLP_DECLSPEC complex<long double> _STLP_CALL pow(const complex<long double>& a,
sl@0
   343
                                                   const complex<long double>& b)
sl@0
   344
{ return powT(a, b); }
sl@0
   345
#endif
sl@0
   346
sl@0
   347
_STLP_END_NAMESPACE