Chris@16
|
1
|
Chris@16
|
2 // Copyright Christopher Kormanyos 2002 - 2011.
|
Chris@16
|
3 // Copyright 2011 John Maddock. Distributed under the Boost
|
Chris@16
|
4 // Distributed under the Boost Software License, Version 1.0.
|
Chris@16
|
5 // (See accompanying file LICENSE_1_0.txt or copy at
|
Chris@16
|
6 // http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
7
|
Chris@16
|
8 // This work is based on an earlier work:
|
Chris@16
|
9 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
|
Chris@16
|
10 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
|
Chris@16
|
11 //
|
Chris@16
|
12 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
|
Chris@16
|
13 //
|
Chris@16
|
14
|
Chris@16
|
15 template <class T>
|
Chris@16
|
16 void hyp0F1(T& result, const T& b, const T& x)
|
Chris@16
|
17 {
|
Chris@16
|
18 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
Chris@16
|
19 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
20
|
Chris@16
|
21 // Compute the series representation of Hypergeometric0F1 taken from
|
Chris@16
|
22 // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/06/01/01/
|
Chris@16
|
23 // There are no checks on input range or parameter boundaries.
|
Chris@16
|
24
|
Chris@16
|
25 T x_pow_n_div_n_fact(x);
|
Chris@16
|
26 T pochham_b (b);
|
Chris@16
|
27 T bp (b);
|
Chris@16
|
28
|
Chris@16
|
29 eval_divide(result, x_pow_n_div_n_fact, pochham_b);
|
Chris@16
|
30 eval_add(result, ui_type(1));
|
Chris@16
|
31
|
Chris@16
|
32 si_type n;
|
Chris@16
|
33
|
Chris@16
|
34 T tol;
|
Chris@16
|
35 tol = ui_type(1);
|
Chris@16
|
36 eval_ldexp(tol, tol, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
|
Chris@16
|
37 eval_multiply(tol, result);
|
Chris@16
|
38 if(eval_get_sign(tol) < 0)
|
Chris@16
|
39 tol.negate();
|
Chris@16
|
40 T term;
|
Chris@16
|
41
|
Chris@101
|
42 static const int series_limit =
|
Chris@16
|
43 boost::multiprecision::detail::digits2<number<T, et_on> >::value < 100
|
Chris@16
|
44 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value;
|
Chris@16
|
45 // Series expansion of hyperg_0f1(; b; x).
|
Chris@16
|
46 for(n = 2; n < series_limit; ++n)
|
Chris@16
|
47 {
|
Chris@16
|
48 eval_multiply(x_pow_n_div_n_fact, x);
|
Chris@16
|
49 eval_divide(x_pow_n_div_n_fact, n);
|
Chris@16
|
50 eval_increment(bp);
|
Chris@16
|
51 eval_multiply(pochham_b, bp);
|
Chris@16
|
52
|
Chris@16
|
53 eval_divide(term, x_pow_n_div_n_fact, pochham_b);
|
Chris@16
|
54 eval_add(result, term);
|
Chris@16
|
55
|
Chris@16
|
56 bool neg_term = eval_get_sign(term) < 0;
|
Chris@16
|
57 if(neg_term)
|
Chris@16
|
58 term.negate();
|
Chris@16
|
59 if(term.compare(tol) <= 0)
|
Chris@16
|
60 break;
|
Chris@16
|
61 }
|
Chris@16
|
62
|
Chris@16
|
63 if(n >= series_limit)
|
Chris@16
|
64 BOOST_THROW_EXCEPTION(std::runtime_error("H0F1 Failed to Converge"));
|
Chris@16
|
65 }
|
Chris@16
|
66
|
Chris@16
|
67
|
Chris@16
|
68 template <class T>
|
Chris@16
|
69 void eval_sin(T& result, const T& x)
|
Chris@16
|
70 {
|
Chris@16
|
71 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sin function is only valid for floating point types.");
|
Chris@16
|
72 if(&result == &x)
|
Chris@16
|
73 {
|
Chris@16
|
74 T temp;
|
Chris@16
|
75 eval_sin(temp, x);
|
Chris@16
|
76 result = temp;
|
Chris@16
|
77 return;
|
Chris@16
|
78 }
|
Chris@16
|
79
|
Chris@16
|
80 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
Chris@16
|
81 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
82 typedef typename mpl::front<typename T::float_types>::type fp_type;
|
Chris@16
|
83
|
Chris@16
|
84 switch(eval_fpclassify(x))
|
Chris@16
|
85 {
|
Chris@16
|
86 case FP_INFINITE:
|
Chris@16
|
87 case FP_NAN:
|
Chris@16
|
88 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
89 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
90 else
|
Chris@16
|
91 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
92 return;
|
Chris@16
|
93 case FP_ZERO:
|
Chris@16
|
94 result = ui_type(0);
|
Chris@16
|
95 return;
|
Chris@16
|
96 default: ;
|
Chris@16
|
97 }
|
Chris@16
|
98
|
Chris@16
|
99 // Local copy of the argument
|
Chris@16
|
100 T xx = x;
|
Chris@16
|
101
|
Chris@16
|
102 // Analyze and prepare the phase of the argument.
|
Chris@16
|
103 // Make a local, positive copy of the argument, xx.
|
Chris@16
|
104 // The argument xx will be reduced to 0 <= xx <= pi/2.
|
Chris@16
|
105 bool b_negate_sin = false;
|
Chris@16
|
106
|
Chris@16
|
107 if(eval_get_sign(x) < 0)
|
Chris@16
|
108 {
|
Chris@16
|
109 xx.negate();
|
Chris@16
|
110 b_negate_sin = !b_negate_sin;
|
Chris@16
|
111 }
|
Chris@16
|
112
|
Chris@16
|
113 T n_pi, t;
|
Chris@16
|
114 // Remove even multiples of pi.
|
Chris@16
|
115 if(xx.compare(get_constant_pi<T>()) > 0)
|
Chris@16
|
116 {
|
Chris@16
|
117 eval_divide(n_pi, xx, get_constant_pi<T>());
|
Chris@16
|
118 eval_trunc(n_pi, n_pi);
|
Chris@16
|
119 t = ui_type(2);
|
Chris@16
|
120 eval_fmod(t, n_pi, t);
|
Chris@16
|
121 const bool b_n_pi_is_even = eval_get_sign(t) == 0;
|
Chris@16
|
122 eval_multiply(n_pi, get_constant_pi<T>());
|
Chris@16
|
123 eval_subtract(xx, n_pi);
|
Chris@16
|
124
|
Chris@16
|
125 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
Chris@16
|
126 BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
|
Chris@16
|
127
|
Chris@16
|
128 // Adjust signs if the multiple of pi is not even.
|
Chris@16
|
129 if(!b_n_pi_is_even)
|
Chris@16
|
130 {
|
Chris@16
|
131 b_negate_sin = !b_negate_sin;
|
Chris@16
|
132 }
|
Chris@16
|
133 }
|
Chris@16
|
134
|
Chris@16
|
135 // Reduce the argument to 0 <= xx <= pi/2.
|
Chris@16
|
136 eval_ldexp(t, get_constant_pi<T>(), -1);
|
Chris@16
|
137 if(xx.compare(t) > 0)
|
Chris@16
|
138 {
|
Chris@16
|
139 eval_subtract(xx, get_constant_pi<T>(), xx);
|
Chris@16
|
140 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
Chris@16
|
141 }
|
Chris@16
|
142
|
Chris@16
|
143 eval_subtract(t, xx);
|
Chris@16
|
144 const bool b_zero = eval_get_sign(xx) == 0;
|
Chris@16
|
145 const bool b_pi_half = eval_get_sign(t) == 0;
|
Chris@16
|
146
|
Chris@16
|
147 // Check if the reduced argument is very close to 0 or pi/2.
|
Chris@16
|
148 const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
|
Chris@16
|
149 const bool b_near_pi_half = t.compare(fp_type(1e-1)) < 0;;
|
Chris@16
|
150
|
Chris@16
|
151 if(b_zero)
|
Chris@16
|
152 {
|
Chris@16
|
153 result = ui_type(0);
|
Chris@16
|
154 }
|
Chris@16
|
155 else if(b_pi_half)
|
Chris@16
|
156 {
|
Chris@16
|
157 result = ui_type(1);
|
Chris@16
|
158 }
|
Chris@16
|
159 else if(b_near_zero)
|
Chris@16
|
160 {
|
Chris@16
|
161 eval_multiply(t, xx, xx);
|
Chris@16
|
162 eval_divide(t, si_type(-4));
|
Chris@16
|
163 T t2;
|
Chris@16
|
164 t2 = fp_type(1.5);
|
Chris@16
|
165 hyp0F1(result, t2, t);
|
Chris@16
|
166 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
Chris@16
|
167 eval_multiply(result, xx);
|
Chris@16
|
168 }
|
Chris@16
|
169 else if(b_near_pi_half)
|
Chris@16
|
170 {
|
Chris@16
|
171 eval_multiply(t, t);
|
Chris@16
|
172 eval_divide(t, si_type(-4));
|
Chris@16
|
173 T t2;
|
Chris@16
|
174 t2 = fp_type(0.5);
|
Chris@16
|
175 hyp0F1(result, t2, t);
|
Chris@16
|
176 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
Chris@16
|
177 }
|
Chris@16
|
178 else
|
Chris@16
|
179 {
|
Chris@16
|
180 // Scale to a small argument for an efficient Taylor series,
|
Chris@16
|
181 // implemented as a hypergeometric function. Use a standard
|
Chris@16
|
182 // divide by three identity a certain number of times.
|
Chris@16
|
183 // Here we use division by 3^9 --> (19683 = 3^9).
|
Chris@16
|
184
|
Chris@16
|
185 static const si_type n_scale = 9;
|
Chris@16
|
186 static const si_type n_three_pow_scale = static_cast<si_type>(19683L);
|
Chris@16
|
187
|
Chris@16
|
188 eval_divide(xx, n_three_pow_scale);
|
Chris@16
|
189
|
Chris@16
|
190 // Now with small arguments, we are ready for a series expansion.
|
Chris@16
|
191 eval_multiply(t, xx, xx);
|
Chris@16
|
192 eval_divide(t, si_type(-4));
|
Chris@16
|
193 T t2;
|
Chris@16
|
194 t2 = fp_type(1.5);
|
Chris@16
|
195 hyp0F1(result, t2, t);
|
Chris@16
|
196 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
Chris@16
|
197 eval_multiply(result, xx);
|
Chris@16
|
198
|
Chris@16
|
199 // Convert back using multiple angle identity.
|
Chris@16
|
200 for(boost::int32_t k = static_cast<boost::int32_t>(0); k < n_scale; k++)
|
Chris@16
|
201 {
|
Chris@16
|
202 // Rescale the cosine value using the multiple angle identity.
|
Chris@16
|
203 eval_multiply(t2, result, ui_type(3));
|
Chris@16
|
204 eval_multiply(t, result, result);
|
Chris@16
|
205 eval_multiply(t, result);
|
Chris@16
|
206 eval_multiply(t, ui_type(4));
|
Chris@16
|
207 eval_subtract(result, t2, t);
|
Chris@16
|
208 }
|
Chris@16
|
209 }
|
Chris@16
|
210
|
Chris@16
|
211 if(b_negate_sin)
|
Chris@16
|
212 result.negate();
|
Chris@16
|
213 }
|
Chris@16
|
214
|
Chris@16
|
215 template <class T>
|
Chris@16
|
216 void eval_cos(T& result, const T& x)
|
Chris@16
|
217 {
|
Chris@16
|
218 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cos function is only valid for floating point types.");
|
Chris@16
|
219 if(&result == &x)
|
Chris@16
|
220 {
|
Chris@16
|
221 T temp;
|
Chris@16
|
222 eval_cos(temp, x);
|
Chris@16
|
223 result = temp;
|
Chris@16
|
224 return;
|
Chris@16
|
225 }
|
Chris@16
|
226
|
Chris@16
|
227 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
Chris@16
|
228 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
229 typedef typename mpl::front<typename T::float_types>::type fp_type;
|
Chris@16
|
230
|
Chris@16
|
231 switch(eval_fpclassify(x))
|
Chris@16
|
232 {
|
Chris@16
|
233 case FP_INFINITE:
|
Chris@16
|
234 case FP_NAN:
|
Chris@16
|
235 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
236 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
237 else
|
Chris@16
|
238 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
239 return;
|
Chris@16
|
240 case FP_ZERO:
|
Chris@16
|
241 result = ui_type(1);
|
Chris@16
|
242 return;
|
Chris@16
|
243 default: ;
|
Chris@16
|
244 }
|
Chris@16
|
245
|
Chris@16
|
246 // Local copy of the argument
|
Chris@16
|
247 T xx = x;
|
Chris@16
|
248
|
Chris@16
|
249 // Analyze and prepare the phase of the argument.
|
Chris@16
|
250 // Make a local, positive copy of the argument, xx.
|
Chris@16
|
251 // The argument xx will be reduced to 0 <= xx <= pi/2.
|
Chris@16
|
252 bool b_negate_cos = false;
|
Chris@16
|
253
|
Chris@16
|
254 if(eval_get_sign(x) < 0)
|
Chris@16
|
255 {
|
Chris@16
|
256 xx.negate();
|
Chris@16
|
257 }
|
Chris@16
|
258
|
Chris@16
|
259 T n_pi, t;
|
Chris@16
|
260 // Remove even multiples of pi.
|
Chris@16
|
261 if(xx.compare(get_constant_pi<T>()) > 0)
|
Chris@16
|
262 {
|
Chris@16
|
263 eval_divide(t, xx, get_constant_pi<T>());
|
Chris@16
|
264 eval_trunc(n_pi, t);
|
Chris@16
|
265 BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
|
Chris@16
|
266 eval_multiply(t, n_pi, get_constant_pi<T>());
|
Chris@16
|
267 BOOST_MATH_INSTRUMENT_CODE(t.str(0, std::ios_base::scientific));
|
Chris@16
|
268 eval_subtract(xx, t);
|
Chris@16
|
269 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
Chris@16
|
270
|
Chris@16
|
271 // Adjust signs if the multiple of pi is not even.
|
Chris@16
|
272 t = ui_type(2);
|
Chris@16
|
273 eval_fmod(t, n_pi, t);
|
Chris@16
|
274 const bool b_n_pi_is_even = eval_get_sign(t) == 0;
|
Chris@16
|
275
|
Chris@16
|
276 if(!b_n_pi_is_even)
|
Chris@16
|
277 {
|
Chris@16
|
278 b_negate_cos = !b_negate_cos;
|
Chris@16
|
279 }
|
Chris@16
|
280 }
|
Chris@16
|
281
|
Chris@16
|
282 // Reduce the argument to 0 <= xx <= pi/2.
|
Chris@16
|
283 eval_ldexp(t, get_constant_pi<T>(), -1);
|
Chris@16
|
284 int com = xx.compare(t);
|
Chris@16
|
285 if(com > 0)
|
Chris@16
|
286 {
|
Chris@16
|
287 eval_subtract(xx, get_constant_pi<T>(), xx);
|
Chris@16
|
288 b_negate_cos = !b_negate_cos;
|
Chris@16
|
289 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
Chris@16
|
290 }
|
Chris@16
|
291
|
Chris@16
|
292 const bool b_zero = eval_get_sign(xx) == 0;
|
Chris@16
|
293 const bool b_pi_half = com == 0;
|
Chris@16
|
294
|
Chris@16
|
295 // Check if the reduced argument is very close to 0.
|
Chris@16
|
296 const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
|
Chris@16
|
297
|
Chris@16
|
298 if(b_zero)
|
Chris@16
|
299 {
|
Chris@16
|
300 result = si_type(1);
|
Chris@16
|
301 }
|
Chris@16
|
302 else if(b_pi_half)
|
Chris@16
|
303 {
|
Chris@16
|
304 result = si_type(0);
|
Chris@16
|
305 }
|
Chris@16
|
306 else if(b_near_zero)
|
Chris@16
|
307 {
|
Chris@16
|
308 eval_multiply(t, xx, xx);
|
Chris@16
|
309 eval_divide(t, si_type(-4));
|
Chris@16
|
310 n_pi = fp_type(0.5f);
|
Chris@16
|
311 hyp0F1(result, n_pi, t);
|
Chris@16
|
312 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
Chris@16
|
313 }
|
Chris@16
|
314 else
|
Chris@16
|
315 {
|
Chris@16
|
316 eval_subtract(t, xx);
|
Chris@16
|
317 eval_sin(result, t);
|
Chris@16
|
318 }
|
Chris@16
|
319 if(b_negate_cos)
|
Chris@16
|
320 result.negate();
|
Chris@16
|
321 }
|
Chris@16
|
322
|
Chris@16
|
323 template <class T>
|
Chris@16
|
324 void eval_tan(T& result, const T& x)
|
Chris@16
|
325 {
|
Chris@16
|
326 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tan function is only valid for floating point types.");
|
Chris@16
|
327 if(&result == &x)
|
Chris@16
|
328 {
|
Chris@16
|
329 T temp;
|
Chris@16
|
330 eval_tan(temp, x);
|
Chris@16
|
331 result = temp;
|
Chris@16
|
332 return;
|
Chris@16
|
333 }
|
Chris@16
|
334 T t;
|
Chris@16
|
335 eval_sin(result, x);
|
Chris@16
|
336 eval_cos(t, x);
|
Chris@16
|
337 eval_divide(result, t);
|
Chris@16
|
338 }
|
Chris@16
|
339
|
Chris@16
|
340 template <class T>
|
Chris@16
|
341 void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x)
|
Chris@16
|
342 {
|
Chris@16
|
343 // Compute the series representation of hyperg_2f1 taken from
|
Chris@16
|
344 // Abramowitz and Stegun 15.1.1.
|
Chris@16
|
345 // There are no checks on input range or parameter boundaries.
|
Chris@16
|
346
|
Chris@16
|
347 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
348
|
Chris@16
|
349 T x_pow_n_div_n_fact(x);
|
Chris@16
|
350 T pochham_a (a);
|
Chris@16
|
351 T pochham_b (b);
|
Chris@16
|
352 T pochham_c (c);
|
Chris@16
|
353 T ap (a);
|
Chris@16
|
354 T bp (b);
|
Chris@16
|
355 T cp (c);
|
Chris@16
|
356
|
Chris@16
|
357 eval_multiply(result, pochham_a, pochham_b);
|
Chris@16
|
358 eval_divide(result, pochham_c);
|
Chris@16
|
359 eval_multiply(result, x_pow_n_div_n_fact);
|
Chris@16
|
360 eval_add(result, ui_type(1));
|
Chris@16
|
361
|
Chris@16
|
362 T lim;
|
Chris@16
|
363 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
|
Chris@16
|
364
|
Chris@16
|
365 if(eval_get_sign(lim) < 0)
|
Chris@16
|
366 lim.negate();
|
Chris@16
|
367
|
Chris@16
|
368 ui_type n;
|
Chris@16
|
369 T term;
|
Chris@16
|
370
|
Chris@16
|
371 static const unsigned series_limit =
|
Chris@16
|
372 boost::multiprecision::detail::digits2<number<T, et_on> >::value < 100
|
Chris@16
|
373 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value;
|
Chris@16
|
374 // Series expansion of hyperg_2f1(a, b; c; x).
|
Chris@16
|
375 for(n = 2; n < series_limit; ++n)
|
Chris@16
|
376 {
|
Chris@16
|
377 eval_multiply(x_pow_n_div_n_fact, x);
|
Chris@16
|
378 eval_divide(x_pow_n_div_n_fact, n);
|
Chris@16
|
379
|
Chris@16
|
380 eval_increment(ap);
|
Chris@16
|
381 eval_multiply(pochham_a, ap);
|
Chris@16
|
382 eval_increment(bp);
|
Chris@16
|
383 eval_multiply(pochham_b, bp);
|
Chris@16
|
384 eval_increment(cp);
|
Chris@16
|
385 eval_multiply(pochham_c, cp);
|
Chris@16
|
386
|
Chris@16
|
387 eval_multiply(term, pochham_a, pochham_b);
|
Chris@16
|
388 eval_divide(term, pochham_c);
|
Chris@16
|
389 eval_multiply(term, x_pow_n_div_n_fact);
|
Chris@16
|
390 eval_add(result, term);
|
Chris@16
|
391
|
Chris@16
|
392 if(eval_get_sign(term) < 0)
|
Chris@16
|
393 term.negate();
|
Chris@16
|
394 if(lim.compare(term) >= 0)
|
Chris@16
|
395 break;
|
Chris@16
|
396 }
|
Chris@16
|
397 if(n > series_limit)
|
Chris@16
|
398 BOOST_THROW_EXCEPTION(std::runtime_error("H2F1 failed to converge."));
|
Chris@16
|
399 }
|
Chris@16
|
400
|
Chris@16
|
401 template <class T>
|
Chris@16
|
402 void eval_asin(T& result, const T& x)
|
Chris@16
|
403 {
|
Chris@16
|
404 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The asin function is only valid for floating point types.");
|
Chris@16
|
405 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
406 typedef typename mpl::front<typename T::float_types>::type fp_type;
|
Chris@16
|
407
|
Chris@16
|
408 if(&result == &x)
|
Chris@16
|
409 {
|
Chris@16
|
410 T t(x);
|
Chris@16
|
411 eval_asin(result, t);
|
Chris@16
|
412 return;
|
Chris@16
|
413 }
|
Chris@16
|
414
|
Chris@16
|
415 switch(eval_fpclassify(x))
|
Chris@16
|
416 {
|
Chris@16
|
417 case FP_NAN:
|
Chris@16
|
418 case FP_INFINITE:
|
Chris@16
|
419 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
420 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
421 else
|
Chris@16
|
422 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
423 return;
|
Chris@16
|
424 case FP_ZERO:
|
Chris@16
|
425 result = ui_type(0);
|
Chris@16
|
426 return;
|
Chris@16
|
427 default: ;
|
Chris@16
|
428 }
|
Chris@16
|
429
|
Chris@16
|
430 const bool b_neg = eval_get_sign(x) < 0;
|
Chris@16
|
431
|
Chris@16
|
432 T xx(x);
|
Chris@16
|
433 if(b_neg)
|
Chris@16
|
434 xx.negate();
|
Chris@16
|
435
|
Chris@16
|
436 int c = xx.compare(ui_type(1));
|
Chris@16
|
437 if(c > 0)
|
Chris@16
|
438 {
|
Chris@16
|
439 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
440 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
441 else
|
Chris@16
|
442 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
443 return;
|
Chris@16
|
444 }
|
Chris@16
|
445 else if(c == 0)
|
Chris@16
|
446 {
|
Chris@16
|
447 result = get_constant_pi<T>();
|
Chris@16
|
448 eval_ldexp(result, result, -1);
|
Chris@16
|
449 if(b_neg)
|
Chris@16
|
450 result.negate();
|
Chris@16
|
451 return;
|
Chris@16
|
452 }
|
Chris@16
|
453
|
Chris@16
|
454 if(xx.compare(fp_type(1e-4)) < 0)
|
Chris@16
|
455 {
|
Chris@16
|
456 // http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
|
Chris@16
|
457 eval_multiply(xx, xx);
|
Chris@16
|
458 T t1, t2;
|
Chris@16
|
459 t1 = fp_type(0.5f);
|
Chris@16
|
460 t2 = fp_type(1.5f);
|
Chris@16
|
461 hyp2F1(result, t1, t1, t2, xx);
|
Chris@16
|
462 eval_multiply(result, x);
|
Chris@16
|
463 return;
|
Chris@16
|
464 }
|
Chris@16
|
465 else if(xx.compare(fp_type(1 - 1e-4f)) > 0)
|
Chris@16
|
466 {
|
Chris@16
|
467 T dx1;
|
Chris@16
|
468 T t1, t2;
|
Chris@16
|
469 eval_subtract(dx1, ui_type(1), xx);
|
Chris@16
|
470 t1 = fp_type(0.5f);
|
Chris@16
|
471 t2 = fp_type(1.5f);
|
Chris@16
|
472 eval_ldexp(dx1, dx1, -1);
|
Chris@16
|
473 hyp2F1(result, t1, t1, t2, dx1);
|
Chris@16
|
474 eval_ldexp(dx1, dx1, 2);
|
Chris@16
|
475 eval_sqrt(t1, dx1);
|
Chris@16
|
476 eval_multiply(result, t1);
|
Chris@16
|
477 eval_ldexp(t1, get_constant_pi<T>(), -1);
|
Chris@16
|
478 result.negate();
|
Chris@16
|
479 eval_add(result, t1);
|
Chris@16
|
480 if(b_neg)
|
Chris@16
|
481 result.negate();
|
Chris@16
|
482 return;
|
Chris@16
|
483 }
|
Chris@16
|
484
|
Chris@16
|
485 // Get initial estimate using standard math function asin.
|
Chris@16
|
486 double dd;
|
Chris@16
|
487 eval_convert_to(&dd, xx);
|
Chris@16
|
488
|
Chris@16
|
489 result = fp_type(std::asin(dd));
|
Chris@16
|
490
|
Chris@101
|
491 unsigned current_digits = std::numeric_limits<double>::digits - 5;
|
Chris@101
|
492 unsigned target_precision = boost::multiprecision::detail::digits2<number<T, et_on> >::value;
|
Chris@101
|
493
|
Chris@16
|
494 // Newton-Raphson iteration
|
Chris@101
|
495 while(current_digits < target_precision)
|
Chris@16
|
496 {
|
Chris@16
|
497 T s, c;
|
Chris@16
|
498 eval_sin(s, result);
|
Chris@16
|
499 eval_cos(c, result);
|
Chris@16
|
500 eval_subtract(s, xx);
|
Chris@16
|
501 eval_divide(s, c);
|
Chris@16
|
502 eval_subtract(result, s);
|
Chris@16
|
503
|
Chris@101
|
504 current_digits *= 2;
|
Chris@101
|
505 /*
|
Chris@16
|
506 T lim;
|
Chris@16
|
507 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
|
Chris@16
|
508 if(eval_get_sign(s) < 0)
|
Chris@16
|
509 s.negate();
|
Chris@16
|
510 if(eval_get_sign(lim) < 0)
|
Chris@16
|
511 lim.negate();
|
Chris@16
|
512 if(lim.compare(s) >= 0)
|
Chris@16
|
513 break;
|
Chris@101
|
514 */
|
Chris@16
|
515 }
|
Chris@16
|
516 if(b_neg)
|
Chris@16
|
517 result.negate();
|
Chris@16
|
518 }
|
Chris@16
|
519
|
Chris@16
|
520 template <class T>
|
Chris@16
|
521 inline void eval_acos(T& result, const T& x)
|
Chris@16
|
522 {
|
Chris@16
|
523 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
|
Chris@16
|
524 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
525
|
Chris@16
|
526 switch(eval_fpclassify(x))
|
Chris@16
|
527 {
|
Chris@16
|
528 case FP_NAN:
|
Chris@16
|
529 case FP_INFINITE:
|
Chris@16
|
530 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
531 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
532 else
|
Chris@16
|
533 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
534 return;
|
Chris@16
|
535 case FP_ZERO:
|
Chris@16
|
536 result = get_constant_pi<T>();
|
Chris@16
|
537 eval_ldexp(result, result, -1); // divide by two.
|
Chris@16
|
538 return;
|
Chris@16
|
539 }
|
Chris@16
|
540
|
Chris@16
|
541 eval_abs(result, x);
|
Chris@16
|
542 int c = result.compare(ui_type(1));
|
Chris@16
|
543
|
Chris@16
|
544 if(c > 0)
|
Chris@16
|
545 {
|
Chris@16
|
546 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
547 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
548 else
|
Chris@16
|
549 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
550 return;
|
Chris@16
|
551 }
|
Chris@16
|
552 else if(c == 0)
|
Chris@16
|
553 {
|
Chris@16
|
554 if(eval_get_sign(x) < 0)
|
Chris@16
|
555 result = get_constant_pi<T>();
|
Chris@16
|
556 else
|
Chris@16
|
557 result = ui_type(0);
|
Chris@16
|
558 return;
|
Chris@16
|
559 }
|
Chris@16
|
560
|
Chris@16
|
561 eval_asin(result, x);
|
Chris@16
|
562 T t;
|
Chris@16
|
563 eval_ldexp(t, get_constant_pi<T>(), -1);
|
Chris@16
|
564 eval_subtract(result, t);
|
Chris@16
|
565 result.negate();
|
Chris@16
|
566 }
|
Chris@16
|
567
|
Chris@16
|
568 template <class T>
|
Chris@16
|
569 void eval_atan(T& result, const T& x)
|
Chris@16
|
570 {
|
Chris@16
|
571 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
|
Chris@16
|
572 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
Chris@16
|
573 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
574 typedef typename mpl::front<typename T::float_types>::type fp_type;
|
Chris@16
|
575
|
Chris@16
|
576 switch(eval_fpclassify(x))
|
Chris@16
|
577 {
|
Chris@16
|
578 case FP_NAN:
|
Chris@16
|
579 result = x;
|
Chris@16
|
580 return;
|
Chris@16
|
581 case FP_ZERO:
|
Chris@16
|
582 result = ui_type(0);
|
Chris@16
|
583 return;
|
Chris@16
|
584 case FP_INFINITE:
|
Chris@16
|
585 if(eval_get_sign(x) < 0)
|
Chris@16
|
586 {
|
Chris@16
|
587 eval_ldexp(result, get_constant_pi<T>(), -1);
|
Chris@16
|
588 result.negate();
|
Chris@16
|
589 }
|
Chris@16
|
590 else
|
Chris@16
|
591 eval_ldexp(result, get_constant_pi<T>(), -1);
|
Chris@16
|
592 return;
|
Chris@16
|
593 default: ;
|
Chris@16
|
594 }
|
Chris@16
|
595
|
Chris@16
|
596 const bool b_neg = eval_get_sign(x) < 0;
|
Chris@16
|
597
|
Chris@16
|
598 T xx(x);
|
Chris@16
|
599 if(b_neg)
|
Chris@16
|
600 xx.negate();
|
Chris@16
|
601
|
Chris@16
|
602 if(xx.compare(fp_type(0.1)) < 0)
|
Chris@16
|
603 {
|
Chris@16
|
604 T t1, t2, t3;
|
Chris@16
|
605 t1 = ui_type(1);
|
Chris@16
|
606 t2 = fp_type(0.5f);
|
Chris@16
|
607 t3 = fp_type(1.5f);
|
Chris@16
|
608 eval_multiply(xx, xx);
|
Chris@16
|
609 xx.negate();
|
Chris@16
|
610 hyp2F1(result, t1, t2, t3, xx);
|
Chris@16
|
611 eval_multiply(result, x);
|
Chris@16
|
612 return;
|
Chris@16
|
613 }
|
Chris@16
|
614
|
Chris@16
|
615 if(xx.compare(fp_type(10)) > 0)
|
Chris@16
|
616 {
|
Chris@16
|
617 T t1, t2, t3;
|
Chris@16
|
618 t1 = fp_type(0.5f);
|
Chris@16
|
619 t2 = ui_type(1u);
|
Chris@16
|
620 t3 = fp_type(1.5f);
|
Chris@16
|
621 eval_multiply(xx, xx);
|
Chris@16
|
622 eval_divide(xx, si_type(-1), xx);
|
Chris@16
|
623 hyp2F1(result, t1, t2, t3, xx);
|
Chris@16
|
624 eval_divide(result, x);
|
Chris@16
|
625 if(!b_neg)
|
Chris@16
|
626 result.negate();
|
Chris@16
|
627 eval_ldexp(t1, get_constant_pi<T>(), -1);
|
Chris@16
|
628 eval_add(result, t1);
|
Chris@16
|
629 if(b_neg)
|
Chris@16
|
630 result.negate();
|
Chris@16
|
631 return;
|
Chris@16
|
632 }
|
Chris@16
|
633
|
Chris@16
|
634
|
Chris@16
|
635 // Get initial estimate using standard math function atan.
|
Chris@16
|
636 fp_type d;
|
Chris@16
|
637 eval_convert_to(&d, xx);
|
Chris@16
|
638 result = fp_type(std::atan(d));
|
Chris@16
|
639
|
Chris@16
|
640 // Newton-Raphson iteration
|
Chris@16
|
641 static const boost::int32_t double_digits10_minus_a_few = std::numeric_limits<double>::digits10 - 3;
|
Chris@16
|
642
|
Chris@16
|
643 T s, c, t;
|
Chris@16
|
644 for(boost::int32_t digits = double_digits10_minus_a_few; digits <= std::numeric_limits<number<T, et_on> >::digits10; digits *= 2)
|
Chris@16
|
645 {
|
Chris@16
|
646 eval_sin(s, result);
|
Chris@16
|
647 eval_cos(c, result);
|
Chris@16
|
648 eval_multiply(t, xx, c);
|
Chris@16
|
649 eval_subtract(t, s);
|
Chris@16
|
650 eval_multiply(s, t, c);
|
Chris@16
|
651 eval_add(result, s);
|
Chris@16
|
652 }
|
Chris@16
|
653 if(b_neg)
|
Chris@16
|
654 result.negate();
|
Chris@16
|
655 }
|
Chris@16
|
656
|
Chris@16
|
657 template <class T>
|
Chris@16
|
658 void eval_atan2(T& result, const T& y, const T& x)
|
Chris@16
|
659 {
|
Chris@16
|
660 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
|
Chris@16
|
661 if(&result == &y)
|
Chris@16
|
662 {
|
Chris@16
|
663 T temp(y);
|
Chris@16
|
664 eval_atan2(result, temp, x);
|
Chris@16
|
665 return;
|
Chris@16
|
666 }
|
Chris@16
|
667 else if(&result == &x)
|
Chris@16
|
668 {
|
Chris@16
|
669 T temp(x);
|
Chris@16
|
670 eval_atan2(result, y, temp);
|
Chris@16
|
671 return;
|
Chris@16
|
672 }
|
Chris@16
|
673
|
Chris@16
|
674 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
Chris@16
|
675
|
Chris@16
|
676 switch(eval_fpclassify(y))
|
Chris@16
|
677 {
|
Chris@16
|
678 case FP_NAN:
|
Chris@16
|
679 result = y;
|
Chris@16
|
680 return;
|
Chris@16
|
681 case FP_ZERO:
|
Chris@16
|
682 {
|
Chris@16
|
683 int c = eval_get_sign(x);
|
Chris@16
|
684 if(c < 0)
|
Chris@16
|
685 result = get_constant_pi<T>();
|
Chris@16
|
686 else if(c >= 0)
|
Chris@16
|
687 result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
|
Chris@16
|
688 return;
|
Chris@16
|
689 }
|
Chris@16
|
690 case FP_INFINITE:
|
Chris@16
|
691 {
|
Chris@16
|
692 if(eval_fpclassify(x) == FP_INFINITE)
|
Chris@16
|
693 {
|
Chris@16
|
694 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
Chris@16
|
695 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
Chris@16
|
696 else
|
Chris@16
|
697 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
Chris@16
|
698 }
|
Chris@16
|
699 else
|
Chris@16
|
700 {
|
Chris@16
|
701 eval_ldexp(result, get_constant_pi<T>(), -1);
|
Chris@16
|
702 if(eval_get_sign(y) < 0)
|
Chris@16
|
703 result.negate();
|
Chris@16
|
704 }
|
Chris@16
|
705 return;
|
Chris@16
|
706 }
|
Chris@16
|
707 }
|
Chris@16
|
708
|
Chris@16
|
709 switch(eval_fpclassify(x))
|
Chris@16
|
710 {
|
Chris@16
|
711 case FP_NAN:
|
Chris@16
|
712 result = x;
|
Chris@16
|
713 return;
|
Chris@16
|
714 case FP_ZERO:
|
Chris@16
|
715 {
|
Chris@16
|
716 eval_ldexp(result, get_constant_pi<T>(), -1);
|
Chris@16
|
717 if(eval_get_sign(y) < 0)
|
Chris@16
|
718 result.negate();
|
Chris@16
|
719 return;
|
Chris@16
|
720 }
|
Chris@16
|
721 case FP_INFINITE:
|
Chris@16
|
722 if(eval_get_sign(x) > 0)
|
Chris@16
|
723 result = ui_type(0);
|
Chris@16
|
724 else
|
Chris@16
|
725 result = get_constant_pi<T>();
|
Chris@16
|
726 if(eval_get_sign(y) < 0)
|
Chris@16
|
727 result.negate();
|
Chris@16
|
728 return;
|
Chris@16
|
729 }
|
Chris@16
|
730
|
Chris@16
|
731 T xx;
|
Chris@16
|
732 eval_divide(xx, y, x);
|
Chris@16
|
733 if(eval_get_sign(xx) < 0)
|
Chris@16
|
734 xx.negate();
|
Chris@16
|
735
|
Chris@16
|
736 eval_atan(result, xx);
|
Chris@16
|
737
|
Chris@16
|
738 // Determine quadrant (sign) based on signs of x, y
|
Chris@16
|
739 const bool y_neg = eval_get_sign(y) < 0;
|
Chris@16
|
740 const bool x_neg = eval_get_sign(x) < 0;
|
Chris@16
|
741
|
Chris@16
|
742 if(y_neg != x_neg)
|
Chris@16
|
743 result.negate();
|
Chris@16
|
744
|
Chris@16
|
745 if(x_neg)
|
Chris@16
|
746 {
|
Chris@16
|
747 if(y_neg)
|
Chris@16
|
748 eval_subtract(result, get_constant_pi<T>());
|
Chris@16
|
749 else
|
Chris@16
|
750 eval_add(result, get_constant_pi<T>());
|
Chris@16
|
751 }
|
Chris@16
|
752 }
|
Chris@16
|
753 template<class T, class A>
|
Chris@16
|
754 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& a)
|
Chris@16
|
755 {
|
Chris@16
|
756 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
|
Chris@16
|
757 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
|
Chris@16
|
758 cast_type c;
|
Chris@16
|
759 c = a;
|
Chris@16
|
760 eval_atan2(result, x, c);
|
Chris@16
|
761 }
|
Chris@16
|
762
|
Chris@16
|
763 template<class T, class A>
|
Chris@16
|
764 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const A& x, const T& a)
|
Chris@16
|
765 {
|
Chris@16
|
766 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
|
Chris@16
|
767 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
|
Chris@16
|
768 cast_type c;
|
Chris@16
|
769 c = x;
|
Chris@16
|
770 eval_atan2(result, c, a);
|
Chris@16
|
771 }
|
Chris@16
|
772
|