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