lib/c++/stlport/complex_trig.cpp

   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 #include "stlport_prefix.h"
  19
  20
  21 // Trigonometric and hyperbolic functions for complex<float>,
  22 // complex<double>, and complex<long double>
  23 #include <complex>
  24 #include <cfloat>
  25 #include <cmath>
  26
  27 _STLP_BEGIN_NAMESPACE
  28
  29
  30 //----------------------------------------------------------------------
  31 // helpers
  32 #if defined (__sgi)
  33   static const union { unsigned int i; float f; } float_ulimit = { 0x42b2d4fc };
  34   static const float float_limit = float_ulimit.f;
  35   static union {
  36     struct { unsigned int h; unsigned int l; } w;
  37     double d;
  38   } double_ulimit = { 0x408633ce, 0x8fb9f87d };
  39   static const double double_limit = double_ulimit.d;
  40   static union {
  41     struct { unsigned int h[2]; unsigned int l[2]; } w;
  42     long double ld;
  43   } ldouble_ulimit = {0x408633ce, 0x8fb9f87e, 0xbd23b659, 0x4e9bd8b1};
  44 #  if !defined (_STLP_NO_LONG_DOUBLE)
  45   static const long double ldouble_limit = ldouble_ulimit.ld;
  46 #  endif
  47 #else
  48 #  if defined (M_LN2) && defined (FLT_MAX_EXP)
  49   static const float float_limit = float(M_LN2 * FLT_MAX_EXP);
  50   static const double double_limit = M_LN2 * DBL_MAX_EXP;
  51 #  else
  52   static const float float_limit = ::log(FLT_MAX);
  53   static const double double_limit = ::log(DBL_MAX);
  54 #  endif
  55 #  if !defined (_STLP_NO_LONG_DOUBLE)
  56 #    if defined (M_LN2l)
  57   static const long double ldouble_limit = M_LN2l * LDBL_MAX_EXP;
  58 #    else
  59   static const long double ldouble_limit = ::log(LDBL_MAX);
  60 #    endif
  61 #  endif
  62 #endif
  63
  64
  65 //----------------------------------------------------------------------
  66 // sin
  67 template <class _Tp>
  68 static complex<_Tp> sinT(const complex<_Tp>& z) {
  69   return complex<_Tp>(::sin(z._M_re) * ::cosh(z._M_im),
  70                       ::cos(z._M_re) * ::sinh(z._M_im));
  71 }
  72
  73 _STLP_DECLSPEC complex<float> _STLP_CALL sin(const complex<float>& z)
  74 { return sinT(z); }
  75
  76 _STLP_DECLSPEC complex<double> _STLP_CALL sin(const complex<double>& z)
  77 { return sinT(z); }
  78
  79 #if !defined (_STLP_NO_LONG_DOUBLE)
  80 _STLP_DECLSPEC complex<long double> _STLP_CALL sin(const complex<long double>& z)
  81 { return sinT(z); }
  82 #endif
  83
  84 //----------------------------------------------------------------------
  85 // cos
  86 template <class _Tp>
  87 static complex<_Tp> cosT(const complex<_Tp>& z) {
  88   return complex<_Tp>(::cos(z._M_re) * ::cosh(z._M_im),
  89                      -::sin(z._M_re) * ::sinh(z._M_im));
  90 }
  91
  92 _STLP_DECLSPEC complex<float> _STLP_CALL cos(const complex<float>& z)
  93 { return cosT(z); }
  94
  95 _STLP_DECLSPEC complex<double> _STLP_CALL cos(const complex<double>& z)
  96 { return cosT(z); }
  97
  98 #if !defined (_STLP_NO_LONG_DOUBLE)
  99 _STLP_DECLSPEC complex<long double> _STLP_CALL cos(const complex<long double>& z)
 100 { return cosT(z); }
 101 #endif
 102
 103 //----------------------------------------------------------------------
 104 // tan
 105 template <class _Tp>
 106 static complex<_Tp> tanT(const complex<_Tp>& z, const _Tp& Tp_limit) {
 107   _Tp re2 = 2.f * z._M_re;
 108   _Tp im2 = 2.f * z._M_im;
 109
 110   if (::abs(im2) > Tp_limit)
 111     return complex<_Tp>(0.f, (im2 > 0 ? 1.f : -1.f));
 112   else {
 113     _Tp den = ::cos(re2) + ::cosh(im2);
 114     return complex<_Tp>(::sin(re2) / den, ::sinh(im2) / den);
 115   }
 116 }
 117
 118 _STLP_DECLSPEC complex<float> _STLP_CALL tan(const complex<float>& z)
 119 { return tanT(z, float_limit); }
 120
 121 _STLP_DECLSPEC complex<double> _STLP_CALL tan(const complex<double>& z)
 122 { return tanT(z, double_limit); }
 123
 124 #if !defined (_STLP_NO_LONG_DOUBLE)
 125 _STLP_DECLSPEC complex<long double> _STLP_CALL tan(const complex<long double>& z)
 126 { return tanT(z, ldouble_limit); }
 127 #endif
 128
 129 //----------------------------------------------------------------------
 130 // sinh
 131 template <class _Tp>
 132 static complex<_Tp> sinhT(const complex<_Tp>& z) {
 133   return complex<_Tp>(::sinh(z._M_re) * ::cos(z._M_im),
 134                       ::cosh(z._M_re) * ::sin(z._M_im));
 135 }
 136
 137 _STLP_DECLSPEC complex<float> _STLP_CALL sinh(const complex<float>& z)
 138 { return sinhT(z); }
 139
 140 _STLP_DECLSPEC complex<double> _STLP_CALL sinh(const complex<double>& z)
 141 { return sinhT(z); }
 142
 143 #if !defined (_STLP_NO_LONG_DOUBLE)
 144 _STLP_DECLSPEC complex<long double> _STLP_CALL sinh(const complex<long double>& z)
 145 { return sinhT(z); }
 146 #endif
 147
 148 //----------------------------------------------------------------------
 149 // cosh
 150 template <class _Tp>
 151 static complex<_Tp> coshT(const complex<_Tp>& z) {
 152   return complex<_Tp>(::cosh(z._M_re) * ::cos(z._M_im),
 153                       ::sinh(z._M_re) * ::sin(z._M_im));
 154 }
 155
 156 _STLP_DECLSPEC complex<float> _STLP_CALL cosh(const complex<float>& z)
 157 { return coshT(z); }
 158
 159 _STLP_DECLSPEC complex<double> _STLP_CALL cosh(const complex<double>& z)
 160 { return coshT(z); }
 161
 162 #if !defined (_STLP_NO_LONG_DOUBLE)
 163 _STLP_DECLSPEC complex<long double> _STLP_CALL cosh(const complex<long double>& z)
 164 { return coshT(z); }
 165 #endif
 166
 167 //----------------------------------------------------------------------
 168 // tanh
 169 template <class _Tp>
 170 static complex<_Tp> tanhT(const complex<_Tp>& z, const _Tp& Tp_limit) {
 171   _Tp re2 = 2.f * z._M_re;
 172   _Tp im2 = 2.f * z._M_im;
 173   if (::abs(re2) > Tp_limit)
 174     return complex<_Tp>((re2 > 0 ? 1.f : -1.f), 0.f);
 175   else {
 176     _Tp den = ::cosh(re2) + ::cos(im2);
 177     return complex<_Tp>(::sinh(re2) / den, ::sin(im2) / den);
 178   }
 179 }
 180
 181 _STLP_DECLSPEC complex<float> _STLP_CALL tanh(const complex<float>& z)
 182 { return tanhT(z, float_limit); }
 183
 184 _STLP_DECLSPEC complex<double> _STLP_CALL tanh(const complex<double>& z)
 185 { return tanhT(z, double_limit); }
 186
 187 #if !defined (_STLP_NO_LONG_DOUBLE)
 188 _STLP_DECLSPEC complex<long double> _STLP_CALL tanh(const complex<long double>& z)
 189 { return tanhT(z, ldouble_limit); }
 190 #endif
 191
 192 _STLP_END_NAMESPACE