import quincy beta 17.1.0

[ceph.git] / ceph / src / boost / boost / math / special_functions / chebyshev.hpp

diff --git a/ceph/src/boost/boost/math/special_functions/chebyshev.hpp b/ceph/src/boost/boost/math/special_functions/chebyshev.hpp
index 0f4d9f2dac003a9df84ba0e3d6def844daea0339..8a870b9267bf05979e2fb660729c351d8ba83b7f 100644 (file)
--- a/ceph/src/boost/boost/math/special_functions/chebyshev.hpp
+++ b/ceph/src/boost/boost/math/special_functions/chebyshev.hpp
@@ -6,6 +6,7 @@
 #ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
 #define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
 #include <cmath>
+#include <boost/math/policies/error_handling.hpp>
 #include <boost/math/constants/constants.hpp>
 
 #if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))
@@ -26,13 +27,15 @@ inline typename tools::promote_args<T1, T2, T3>::type chebyshev_next(T1 const &
 
 namespace detail {
 
-template<class Real, bool second>
-inline Real chebyshev_imp(unsigned n, Real const & x)
+template<class Real, bool second, class Policy>
+inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)
 {
 #ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
     using std::acosh;
+#define BOOST_MATH_ACOSH_POLICY
 #else
    using boost::math::acosh;
+#define BOOST_MATH_ACOSH_POLICY , Policy()
 #endif
     using std::cosh;
     using std::pow;
@@ -52,17 +55,17 @@ inline Real chebyshev_imp(unsigned n, Real const & x)
     {
         if (x > 1)
         {
-            return cosh(n*acosh(x));
+            return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));
         }
         if (x < -1)
         {
             if (n & 1)
             {
-                return -cosh(n*acosh(-x));
+                return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
             }
             else
             {
-                return cosh(n*acosh(-x));
+                return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
             }
         }
         T1 = x;
@@ -90,7 +93,14 @@ chebyshev_t(unsigned n, Real const & x, const Policy&)
 {
    typedef typename tools::promote_args<Real>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, false>(n, static_cast<value_type>(x)), "boost::math::chebyshev_t<%1%>(unsigned, %1%)");
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, false>(n, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_t<%1%>(unsigned, %1%)");
 }
 
 template<class Real>
@@ -105,7 +115,14 @@ chebyshev_u(unsigned n, Real const & x, const Policy&)
 {
    typedef typename tools::promote_args<Real>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
-   return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, true>(n, static_cast<value_type>(x)), "boost::math::chebyshev_u<%1%>(unsigned, %1%)");
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, true>(n, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_u<%1%>(unsigned, %1%)");
 }
 
 template<class Real>
@@ -120,11 +137,17 @@ chebyshev_t_prime(unsigned n, Real const & x, const Policy&)
 {
    typedef typename tools::promote_args<Real>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
    if (n == 0)
    {
       return result_type(0);
    }
-   return policies::checked_narrowing_cast<result_type, Policy>(n * detail::chebyshev_imp<value_type, true>(n - 1, static_cast<value_type>(x)), "boost::math::chebyshev_t_prime<%1%>(unsigned, %1%)");
+   return policies::checked_narrowing_cast<result_type, Policy>(n * detail::chebyshev_imp<value_type, true>(n - 1, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_t_prime<%1%>(unsigned, %1%)");
 }
 
 template<class Real>
@@ -165,5 +188,100 @@ inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, co
 }
 
 
+
+namespace detail {
+template<class Real>
+inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
+{
+    Real t;
+    Real u;
+    // This cutoff is not super well defined, but it's a good estimate.
+    // See "An Error Analysis of the Modified Clenshaw Method for Evaluating Chebyshev and Fourier Series"
+    // J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379–391
+    // https://doi.org/10.1093/imamat/20.3.379
+    const Real cutoff = 0.6;
+    if (x - a < b - x)
+    {
+        u = 2*(x-a)/(b-a);
+        t = u - 1;
+        if (t > -cutoff)
+        {
+            Real b2 = 0;
+            Real b1 = c[length -1];
+            for(size_t j = length - 2; j >= 1; --j)
+            {
+                Real tmp = 2*t*b1 - b2 + c[j];
+                b2 = b1;
+                b1 = tmp;
+            }
+            return t*b1 - b2 + c[0]/2;
+        }
+        else
+        {
+            Real b = c[length -1];
+            Real d = b;
+            Real b2 = 0;
+            for (size_t r = length - 2; r >= 1; --r)
+            {
+                d = 2*u*b - d + c[r];
+                b2 = b;
+                b = d - b;
+            }
+            return t*b - b2 + c[0]/2;
+        }
+    }
+    else
+    {
+        u = -2*(b-x)/(b-a);
+        t = u + 1;
+        if (t < cutoff)
+        {
+            Real b2 = 0;
+            Real b1 = c[length -1];
+            for(size_t j = length - 2; j >= 1; --j)
+            {
+                Real tmp = 2*t*b1 - b2 + c[j];
+                b2 = b1;
+                b1 = tmp;
+            }
+            return t*b1 - b2 + c[0]/2;
+        }
+        else
+        {
+            Real b = c[length -1];
+            Real d = b;
+            Real b2 = 0;
+            for (size_t r = length - 2; r >= 1; --r)
+            {
+                d = 2*u*b + d + c[r];
+                b2 = b;
+                b = d + b;
+            }
+            return t*b - b2 + c[0]/2;
+        }
+    }
+}
+
+} // namespace detail
+
+template<class Real>
+inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
+{
+    if (x < a || x > b)
+    {
+        throw std::domain_error("x in [a, b] is required.");
+    }
+    if (length < 2)
+    {
+        if (length == 0)
+        {
+            return 0;
+        }
+        return c[0]/2;
+    }
+    return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);
+}
+
+
 }}
 #endif