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