sinc.hpp 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124
  1. // boost sinc.hpp header file
  2. // (C) Copyright Hubert Holin 2001.
  3. // Distributed under the Boost Software License, Version 1.0. (See
  4. // accompanying file LICENSE_1_0.txt or copy at
  5. // http://www.boost.org/LICENSE_1_0.txt)
  6. // See http://www.boost.org for updates, documentation, and revision history.
  7. #ifndef BOOST_SINC_HPP
  8. #define BOOST_SINC_HPP
  9. #ifdef _MSC_VER
  10. #pragma once
  11. #endif
  12. #include <boost/math/tools/config.hpp>
  13. #include <boost/math/tools/precision.hpp>
  14. #include <boost/math/policies/policy.hpp>
  15. #include <boost/math/special_functions/math_fwd.hpp>
  16. #include <boost/math/special_functions/fpclassify.hpp>
  17. #include <limits>
  18. #include <string>
  19. #include <stdexcept>
  20. #include <cmath>
  21. // These are the the "Sinus Cardinal" functions.
  22. namespace boost
  23. {
  24. namespace math
  25. {
  26. namespace detail
  27. {
  28. // This is the "Sinus Cardinal" of index Pi.
  29. template<typename T>
  30. inline T sinc_pi_imp(const T x)
  31. {
  32. BOOST_MATH_STD_USING
  33. if ((boost::math::isinf)(x))
  34. {
  35. return 0;
  36. }
  37. else if (abs(x) >= 3.3 * tools::forth_root_epsilon<T>())
  38. {
  39. return(sin(x)/x);
  40. }
  41. else
  42. {
  43. // |x| < (eps*120)^(1/4)
  44. return 1 - x * x / 6;
  45. }
  46. }
  47. } // namespace detail
  48. template <class T>
  49. inline typename tools::promote_args<T>::type sinc_pi(T x)
  50. {
  51. typedef typename tools::promote_args<T>::type result_type;
  52. return detail::sinc_pi_imp(static_cast<result_type>(x));
  53. }
  54. template <class T, class Policy>
  55. inline typename tools::promote_args<T>::type sinc_pi(T x, const Policy&)
  56. {
  57. typedef typename tools::promote_args<T>::type result_type;
  58. return detail::sinc_pi_imp(static_cast<result_type>(x));
  59. }
  60. template<typename T, template<typename> class U>
  61. inline U<T> sinc_pi(const U<T> x)
  62. {
  63. BOOST_MATH_STD_USING
  64. using ::std::numeric_limits;
  65. T const taylor_0_bound = tools::epsilon<T>();
  66. T const taylor_2_bound = tools::root_epsilon<T>();
  67. T const taylor_n_bound = tools::forth_root_epsilon<T>();
  68. if (abs(x) >= taylor_n_bound)
  69. {
  70. return(sin(x)/x);
  71. }
  72. else
  73. {
  74. // approximation by taylor series in x at 0 up to order 0
  75. #ifdef __MWERKS__
  76. U<T> result = static_cast<U<T> >(1);
  77. #else
  78. U<T> result = U<T>(1);
  79. #endif
  80. if (abs(x) >= taylor_0_bound)
  81. {
  82. U<T> x2 = x*x;
  83. // approximation by taylor series in x at 0 up to order 2
  84. result -= x2/static_cast<T>(6);
  85. if (abs(x) >= taylor_2_bound)
  86. {
  87. // approximation by taylor series in x at 0 up to order 4
  88. result += (x2*x2)/static_cast<T>(120);
  89. }
  90. }
  91. return(result);
  92. }
  93. }
  94. template<typename T, template<typename> class U, class Policy>
  95. inline U<T> sinc_pi(const U<T> x, const Policy&)
  96. {
  97. return sinc_pi(x);
  98. }
  99. }
  100. }
  101. #endif /* BOOST_SINC_HPP */