diff options
author | DongHun Kwak <dh0128.kwak@samsung.com> | 2019-12-05 15:11:01 +0900 |
---|---|---|
committer | DongHun Kwak <dh0128.kwak@samsung.com> | 2019-12-05 15:11:01 +0900 |
commit | 3fdc3e5ee96dca5b11d1694975a65200787eab86 (patch) | |
tree | 5c1733853892b8397d67706fa453a9bd978d2102 /boost/math/special_functions | |
parent | 88e602c57797660ebe0f9e15dbd64c1ff16dead3 (diff) | |
download | boost-3fdc3e5ee96dca5b11d1694975a65200787eab86.tar.gz boost-3fdc3e5ee96dca5b11d1694975a65200787eab86.tar.bz2 boost-3fdc3e5ee96dca5b11d1694975a65200787eab86.zip |
Imported Upstream version 1.66.0upstream/1.66.0
Diffstat (limited to 'boost/math/special_functions')
-rw-r--r-- | boost/math/special_functions/chebyshev.hpp | 169 | ||||
-rw-r--r-- | boost/math/special_functions/chebyshev_transform.hpp | 247 | ||||
-rw-r--r-- | boost/math/special_functions/detail/bernoulli_details.hpp | 58 | ||||
-rw-r--r-- | boost/math/special_functions/detail/bessel_i0.hpp | 2 | ||||
-rw-r--r-- | boost/math/special_functions/detail/bessel_i1.hpp | 7 | ||||
-rw-r--r-- | boost/math/special_functions/legendre.hpp | 2 | ||||
-rw-r--r-- | boost/math/special_functions/legendre_stieltjes.hpp | 6 | ||||
-rw-r--r-- | boost/math/special_functions/math_fwd.hpp | 37 | ||||
-rw-r--r-- | boost/math/special_functions/next.hpp | 2 |
9 files changed, 466 insertions, 64 deletions
diff --git a/boost/math/special_functions/chebyshev.hpp b/boost/math/special_functions/chebyshev.hpp new file mode 100644 index 0000000000..0f4d9f2dac --- /dev/null +++ b/boost/math/special_functions/chebyshev.hpp @@ -0,0 +1,169 @@ +// (C) Copyright Nick Thompson 2017. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_HPP +#define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP +#include <cmath> +#include <boost/math/constants/constants.hpp> + +#if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610)) +# define BOOST_MATH_CHEB_USE_STD_ACOSH +#endif + +#ifndef BOOST_MATH_CHEB_USE_STD_ACOSH +# include <boost/math/special_functions/acosh.hpp> +#endif + +namespace boost { namespace math { + +template<class T1, class T2, class T3> +inline typename tools::promote_args<T1, T2, T3>::type chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1) +{ + return 2*x*Tn - Tn_1; +} + +namespace detail { + +template<class Real, bool second> +inline Real chebyshev_imp(unsigned n, Real const & x) +{ +#ifdef BOOST_MATH_CHEB_USE_STD_ACOSH + using std::acosh; +#else + using boost::math::acosh; +#endif + using std::cosh; + using std::pow; + using std::sqrt; + Real T0 = 1; + Real T1; + if (second) + { + if (x > 1 || x < -1) + { + Real t = sqrt(x*x -1); + return static_cast<Real>((pow(x+t, (int)(n+1)) - pow(x-t, (int)(n+1)))/(2*t)); + } + T1 = 2*x; + } + else + { + if (x > 1) + { + return cosh(n*acosh(x)); + } + if (x < -1) + { + if (n & 1) + { + return -cosh(n*acosh(-x)); + } + else + { + return cosh(n*acosh(-x)); + } + } + T1 = x; + } + + if (n == 0) + { + return T0; + } + + unsigned l = 1; + while(l < n) + { + std::swap(T0, T1); + T1 = boost::math::chebyshev_next(x, T0, T1); + ++l; + } + return T1; +} +} // namespace detail + +template <class Real, class Policy> +inline typename tools::promote_args<Real>::type +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%)"); +} + +template<class Real> +inline typename tools::promote_args<Real>::type chebyshev_t(unsigned n, Real const & x) +{ + return chebyshev_t(n, x, policies::policy<>()); +} + +template <class Real, class Policy> +inline typename tools::promote_args<Real>::type +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%)"); +} + +template<class Real> +inline typename tools::promote_args<Real>::type chebyshev_u(unsigned n, Real const & x) +{ + return chebyshev_u(n, x, policies::policy<>()); +} + +template <class Real, class Policy> +inline typename tools::promote_args<Real>::type +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; + 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%)"); +} + +template<class Real> +inline typename tools::promote_args<Real>::type chebyshev_t_prime(unsigned n, Real const & x) +{ + return chebyshev_t_prime(n, x, policies::policy<>()); +} + +/* + * This is Algorithm 3.1 of + * Gil, Amparo, Javier Segura, and Nico M. Temme. + * Numerical methods for special functions. + * Society for Industrial and Applied Mathematics, 2007. + * https://www.siam.org/books/ot99/OT99SampleChapter.pdf + * However, our definition of c0 differs by a factor of 1/2, as stated in the docs. . . + */ +template<class Real, class T2> +inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x) +{ + using boost::math::constants::half; + if (length < 2) + { + if (length == 0) + { + return 0; + } + return c[0]/2; + } + Real b2 = 0; + Real b1 = c[length -1]; + for(size_t j = length - 2; j >= 1; --j) + { + Real tmp = 2*x*b1 - b2 + c[j]; + b2 = b1; + b1 = tmp; + } + return x*b1 - b2 + half<Real>()*c[0]; +} + + +}} +#endif diff --git a/boost/math/special_functions/chebyshev_transform.hpp b/boost/math/special_functions/chebyshev_transform.hpp new file mode 100644 index 0000000000..d4ce106d72 --- /dev/null +++ b/boost/math/special_functions/chebyshev_transform.hpp @@ -0,0 +1,247 @@ +// (C) Copyright Nick Thompson 2017. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP +#define BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP +#include <cmath> +#include <type_traits> +#include <fftw3.h> +#include <boost/math/constants/constants.hpp> +#include <boost/math/special_functions/chebyshev.hpp> + +#ifdef BOOST_HAS_FLOAT128 +#include <quadmath.h> +#endif + +namespace boost { namespace math { + +namespace detail{ + +template <class T> +struct fftw_cos_transform; + +template<> +struct fftw_cos_transform<double> +{ + fftw_cos_transform(int n, double* data1, double* data2) + { + plan = fftw_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); + } + ~fftw_cos_transform() + { + fftw_destroy_plan(plan); + } + void execute(double* data1, double* data2) + { + fftw_execute_r2r(plan, data1, data2); + } + static double cos(double x) { return std::cos(x); } + static double fabs(double x) { return std::fabs(x); } +private: + fftw_plan plan; +}; + +template<> +struct fftw_cos_transform<float> +{ + fftw_cos_transform(int n, float* data1, float* data2) + { + plan = fftwf_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); + } + ~fftw_cos_transform() + { + fftwf_destroy_plan(plan); + } + void execute(float* data1, float* data2) + { + fftwf_execute_r2r(plan, data1, data2); + } + static float cos(float x) { return std::cos(x); } + static float fabs(float x) { return std::fabs(x); } +private: + fftwf_plan plan; +}; + +template<> +struct fftw_cos_transform<long double> +{ + fftw_cos_transform(int n, long double* data1, long double* data2) + { + plan = fftwl_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); + } + ~fftw_cos_transform() + { + fftwl_destroy_plan(plan); + } + void execute(long double* data1, long double* data2) + { + fftwl_execute_r2r(plan, data1, data2); + } + static long double cos(long double x) { return std::cos(x); } + static long double fabs(long double x) { return std::fabs(x); } +private: + fftwl_plan plan; +}; +#ifdef BOOST_HAS_FLOAT128 +template<> +struct fftw_cos_transform<__float128> +{ + fftw_cos_transform(int n, __float128* data1, __float128* data2) + { + plan = fftwq_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); + } + ~fftw_cos_transform() + { + fftwq_destroy_plan(plan); + } + void execute(__float128* data1, __float128* data2) + { + fftwq_execute_r2r(plan, data1, data2); + } + static __float128 cos(__float128 x) { return cosq(x); } + static __float128 fabs(__float128 x) { return fabsq(x); } +private: + fftwq_plan plan; +}; + +#endif +} + +template<class Real> +class chebyshev_transform +{ +public: + template<class F> + chebyshev_transform(const F& f, Real a, Real b, + Real tol = 500 * std::numeric_limits<Real>::epsilon(), + size_t max_refinements = 15) : m_a(a), m_b(b) + { + if (a >= b) + { + throw std::domain_error("a < b is required.\n"); + } + using boost::math::constants::half; + using boost::math::constants::pi; + using std::cos; + using std::abs; + Real bma = (b-a)*half<Real>(); + Real bpa = (b+a)*half<Real>(); + size_t n = 256; + std::vector<Real> vf; + + size_t refinements = 0; + while(refinements < max_refinements) + { + vf.resize(n); + m_coeffs.resize(n); + + detail::fftw_cos_transform<Real> plan(static_cast<int>(n), vf.data(), m_coeffs.data()); + Real inv_n = 1/static_cast<Real>(n); + for(size_t j = 0; j < n/2; ++j) + { + // Use symmetry cos((j+1/2)pi/n) = - cos((n-1-j+1/2)pi/n) + Real y = detail::fftw_cos_transform<Real>::cos(pi<Real>()*(j+half<Real>())*inv_n); + vf[j] = f(y*bma + bpa)*inv_n; + vf[n-1-j]= f(bpa-y*bma)*inv_n; + } + + plan.execute(vf.data(), m_coeffs.data()); + Real max_coeff = 0; + for (auto const & coeff : m_coeffs) + { + if (detail::fftw_cos_transform<Real>::fabs(coeff) > max_coeff) + { + max_coeff = detail::fftw_cos_transform<Real>::fabs(coeff); + } + } + size_t j = m_coeffs.size() - 1; + while (abs(m_coeffs[j])/max_coeff < tol) + { + --j; + } + // If ten coefficients are eliminated, the we say we've done all + // we need to do: + if (n - j > 10) + { + m_coeffs.resize(j+1); + return; + } + + n *= 2; + ++refinements; + } + } + + Real operator()(Real x) const + { + using boost::math::constants::half; + if (x > m_b || x < m_a) + { + throw std::domain_error("x not in [a, b]\n"); + } + + Real z = (2*x - m_a - m_b)/(m_b - m_a); + return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), z); + } + + // Integral over entire domain [a, b] + Real integrate() const + { + Real Q = m_coeffs[0]/2; + for(size_t j = 2; j < m_coeffs.size(); j += 2) + { + Q += -m_coeffs[j]/((j+1)*(j-1)); + } + return (m_b - m_a)*Q; + } + + const std::vector<Real>& coefficients() const + { + return m_coeffs; + } + + Real prime(Real x) const + { + Real z = (2*x - m_a - m_b)/(m_b - m_a); + Real dzdx = 2/(m_b - m_a); + if (m_coeffs.size() < 2) + { + return 0; + } + Real b2 = 0; + Real d2 = 0; + Real b1 = m_coeffs[m_coeffs.size() -1]; + Real d1 = 0; + for(size_t j = m_coeffs.size() - 2; j >= 1; --j) + { + Real tmp1 = 2*z*b1 - b2 + m_coeffs[j]; + Real tmp2 = 2*z*d1 - d2 + 2*b1; + b2 = b1; + b1 = tmp1; + + d2 = d1; + d1 = tmp2; + } + return dzdx*(z*d1 - d2 + b1); + } + + void print_coefficients() const + { + std::cout << "{"; + for(auto const & coeff : m_coeffs) { + std::cout << coeff << ", "; + } + std::cout << "}\n"; + } + + +private: + std::vector<Real> m_coeffs; + Real m_a; + Real m_b; +}; + +}} +#endif diff --git a/boost/math/special_functions/detail/bernoulli_details.hpp b/boost/math/special_functions/detail/bernoulli_details.hpp index ce95034439..41a59e53c6 100644 --- a/boost/math/special_functions/detail/bernoulli_details.hpp +++ b/boost/math/special_functions/detail/bernoulli_details.hpp @@ -9,67 +9,11 @@ #include <boost/config.hpp> #include <boost/detail/lightweight_mutex.hpp> +#include <boost/math/tools/atomic.hpp> #include <boost/utility/enable_if.hpp> #include <boost/math/tools/toms748_solve.hpp> #include <vector> -#ifdef BOOST_HAS_THREADS - -#ifndef BOOST_NO_CXX11_HDR_ATOMIC -# include <atomic> -# define BOOST_MATH_ATOMIC_NS std -#if ATOMIC_INT_LOCK_FREE == 2 -typedef std::atomic<int> atomic_counter_type; -typedef int atomic_integer_type; -#elif ATOMIC_SHORT_LOCK_FREE == 2 -typedef std::atomic<short> atomic_counter_type; -typedef short atomic_integer_type; -#elif ATOMIC_LONG_LOCK_FREE == 2 -typedef std::atomic<long> atomic_counter_type; -typedef long atomic_integer_type; -#elif ATOMIC_LLONG_LOCK_FREE == 2 -typedef std::atomic<long long> atomic_counter_type; -typedef long long atomic_integer_type; -#else -# define BOOST_MATH_NO_ATOMIC_INT -#endif - -#else // BOOST_NO_CXX11_HDR_ATOMIC -// -// We need Boost.Atomic, but on any platform that supports auto-linking we do -// not need to link against a separate library: -// -#define BOOST_ATOMIC_NO_LIB -#include <boost/atomic.hpp> -# define BOOST_MATH_ATOMIC_NS boost - -namespace boost{ namespace math{ namespace detail{ - -// -// We need a type to use as an atomic counter: -// -#if BOOST_ATOMIC_INT_LOCK_FREE == 2 -typedef boost::atomic<int> atomic_counter_type; -typedef int atomic_integer_type; -#elif BOOST_ATOMIC_SHORT_LOCK_FREE == 2 -typedef boost::atomic<short> atomic_counter_type; -typedef short atomic_integer_type; -#elif BOOST_ATOMIC_LONG_LOCK_FREE == 2 -typedef boost::atomic<long> atomic_counter_type; -typedef long atomic_integer_type; -#elif BOOST_ATOMIC_LLONG_LOCK_FREE == 2 -typedef boost::atomic<long long> atomic_counter_type; -typedef long long atomic_integer_type; -#else -# define BOOST_MATH_NO_ATOMIC_INT -#endif - -}}} // namespaces - -#endif // BOOST_NO_CXX11_HDR_ATOMIC - -#endif // BOOST_HAS_THREADS - namespace boost{ namespace math{ namespace detail{ // // Asymptotic expansion for B2n due to diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp index 7738229454..c70f42e57a 100644 --- a/boost/math/special_functions/detail/bessel_i0.hpp +++ b/boost/math/special_functions/detail/bessel_i0.hpp @@ -17,7 +17,7 @@ // Modified Bessel function of the first kind of order zero // we use the approximating forms derived in: -// "Rational Approximations for the Modified Bessel Function of the First Kind – I0(x) for Computations with Double Precision" +// "Rational Approximations for the Modified Bessel Function of the First Kind - I0(x) for Computations with Double Precision" // by Pavel Holoborodko, // see http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision // The actual coefficients used are our own, and extend Pavel's work to precision's other than double. diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp index a6364a6cc9..3c288d72e5 100644 --- a/boost/math/special_functions/detail/bessel_i1.hpp +++ b/boost/math/special_functions/detail/bessel_i1.hpp @@ -1,6 +1,11 @@ +// Copyright (c) 2017 John Maddock +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + // Modified Bessel function of the first kind of order zero // we use the approximating forms derived in: -// "Rational Approximations for the Modified Bessel Function of the First Kind – I1(x) for Computations with Double Precision" +// "Rational Approximations for the Modified Bessel Function of the First Kind - I1(x) for Computations with Double Precision" // by Pavel Holoborodko, // see http://www.advanpix.com/2015/11/12/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i1-for-computations-with-double-precision/ // The actual coefficients used are our own, and extend Pavel's work to precision's other than double. diff --git a/boost/math/special_functions/legendre.hpp b/boost/math/special_functions/legendre.hpp index 6028b377d3..2bfc13efb9 100644 --- a/boost/math/special_functions/legendre.hpp +++ b/boost/math/special_functions/legendre.hpp @@ -192,7 +192,7 @@ std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol) T cos_nk = cos(theta_nk); T upper_bound = cos_nk; // First guess follows from: - // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93–97; + // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93-97; T inv_n_sq = 1/static_cast<T>(n*n); T sin_nk = sin(theta_nk); T x_nk_guess = (1 - inv_n_sq/static_cast<T>(8) + inv_n_sq /static_cast<T>(8*n) - (inv_n_sq*inv_n_sq/384)*(39 - 28 / (sin_nk*sin_nk) ) )*cos_nk; diff --git a/boost/math/special_functions/legendre_stieltjes.hpp b/boost/math/special_functions/legendre_stieltjes.hpp index 5d68198870..1f30ddc3b7 100644 --- a/boost/math/special_functions/legendre_stieltjes.hpp +++ b/boost/math/special_functions/legendre_stieltjes.hpp @@ -34,7 +34,7 @@ public: { throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n"); } - m_m = m; + m_m = static_cast<int>(m); std::ptrdiff_t n = m - 1; std::ptrdiff_t q; std::ptrdiff_t r; @@ -141,11 +141,11 @@ public: { if(m_m & 1) { - Em_prime += m_a[i]*detail::legendre_p_prime_imp(2*i - 1, x, policies::policy<>()); + Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 1), x, policies::policy<>()); } else { - Em_prime += m_a[i]*detail::legendre_p_prime_imp(2*i - 2, x, policies::policy<>()); + Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 2), x, policies::policy<>()); } } return Em_prime; diff --git a/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp index 9becef9f6d..03c55bd1a4 100644 --- a/boost/math/special_functions/math_fwd.hpp +++ b/boost/math/special_functions/math_fwd.hpp @@ -263,6 +263,30 @@ namespace boost typename tools::promote_args<T1, T2, T3>::type hermite_next(unsigned n, T1 x, T2 Hn, T3 Hnm1); + template<class T1, class T2, class T3> + typename tools::promote_args<T1, T2, T3>::type chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1); + + template <class Real, class Policy> + typename tools::promote_args<Real>::type + chebyshev_t(unsigned n, Real const & x, const Policy&); + template<class Real> + typename tools::promote_args<Real>::type chebyshev_t(unsigned n, Real const & x); + + template <class Real, class Policy> + typename tools::promote_args<Real>::type + chebyshev_u(unsigned n, Real const & x, const Policy&); + template<class Real> + typename tools::promote_args<Real>::type chebyshev_u(unsigned n, Real const & x); + + template <class Real, class Policy> + typename tools::promote_args<Real>::type + chebyshev_t_prime(unsigned n, Real const & x, const Policy&); + template<class Real> + typename tools::promote_args<Real>::type chebyshev_t_prime(unsigned n, Real const & x); + + template<class Real, class T2> + Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x); + template <class T1, class T2> std::complex<typename tools::promote_args<T1, T2>::type> spherical_harmonic(unsigned n, int m, T1 theta, T2 phi); @@ -1190,6 +1214,19 @@ namespace boost \ using boost::math::hermite_next;\ \ + using boost::math::chebyshev_next;\ +\ + template<class Real>\ + Real chebyshev_t(unsigned n, Real const & x){ return ::boost::math::chebyshev_t(n, x, Policy()); }\ +\ + template<class Real>\ + Real chebyshev_u(unsigned n, Real const & x){ return ::boost::math::chebyshev_u(n, x, Policy()); }\ +\ + template<class Real>\ + Real chebyshev_t_prime(unsigned n, Real const & x){ return ::boost::math::chebyshev_t_prime(n, x, Policy()); }\ +\ + using ::boost::math::chebyshev_clenshaw_recurrence;\ +\ template <class T1, class T2>\ inline std::complex<typename boost::math::tools::promote_args<T1, T2>::type> \ spherical_harmonic(unsigned n, int m, T1 theta, T2 phi){ return boost::math::spherical_harmonic(n, m, theta, phi, Policy()); }\ diff --git a/boost/math/special_functions/next.hpp b/boost/math/special_functions/next.hpp index f23ddaffe3..606b356542 100644 --- a/boost/math/special_functions/next.hpp +++ b/boost/math/special_functions/next.hpp @@ -29,7 +29,7 @@ namespace boost{ namespace math{ namespace concepts { - struct real_concept; + class real_concept; struct std_real_concept; } |