|
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 |