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