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