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)
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
12 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
17 #pragma warning(disable:6326) // comparison of two constants
21 void hyp0F1(T& result, const T& b, const T& x)
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;
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.
30 T x_pow_n_div_n_fact(x);
34 eval_divide(result, x_pow_n_div_n_fact, pochham_b);
35 eval_add(result, 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)
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)
53 eval_multiply(x_pow_n_div_n_fact, x);
54 eval_divide(x_pow_n_div_n_fact, n);
56 eval_multiply(pochham_b, bp);
58 eval_divide(term, x_pow_n_div_n_fact, pochham_b);
59 eval_add(result, term);
61 bool neg_term = eval_get_sign(term) < 0;
64 if(term.compare(tol) <= 0)
69 BOOST_THROW_EXCEPTION(std::runtime_error("H0F1 Failed to Converge"));
74 void eval_sin(T& result, const T& x)
76 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sin function is only valid for floating point types.");
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;
89 switch(eval_fpclassify(x))
93 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
94 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
96 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
104 // Local copy of the argument
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;
112 if(eval_get_sign(x) < 0)
115 b_negate_sin = !b_negate_sin;
119 // Remove even multiples of pi.
120 if(xx.compare(get_constant_pi<T>()) > 0)
122 eval_divide(n_pi, xx, get_constant_pi<T>());
123 eval_trunc(n_pi, n_pi);
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);
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));
133 // Adjust signs if the multiple of pi is not even.
136 b_negate_sin = !b_negate_sin;
140 // Reduce the argument to 0 <= xx <= pi/2.
141 eval_ldexp(t, get_constant_pi<T>(), -1);
142 if(xx.compare(t) > 0)
144 eval_subtract(xx, get_constant_pi<T>(), xx);
145 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
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;
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;;
166 eval_multiply(t, xx, xx);
167 eval_divide(t, si_type(-4));
170 hyp0F1(result, t2, t);
171 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
172 eval_multiply(result, xx);
174 else if(b_near_pi_half)
177 eval_divide(t, si_type(-4));
180 hyp0F1(result, t2, t);
181 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
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).
190 static const si_type n_scale = 9;
191 static const si_type n_three_pow_scale = static_cast<si_type>(19683L);
193 eval_divide(xx, n_three_pow_scale);
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));
200 hyp0F1(result, t2, t);
201 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
202 eval_multiply(result, xx);
204 // Convert back using multiple angle identity.
205 for(boost::int32_t k = static_cast<boost::int32_t>(0); k < n_scale; k++)
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);
221 void eval_cos(T& result, const T& x)
223 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cos function is only valid for floating point types.");
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;
236 switch(eval_fpclassify(x))
240 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
241 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
243 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
251 // Local copy of the argument
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;
259 if(eval_get_sign(x) < 0)
265 // Remove even multiples of pi.
266 if(xx.compare(get_constant_pi<T>()) > 0)
268 eval_divide(t, xx, get_constant_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));
276 // Adjust signs if the multiple of pi is not even.
278 eval_fmod(t, n_pi, t);
279 const bool b_n_pi_is_even = eval_get_sign(t) == 0;
283 b_negate_cos = !b_negate_cos;
287 // Reduce the argument to 0 <= xx <= pi/2.
288 eval_ldexp(t, get_constant_pi<T>(), -1);
289 int com = xx.compare(t);
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));
297 const bool b_zero = eval_get_sign(xx) == 0;
298 const bool b_pi_half = com == 0;
300 // Check if the reduced argument is very close to 0.
301 const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
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));
321 eval_subtract(t, xx);
329 void eval_tan(T& result, const T& x)
331 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tan function is only valid for floating point types.");
342 eval_divide(result, t);
346 void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x)
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.
352 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
354 T x_pow_n_div_n_fact(x);
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));
368 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
370 if(eval_get_sign(lim) < 0)
376 const unsigned series_limit =
377 boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
378 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
379 // Series expansion of hyperg_2f1(a, b; c; x).
380 for(n = 2; n < series_limit; ++n)
382 eval_multiply(x_pow_n_div_n_fact, x);
383 eval_divide(x_pow_n_div_n_fact, n);
386 eval_multiply(pochham_a, ap);
388 eval_multiply(pochham_b, bp);
390 eval_multiply(pochham_c, cp);
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);
397 if(eval_get_sign(term) < 0)
399 if(lim.compare(term) >= 0)
403 BOOST_THROW_EXCEPTION(std::runtime_error("H2F1 failed to converge."));
407 void eval_asin(T& result, const T& x)
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;
416 eval_asin(result, t);
420 switch(eval_fpclassify(x))
424 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
425 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
427 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
435 const bool b_neg = eval_get_sign(x) < 0;
441 int c = xx.compare(ui_type(1));
444 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
445 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
447 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
452 result = get_constant_pi<T>();
453 eval_ldexp(result, result, -1);
459 if(xx.compare(fp_type(1e-4)) < 0)
461 // http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
462 eval_multiply(xx, xx);
466 hyp2F1(result, t1, t1, t2, xx);
467 eval_multiply(result, x);
470 else if(xx.compare(fp_type(1 - 1e-4f)) > 0)
474 eval_subtract(dx1, ui_type(1), xx);
477 eval_ldexp(dx1, dx1, -1);
478 hyp2F1(result, t1, t1, t2, dx1);
479 eval_ldexp(dx1, dx1, 2);
481 eval_multiply(result, t1);
482 eval_ldexp(t1, get_constant_pi<T>(), -1);
484 eval_add(result, t1);
489 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
490 typedef typename boost::multiprecision::detail::canonical<long double, T>::type guess_type;
492 typedef fp_type guess_type;
494 // Get initial estimate using standard math function asin.
496 eval_convert_to(&dd, xx);
498 result = (guess_type)(std::asin(dd));
500 // Newton-Raphson iteration, we should double our precision with each iteration,
501 // in practice this seems to not quite work in all cases... so terminate when we
502 // have at least 2/3 of the digits correct on the assumption that the correction
503 // we've just added will finish the job...
505 boost::intmax_t current_precision = eval_ilogb(result);
506 boost::intmax_t target_precision = current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3;
508 // Newton-Raphson iteration
509 while(current_precision > target_precision)
512 eval_sin(sine, result);
513 eval_cos(cosine, result);
514 eval_subtract(sine, xx);
515 eval_divide(sine, cosine);
516 eval_subtract(result, sine);
517 current_precision = eval_ilogb(sine);
519 if(current_precision == FP_ILOGB0)
528 inline void eval_acos(T& result, const T& x)
530 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
531 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
533 switch(eval_fpclassify(x))
537 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
538 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
540 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
543 result = get_constant_pi<T>();
544 eval_ldexp(result, result, -1); // divide by two.
549 int c = result.compare(ui_type(1));
553 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
554 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
556 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
561 if(eval_get_sign(x) < 0)
562 result = get_constant_pi<T>();
568 eval_asin(result, x);
570 eval_ldexp(t, get_constant_pi<T>(), -1);
571 eval_subtract(result, t);
576 void eval_atan(T& result, const T& x)
578 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
579 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
580 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
581 typedef typename mpl::front<typename T::float_types>::type fp_type;
583 switch(eval_fpclassify(x))
592 if(eval_get_sign(x) < 0)
594 eval_ldexp(result, get_constant_pi<T>(), -1);
598 eval_ldexp(result, get_constant_pi<T>(), -1);
603 const bool b_neg = eval_get_sign(x) < 0;
609 if(xx.compare(fp_type(0.1)) < 0)
615 eval_multiply(xx, xx);
617 hyp2F1(result, t1, t2, t3, xx);
618 eval_multiply(result, x);
622 if(xx.compare(fp_type(10)) > 0)
628 eval_multiply(xx, xx);
629 eval_divide(xx, si_type(-1), xx);
630 hyp2F1(result, t1, t2, t3, xx);
631 eval_divide(result, x);
634 eval_ldexp(t1, get_constant_pi<T>(), -1);
635 eval_add(result, t1);
642 // Get initial estimate using standard math function atan.
644 eval_convert_to(&d, xx);
645 result = fp_type(std::atan(d));
647 // Newton-Raphson iteration, we should double our precision with each iteration,
648 // in practice this seems to not quite work in all cases... so terminate when we
649 // have at least 2/3 of the digits correct on the assumption that the correction
650 // we've just added will finish the job...
652 boost::intmax_t current_precision = eval_ilogb(result);
653 boost::intmax_t target_precision = current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3;
656 while(current_precision > target_precision)
660 eval_multiply(t, xx, c);
662 eval_multiply(s, t, c);
664 current_precision = eval_ilogb(s);
666 if(current_precision == FP_ILOGB0)
675 void eval_atan2(T& result, const T& y, const T& x)
677 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
681 eval_atan2(result, temp, x);
684 else if(&result == &x)
687 eval_atan2(result, y, temp);
691 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
693 switch(eval_fpclassify(y))
700 int c = eval_get_sign(x);
702 result = get_constant_pi<T>();
704 result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
709 if(eval_fpclassify(x) == FP_INFINITE)
711 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
712 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
714 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
718 eval_ldexp(result, get_constant_pi<T>(), -1);
719 if(eval_get_sign(y) < 0)
726 switch(eval_fpclassify(x))
733 eval_ldexp(result, get_constant_pi<T>(), -1);
734 if(eval_get_sign(y) < 0)
739 if(eval_get_sign(x) > 0)
742 result = get_constant_pi<T>();
743 if(eval_get_sign(y) < 0)
749 eval_divide(xx, y, x);
750 if(eval_get_sign(xx) < 0)
753 eval_atan(result, xx);
755 // Determine quadrant (sign) based on signs of x, y
756 const bool y_neg = eval_get_sign(y) < 0;
757 const bool x_neg = eval_get_sign(x) < 0;
765 eval_subtract(result, get_constant_pi<T>());
767 eval_add(result, get_constant_pi<T>());
770 template<class T, class A>
771 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& a)
773 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
774 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
777 eval_atan2(result, x, c);
780 template<class T, class A>
781 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const A& x, const T& a)
783 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
784 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
787 eval_atan2(result, c, a);