sinc.hpp 4.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145
  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/config/no_tr1/cmath.hpp>
  17. #include <boost/limits.hpp>
  18. #include <string>
  19. #include <stdexcept>
  20. #include <boost/config.hpp>
  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. T const taylor_0_bound = tools::epsilon<T>();
  34. T const taylor_2_bound = tools::root_epsilon<T>();
  35. T const taylor_n_bound = tools::forth_root_epsilon<T>();
  36. if (abs(x) >= taylor_n_bound)
  37. {
  38. return(sin(x)/x);
  39. }
  40. else
  41. {
  42. // approximation by taylor series in x at 0 up to order 0
  43. T result = static_cast<T>(1);
  44. if (abs(x) >= taylor_0_bound)
  45. {
  46. T x2 = x*x;
  47. // approximation by taylor series in x at 0 up to order 2
  48. result -= x2/static_cast<T>(6);
  49. if (abs(x) >= taylor_2_bound)
  50. {
  51. // approximation by taylor series in x at 0 up to order 4
  52. result += (x2*x2)/static_cast<T>(120);
  53. }
  54. }
  55. return(result);
  56. }
  57. }
  58. } // namespace detail
  59. template <class T>
  60. inline typename tools::promote_args<T>::type sinc_pi(T x)
  61. {
  62. typedef typename tools::promote_args<T>::type result_type;
  63. return detail::sinc_pi_imp(static_cast<result_type>(x));
  64. }
  65. template <class T, class Policy>
  66. inline typename tools::promote_args<T>::type sinc_pi(T x, const Policy&)
  67. {
  68. typedef typename tools::promote_args<T>::type result_type;
  69. return detail::sinc_pi_imp(static_cast<result_type>(x));
  70. }
  71. #ifndef BOOST_NO_TEMPLATE_TEMPLATES
  72. template<typename T, template<typename> class U>
  73. inline U<T> sinc_pi(const U<T> x)
  74. {
  75. BOOST_MATH_STD_USING
  76. using ::std::numeric_limits;
  77. T const taylor_0_bound = tools::epsilon<T>();
  78. T const taylor_2_bound = tools::root_epsilon<T>();
  79. T const taylor_n_bound = tools::forth_root_epsilon<T>();
  80. if (abs(x) >= taylor_n_bound)
  81. {
  82. return(sin(x)/x);
  83. }
  84. else
  85. {
  86. // approximation by taylor series in x at 0 up to order 0
  87. #ifdef __MWERKS__
  88. U<T> result = static_cast<U<T> >(1);
  89. #else
  90. U<T> result = U<T>(1);
  91. #endif
  92. if (abs(x) >= taylor_0_bound)
  93. {
  94. U<T> x2 = x*x;
  95. // approximation by taylor series in x at 0 up to order 2
  96. result -= x2/static_cast<T>(6);
  97. if (abs(x) >= taylor_2_bound)
  98. {
  99. // approximation by taylor series in x at 0 up to order 4
  100. result += (x2*x2)/static_cast<T>(120);
  101. }
  102. }
  103. return(result);
  104. }
  105. }
  106. template<typename T, template<typename> class U, class Policy>
  107. inline U<T> sinc_pi(const U<T> x, const Policy&)
  108. {
  109. return sinc_pi(x);
  110. }
  111. #endif /* BOOST_NO_TEMPLATE_TEMPLATES */
  112. }
  113. }
  114. #endif /* BOOST_SINC_HPP */