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 static 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>
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 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>
eval_cos(T & result,const T & x)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>
eval_tan(T & result,const T & x)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>
hyp2F1(T & result,const T & a,const T & b,const T & c,const T & x)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 static 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>
eval_asin(T & result,const T & x)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
490 // Get initial estimate using standard math function asin.
491 double dd;
492 eval_convert_to(&dd, xx);
493
494 result = fp_type(std::asin(dd));
495
496 unsigned current_digits = std::numeric_limits<double>::digits - 5;
497 unsigned target_precision = boost::multiprecision::detail::digits2<number<T, et_on> >::value;
498
499 // Newton-Raphson iteration
500 while(current_digits < target_precision)
501 {
502 T sine, cosine;
503 eval_sin(sine, result);
504 eval_cos(cosine, result);
505 eval_subtract(sine, xx);
506 eval_divide(sine, cosine);
507 eval_subtract(result, sine);
508
509 current_digits *= 2;
510 /*
511 T lim;
512 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
513 if(eval_get_sign(s) < 0)
514 s.negate();
515 if(eval_get_sign(lim) < 0)
516 lim.negate();
517 if(lim.compare(s) >= 0)
518 break;
519 */
520 }
521 if(b_neg)
522 result.negate();
523 }
524
525 template <class T>
eval_acos(T & result,const T & x)526 inline void eval_acos(T& result, const T& x)
527 {
528 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
529 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
530
531 switch(eval_fpclassify(x))
532 {
533 case FP_NAN:
534 case FP_INFINITE:
535 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
536 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
537 else
538 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
539 return;
540 case FP_ZERO:
541 result = get_constant_pi<T>();
542 eval_ldexp(result, result, -1); // divide by two.
543 return;
544 }
545
546 eval_abs(result, x);
547 int c = result.compare(ui_type(1));
548
549 if(c > 0)
550 {
551 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
552 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
553 else
554 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
555 return;
556 }
557 else if(c == 0)
558 {
559 if(eval_get_sign(x) < 0)
560 result = get_constant_pi<T>();
561 else
562 result = ui_type(0);
563 return;
564 }
565
566 eval_asin(result, x);
567 T t;
568 eval_ldexp(t, get_constant_pi<T>(), -1);
569 eval_subtract(result, t);
570 result.negate();
571 }
572
573 template <class T>
eval_atan(T & result,const T & x)574 void eval_atan(T& result, const T& x)
575 {
576 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
577 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
578 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
579 typedef typename mpl::front<typename T::float_types>::type fp_type;
580
581 switch(eval_fpclassify(x))
582 {
583 case FP_NAN:
584 result = x;
585 return;
586 case FP_ZERO:
587 result = ui_type(0);
588 return;
589 case FP_INFINITE:
590 if(eval_get_sign(x) < 0)
591 {
592 eval_ldexp(result, get_constant_pi<T>(), -1);
593 result.negate();
594 }
595 else
596 eval_ldexp(result, get_constant_pi<T>(), -1);
597 return;
598 default: ;
599 }
600
601 const bool b_neg = eval_get_sign(x) < 0;
602
603 T xx(x);
604 if(b_neg)
605 xx.negate();
606
607 if(xx.compare(fp_type(0.1)) < 0)
608 {
609 T t1, t2, t3;
610 t1 = ui_type(1);
611 t2 = fp_type(0.5f);
612 t3 = fp_type(1.5f);
613 eval_multiply(xx, xx);
614 xx.negate();
615 hyp2F1(result, t1, t2, t3, xx);
616 eval_multiply(result, x);
617 return;
618 }
619
620 if(xx.compare(fp_type(10)) > 0)
621 {
622 T t1, t2, t3;
623 t1 = fp_type(0.5f);
624 t2 = ui_type(1u);
625 t3 = fp_type(1.5f);
626 eval_multiply(xx, xx);
627 eval_divide(xx, si_type(-1), xx);
628 hyp2F1(result, t1, t2, t3, xx);
629 eval_divide(result, x);
630 if(!b_neg)
631 result.negate();
632 eval_ldexp(t1, get_constant_pi<T>(), -1);
633 eval_add(result, t1);
634 if(b_neg)
635 result.negate();
636 return;
637 }
638
639
640 // Get initial estimate using standard math function atan.
641 fp_type d;
642 eval_convert_to(&d, xx);
643 result = fp_type(std::atan(d));
644
645 // Newton-Raphson iteration
646 static const boost::int32_t double_digits10_minus_a_few = std::numeric_limits<double>::digits10 - 3;
647
648 T s, c, t;
649 for(boost::int32_t digits = double_digits10_minus_a_few; digits <= std::numeric_limits<number<T, et_on> >::digits10; digits *= 2)
650 {
651 eval_sin(s, result);
652 eval_cos(c, result);
653 eval_multiply(t, xx, c);
654 eval_subtract(t, s);
655 eval_multiply(s, t, c);
656 eval_add(result, s);
657 }
658 if(b_neg)
659 result.negate();
660 }
661
662 template <class T>
eval_atan2(T & result,const T & y,const T & x)663 void eval_atan2(T& result, const T& y, const T& x)
664 {
665 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
666 if(&result == &y)
667 {
668 T temp(y);
669 eval_atan2(result, temp, x);
670 return;
671 }
672 else if(&result == &x)
673 {
674 T temp(x);
675 eval_atan2(result, y, temp);
676 return;
677 }
678
679 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
680
681 switch(eval_fpclassify(y))
682 {
683 case FP_NAN:
684 result = y;
685 return;
686 case FP_ZERO:
687 {
688 int c = eval_get_sign(x);
689 if(c < 0)
690 result = get_constant_pi<T>();
691 else if(c >= 0)
692 result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
693 return;
694 }
695 case FP_INFINITE:
696 {
697 if(eval_fpclassify(x) == FP_INFINITE)
698 {
699 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
700 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
701 else
702 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
703 }
704 else
705 {
706 eval_ldexp(result, get_constant_pi<T>(), -1);
707 if(eval_get_sign(y) < 0)
708 result.negate();
709 }
710 return;
711 }
712 }
713
714 switch(eval_fpclassify(x))
715 {
716 case FP_NAN:
717 result = x;
718 return;
719 case FP_ZERO:
720 {
721 eval_ldexp(result, get_constant_pi<T>(), -1);
722 if(eval_get_sign(y) < 0)
723 result.negate();
724 return;
725 }
726 case FP_INFINITE:
727 if(eval_get_sign(x) > 0)
728 result = ui_type(0);
729 else
730 result = get_constant_pi<T>();
731 if(eval_get_sign(y) < 0)
732 result.negate();
733 return;
734 }
735
736 T xx;
737 eval_divide(xx, y, x);
738 if(eval_get_sign(xx) < 0)
739 xx.negate();
740
741 eval_atan(result, xx);
742
743 // Determine quadrant (sign) based on signs of x, y
744 const bool y_neg = eval_get_sign(y) < 0;
745 const bool x_neg = eval_get_sign(x) < 0;
746
747 if(y_neg != x_neg)
748 result.negate();
749
750 if(x_neg)
751 {
752 if(y_neg)
753 eval_subtract(result, get_constant_pi<T>());
754 else
755 eval_add(result, get_constant_pi<T>());
756 }
757 }
758 template<class T, class A>
eval_atan2(T & result,const T & x,const A & a)759 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& a)
760 {
761 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
762 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
763 cast_type c;
764 c = a;
765 eval_atan2(result, x, c);
766 }
767
768 template<class T, class A>
eval_atan2(T & result,const A & x,const T & a)769 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const A& x, const T& a)
770 {
771 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
772 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
773 cast_type c;
774 c = x;
775 eval_atan2(result, c, a);
776 }
777
778 #ifdef BOOST_MSVC
779 #pragma warning(pop)
780 #endif
781