williamr@2
|
1 |
// (C) Copyright John Maddock 2005.
|
williamr@2
|
2 |
// Distributed under the Boost Software License, Version 1.0. (See accompanying
|
williamr@2
|
3 |
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
williamr@2
|
4 |
|
williamr@2
|
5 |
#ifndef BOOST_MATH_COMPLEX_ASIN_INCLUDED
|
williamr@2
|
6 |
#define BOOST_MATH_COMPLEX_ASIN_INCLUDED
|
williamr@2
|
7 |
|
williamr@2
|
8 |
#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
|
williamr@2
|
9 |
# include <boost/math/complex/details.hpp>
|
williamr@2
|
10 |
#endif
|
williamr@2
|
11 |
#ifndef BOOST_MATH_LOG1P_INCLUDED
|
williamr@2
|
12 |
# include <boost/math/special_functions/log1p.hpp>
|
williamr@2
|
13 |
#endif
|
williamr@2
|
14 |
#include <boost/assert.hpp>
|
williamr@2
|
15 |
|
williamr@2
|
16 |
#ifdef BOOST_NO_STDC_NAMESPACE
|
williamr@2
|
17 |
namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
|
williamr@2
|
18 |
#endif
|
williamr@2
|
19 |
|
williamr@2
|
20 |
namespace boost{ namespace math{
|
williamr@2
|
21 |
|
williamr@2
|
22 |
template<class T>
|
williamr@2
|
23 |
inline std::complex<T> asin(const std::complex<T>& z)
|
williamr@2
|
24 |
{
|
williamr@2
|
25 |
//
|
williamr@2
|
26 |
// This implementation is a transcription of the pseudo-code in:
|
williamr@2
|
27 |
//
|
williamr@2
|
28 |
// "Implementing the complex Arcsine and Arccosine Functions using Exception Handling."
|
williamr@2
|
29 |
// T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
|
williamr@2
|
30 |
// ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
|
williamr@2
|
31 |
//
|
williamr@2
|
32 |
|
williamr@2
|
33 |
//
|
williamr@2
|
34 |
// These static constants should really be in a maths constants library:
|
williamr@2
|
35 |
//
|
williamr@2
|
36 |
static const T one = static_cast<T>(1);
|
williamr@2
|
37 |
//static const T two = static_cast<T>(2);
|
williamr@2
|
38 |
static const T half = static_cast<T>(0.5L);
|
williamr@2
|
39 |
static const T a_crossover = static_cast<T>(1.5L);
|
williamr@2
|
40 |
static const T b_crossover = static_cast<T>(0.6417L);
|
williamr@2
|
41 |
//static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
|
williamr@2
|
42 |
static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
|
williamr@2
|
43 |
static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);
|
williamr@2
|
44 |
static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);
|
williamr@2
|
45 |
|
williamr@2
|
46 |
//
|
williamr@2
|
47 |
// Get real and imaginary parts, discard the signs as we can
|
williamr@2
|
48 |
// figure out the sign of the result later:
|
williamr@2
|
49 |
//
|
williamr@2
|
50 |
T x = std::fabs(z.real());
|
williamr@2
|
51 |
T y = std::fabs(z.imag());
|
williamr@2
|
52 |
T real, imag; // our results
|
williamr@2
|
53 |
|
williamr@2
|
54 |
//
|
williamr@2
|
55 |
// Begin by handling the special cases for infinities and nan's
|
williamr@2
|
56 |
// specified in C99, most of this is handled by the regular logic
|
williamr@2
|
57 |
// below, but handling it as a special case prevents overflow/underflow
|
williamr@2
|
58 |
// arithmetic which may trip up some machines:
|
williamr@2
|
59 |
//
|
williamr@2
|
60 |
if(detail::test_is_nan(x))
|
williamr@2
|
61 |
{
|
williamr@2
|
62 |
if(detail::test_is_nan(y))
|
williamr@2
|
63 |
return std::complex<T>(x, x);
|
williamr@2
|
64 |
if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
|
williamr@2
|
65 |
{
|
williamr@2
|
66 |
real = x;
|
williamr@2
|
67 |
imag = std::numeric_limits<T>::infinity();
|
williamr@2
|
68 |
}
|
williamr@2
|
69 |
else
|
williamr@2
|
70 |
return std::complex<T>(x, x);
|
williamr@2
|
71 |
}
|
williamr@2
|
72 |
else if(detail::test_is_nan(y))
|
williamr@2
|
73 |
{
|
williamr@2
|
74 |
if(x == 0)
|
williamr@2
|
75 |
{
|
williamr@2
|
76 |
real = 0;
|
williamr@2
|
77 |
imag = y;
|
williamr@2
|
78 |
}
|
williamr@2
|
79 |
else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
|
williamr@2
|
80 |
{
|
williamr@2
|
81 |
real = y;
|
williamr@2
|
82 |
imag = std::numeric_limits<T>::infinity();
|
williamr@2
|
83 |
}
|
williamr@2
|
84 |
else
|
williamr@2
|
85 |
return std::complex<T>(y, y);
|
williamr@2
|
86 |
}
|
williamr@2
|
87 |
else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
|
williamr@2
|
88 |
{
|
williamr@2
|
89 |
if(y == std::numeric_limits<T>::infinity())
|
williamr@2
|
90 |
{
|
williamr@2
|
91 |
real = quarter_pi;
|
williamr@2
|
92 |
imag = std::numeric_limits<T>::infinity();
|
williamr@2
|
93 |
}
|
williamr@2
|
94 |
else
|
williamr@2
|
95 |
{
|
williamr@2
|
96 |
real = half_pi;
|
williamr@2
|
97 |
imag = std::numeric_limits<T>::infinity();
|
williamr@2
|
98 |
}
|
williamr@2
|
99 |
}
|
williamr@2
|
100 |
else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
|
williamr@2
|
101 |
{
|
williamr@2
|
102 |
real = 0;
|
williamr@2
|
103 |
imag = std::numeric_limits<T>::infinity();
|
williamr@2
|
104 |
}
|
williamr@2
|
105 |
else
|
williamr@2
|
106 |
{
|
williamr@2
|
107 |
//
|
williamr@2
|
108 |
// special case for real numbers:
|
williamr@2
|
109 |
//
|
williamr@2
|
110 |
if((y == 0) && (x <= one))
|
williamr@2
|
111 |
return std::complex<T>(std::asin(z.real()));
|
williamr@2
|
112 |
//
|
williamr@2
|
113 |
// Figure out if our input is within the "safe area" identified by Hull et al.
|
williamr@2
|
114 |
// This would be more efficient with portable floating point exception handling;
|
williamr@2
|
115 |
// fortunately the quantities M and u identified by Hull et al (figure 3),
|
williamr@2
|
116 |
// match with the max and min methods of numeric_limits<T>.
|
williamr@2
|
117 |
//
|
williamr@2
|
118 |
T safe_max = detail::safe_max(static_cast<T>(8));
|
williamr@2
|
119 |
T safe_min = detail::safe_min(static_cast<T>(4));
|
williamr@2
|
120 |
|
williamr@2
|
121 |
T xp1 = one + x;
|
williamr@2
|
122 |
T xm1 = x - one;
|
williamr@2
|
123 |
|
williamr@2
|
124 |
if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
|
williamr@2
|
125 |
{
|
williamr@2
|
126 |
T yy = y * y;
|
williamr@2
|
127 |
T r = std::sqrt(xp1*xp1 + yy);
|
williamr@2
|
128 |
T s = std::sqrt(xm1*xm1 + yy);
|
williamr@2
|
129 |
T a = half * (r + s);
|
williamr@2
|
130 |
T b = x / a;
|
williamr@2
|
131 |
|
williamr@2
|
132 |
if(b <= b_crossover)
|
williamr@2
|
133 |
{
|
williamr@2
|
134 |
real = std::asin(b);
|
williamr@2
|
135 |
}
|
williamr@2
|
136 |
else
|
williamr@2
|
137 |
{
|
williamr@2
|
138 |
T apx = a + x;
|
williamr@2
|
139 |
if(x <= one)
|
williamr@2
|
140 |
{
|
williamr@2
|
141 |
real = std::atan(x/std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1))));
|
williamr@2
|
142 |
}
|
williamr@2
|
143 |
else
|
williamr@2
|
144 |
{
|
williamr@2
|
145 |
real = std::atan(x/(y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1)))));
|
williamr@2
|
146 |
}
|
williamr@2
|
147 |
}
|
williamr@2
|
148 |
|
williamr@2
|
149 |
if(a <= a_crossover)
|
williamr@2
|
150 |
{
|
williamr@2
|
151 |
T am1;
|
williamr@2
|
152 |
if(x < one)
|
williamr@2
|
153 |
{
|
williamr@2
|
154 |
am1 = half * (yy/(r + xp1) + yy/(s - xm1));
|
williamr@2
|
155 |
}
|
williamr@2
|
156 |
else
|
williamr@2
|
157 |
{
|
williamr@2
|
158 |
am1 = half * (yy/(r + xp1) + (s + xm1));
|
williamr@2
|
159 |
}
|
williamr@2
|
160 |
imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
|
williamr@2
|
161 |
}
|
williamr@2
|
162 |
else
|
williamr@2
|
163 |
{
|
williamr@2
|
164 |
imag = std::log(a + std::sqrt(a*a - one));
|
williamr@2
|
165 |
}
|
williamr@2
|
166 |
}
|
williamr@2
|
167 |
else
|
williamr@2
|
168 |
{
|
williamr@2
|
169 |
//
|
williamr@2
|
170 |
// This is the Hull et al exception handling code from Fig 3 of their paper:
|
williamr@2
|
171 |
//
|
williamr@2
|
172 |
if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
|
williamr@2
|
173 |
{
|
williamr@2
|
174 |
if(x < one)
|
williamr@2
|
175 |
{
|
williamr@2
|
176 |
real = std::asin(x);
|
williamr@2
|
177 |
imag = y / std::sqrt(xp1*xm1);
|
williamr@2
|
178 |
}
|
williamr@2
|
179 |
else
|
williamr@2
|
180 |
{
|
williamr@2
|
181 |
real = half_pi;
|
williamr@2
|
182 |
if(((std::numeric_limits<T>::max)() / xp1) > xm1)
|
williamr@2
|
183 |
{
|
williamr@2
|
184 |
// xp1 * xm1 won't overflow:
|
williamr@2
|
185 |
imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
|
williamr@2
|
186 |
}
|
williamr@2
|
187 |
else
|
williamr@2
|
188 |
{
|
williamr@2
|
189 |
imag = log_two + std::log(x);
|
williamr@2
|
190 |
}
|
williamr@2
|
191 |
}
|
williamr@2
|
192 |
}
|
williamr@2
|
193 |
else if(y <= safe_min)
|
williamr@2
|
194 |
{
|
williamr@2
|
195 |
// There is an assumption in Hull et al's analysis that
|
williamr@2
|
196 |
// if we get here then x == 1. This is true for all "good"
|
williamr@2
|
197 |
// machines where :
|
williamr@2
|
198 |
//
|
williamr@2
|
199 |
// E^2 > 8*sqrt(u); with:
|
williamr@2
|
200 |
//
|
williamr@2
|
201 |
// E = std::numeric_limits<T>::epsilon()
|
williamr@2
|
202 |
// u = (std::numeric_limits<T>::min)()
|
williamr@2
|
203 |
//
|
williamr@2
|
204 |
// Hull et al provide alternative code for "bad" machines
|
williamr@2
|
205 |
// but we have no way to test that here, so for now just assert
|
williamr@2
|
206 |
// on the assumption:
|
williamr@2
|
207 |
//
|
williamr@2
|
208 |
BOOST_ASSERT(x == 1);
|
williamr@2
|
209 |
real = half_pi - std::sqrt(y);
|
williamr@2
|
210 |
imag = std::sqrt(y);
|
williamr@2
|
211 |
}
|
williamr@2
|
212 |
else if(std::numeric_limits<T>::epsilon() * y - one >= x)
|
williamr@2
|
213 |
{
|
williamr@2
|
214 |
real = x/y; // This can underflow!
|
williamr@2
|
215 |
imag = log_two + std::log(y);
|
williamr@2
|
216 |
}
|
williamr@2
|
217 |
else if(x > one)
|
williamr@2
|
218 |
{
|
williamr@2
|
219 |
real = std::atan(x/y);
|
williamr@2
|
220 |
T xoy = x/y;
|
williamr@2
|
221 |
imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
|
williamr@2
|
222 |
}
|
williamr@2
|
223 |
else
|
williamr@2
|
224 |
{
|
williamr@2
|
225 |
T a = std::sqrt(one + y*y);
|
williamr@2
|
226 |
real = x/a; // This can underflow!
|
williamr@2
|
227 |
imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
|
williamr@2
|
228 |
}
|
williamr@2
|
229 |
}
|
williamr@2
|
230 |
}
|
williamr@2
|
231 |
|
williamr@2
|
232 |
//
|
williamr@2
|
233 |
// Finish off by working out the sign of the result:
|
williamr@2
|
234 |
//
|
williamr@2
|
235 |
if(z.real() < 0)
|
williamr@2
|
236 |
real = -real;
|
williamr@2
|
237 |
if(z.imag() < 0)
|
williamr@2
|
238 |
imag = -imag;
|
williamr@2
|
239 |
|
williamr@2
|
240 |
return std::complex<T>(real, imag);
|
williamr@2
|
241 |
}
|
williamr@2
|
242 |
|
williamr@2
|
243 |
} } // namespaces
|
williamr@2
|
244 |
|
williamr@2
|
245 |
#endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED
|