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