Chris@16: //  boost sinc.hpp header file
Chris@16: 
Chris@16: //  (C) Copyright Hubert Holin 2001.
Chris@16: //  Distributed under the Boost Software License, Version 1.0. (See
Chris@16: //  accompanying file LICENSE_1_0.txt or copy at
Chris@16: //  http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: // See http://www.boost.org for updates, documentation, and revision history.
Chris@16: 
Chris@16: #ifndef BOOST_SINC_HPP
Chris@16: #define BOOST_SINC_HPP
Chris@16: 
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/math/tools/config.hpp>
Chris@16: #include <boost/math/tools/precision.hpp>
Chris@16: #include <boost/math/policies/policy.hpp>
Chris@16: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <string>
Chris@16: #include <stdexcept>
Chris@16: 
Chris@16: 
Chris@16: #include <boost/config.hpp>
Chris@16: 
Chris@16: 
Chris@16: // These are the the "Sinus Cardinal" functions.
Chris@16: 
Chris@16: namespace boost
Chris@16: {
Chris@16:     namespace math
Chris@16:     {
Chris@16:        namespace detail
Chris@16:        {
Chris@16:         // This is the "Sinus Cardinal" of index Pi.
Chris@16: 
Chris@16:         template<typename T>
Chris@16:         inline T    sinc_pi_imp(const T x)
Chris@16:         {
Chris@16:             BOOST_MATH_STD_USING
Chris@16: 
Chris@16:             T const    taylor_0_bound = tools::epsilon<T>();
Chris@16:             T const    taylor_2_bound = tools::root_epsilon<T>();
Chris@16:             T const    taylor_n_bound = tools::forth_root_epsilon<T>();
Chris@16: 
Chris@16:             if    (abs(x) >= taylor_n_bound)
Chris@16:             {
Chris@16:                 return(sin(x)/x);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 // approximation by taylor series in x at 0 up to order 0
Chris@16:                 T    result = static_cast<T>(1);
Chris@16: 
Chris@16:                 if    (abs(x) >= taylor_0_bound)
Chris@16:                 {
Chris@16:                     T    x2 = x*x;
Chris@16: 
Chris@16:                     // approximation by taylor series in x at 0 up to order 2
Chris@16:                     result -= x2/static_cast<T>(6);
Chris@16: 
Chris@16:                     if    (abs(x) >= taylor_2_bound)
Chris@16:                     {
Chris@16:                         // approximation by taylor series in x at 0 up to order 4
Chris@16:                         result += (x2*x2)/static_cast<T>(120);
Chris@16:                     }
Chris@16:                 }
Chris@16: 
Chris@16:                 return(result);
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:        } // namespace detail
Chris@16: 
Chris@16:        template <class T>
Chris@16:        inline typename tools::promote_args<T>::type sinc_pi(T x)
Chris@16:        {
Chris@16:           typedef typename tools::promote_args<T>::type result_type;
Chris@16:           return detail::sinc_pi_imp(static_cast<result_type>(x));
Chris@16:        }
Chris@16: 
Chris@16:        template <class T, class Policy>
Chris@16:        inline typename tools::promote_args<T>::type sinc_pi(T x, const Policy&)
Chris@16:        {
Chris@16:           typedef typename tools::promote_args<T>::type result_type;
Chris@16:           return detail::sinc_pi_imp(static_cast<result_type>(x));
Chris@16:        }
Chris@16: 
Chris@16: #ifndef    BOOST_NO_TEMPLATE_TEMPLATES
Chris@16:         template<typename T, template<typename> class U>
Chris@16:         inline U<T>    sinc_pi(const U<T> x)
Chris@16:         {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using    ::std::numeric_limits;
Chris@16: 
Chris@16:             T const    taylor_0_bound = tools::epsilon<T>();
Chris@16:             T const    taylor_2_bound = tools::root_epsilon<T>();
Chris@16:             T const    taylor_n_bound = tools::forth_root_epsilon<T>();
Chris@16: 
Chris@16:             if    (abs(x) >= taylor_n_bound)
Chris@16:             {
Chris@16:                 return(sin(x)/x);
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 // approximation by taylor series in x at 0 up to order 0
Chris@16: #ifdef __MWERKS__
Chris@16:                 U<T>    result = static_cast<U<T> >(1);
Chris@16: #else
Chris@16:                 U<T>    result = U<T>(1);
Chris@16: #endif
Chris@16: 
Chris@16:                 if    (abs(x) >= taylor_0_bound)
Chris@16:                 {
Chris@16:                     U<T>    x2 = x*x;
Chris@16: 
Chris@16:                     // approximation by taylor series in x at 0 up to order 2
Chris@16:                     result -= x2/static_cast<T>(6);
Chris@16: 
Chris@16:                     if    (abs(x) >= taylor_2_bound)
Chris@16:                     {
Chris@16:                         // approximation by taylor series in x at 0 up to order 4
Chris@16:                         result += (x2*x2)/static_cast<T>(120);
Chris@16:                     }
Chris@16:                 }
Chris@16: 
Chris@16:                 return(result);
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         template<typename T, template<typename> class U, class Policy>
Chris@16:         inline U<T>    sinc_pi(const U<T> x, const Policy&)
Chris@16:         {
Chris@16:            return sinc_pi(x);
Chris@16:         }
Chris@16: #endif    /* BOOST_NO_TEMPLATE_TEMPLATES */
Chris@16:     }
Chris@16: }
Chris@16: 
Chris@16: #endif /* BOOST_SINC_HPP */
Chris@16: