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>
sinT(const complex<_Tp> & z)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
sin(const complex<float> & z)73 _STLP_DECLSPEC complex<float> _STLP_CALL sin(const complex<float>& z)
74 { return sinT(z); }
75
sin(const complex<double> & z)76 _STLP_DECLSPEC complex<double> _STLP_CALL sin(const complex<double>& z)
77 { return sinT(z); }
78
79 #if !defined (_STLP_NO_LONG_DOUBLE)
sin(const complex<long double> & z)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>
cosT(const complex<_Tp> & z)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
cos(const complex<float> & z)92 _STLP_DECLSPEC complex<float> _STLP_CALL cos(const complex<float>& z)
93 { return cosT(z); }
94
cos(const complex<double> & z)95 _STLP_DECLSPEC complex<double> _STLP_CALL cos(const complex<double>& z)
96 { return cosT(z); }
97
98 #if !defined (_STLP_NO_LONG_DOUBLE)
cos(const complex<long double> & z)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>
tanT(const complex<_Tp> & z,const _Tp & Tp_limit)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
tan(const complex<float> & z)118 _STLP_DECLSPEC complex<float> _STLP_CALL tan(const complex<float>& z)
119 { return tanT(z, float_limit); }
120
tan(const complex<double> & z)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)
tan(const complex<long double> & z)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>
sinhT(const complex<_Tp> & z)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
sinh(const complex<float> & z)137 _STLP_DECLSPEC complex<float> _STLP_CALL sinh(const complex<float>& z)
138 { return sinhT(z); }
139
sinh(const complex<double> & z)140 _STLP_DECLSPEC complex<double> _STLP_CALL sinh(const complex<double>& z)
141 { return sinhT(z); }
142
143 #if !defined (_STLP_NO_LONG_DOUBLE)
sinh(const complex<long double> & z)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>
coshT(const complex<_Tp> & z)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
cosh(const complex<float> & z)156 _STLP_DECLSPEC complex<float> _STLP_CALL cosh(const complex<float>& z)
157 { return coshT(z); }
158
cosh(const complex<double> & z)159 _STLP_DECLSPEC complex<double> _STLP_CALL cosh(const complex<double>& z)
160 { return coshT(z); }
161
162 #if !defined (_STLP_NO_LONG_DOUBLE)
cosh(const complex<long double> & z)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>
tanhT(const complex<_Tp> & z,const _Tp & Tp_limit)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
tanh(const complex<float> & z)181 _STLP_DECLSPEC complex<float> _STLP_CALL tanh(const complex<float>& z)
182 { return tanhT(z, float_limit); }
183
tanh(const complex<double> & z)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)
tanh(const complex<long double> & z)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
193