//  boost sinc.hpp header file

//  (C) Copyright Hubert Holin 2001.

//  Distributed under 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)

// See http://www.boost.org for updates, documentation, and revision history.

#ifndef BOOST_SINC_HPP

#define BOOST_SINC_HPP

#include <cmath>

#include <boost/limits.hpp>

#include <string>

#include <stdexcept>

#include <boost/config.hpp>

// These are the the "Sinus Cardinal" functions.

namespace boost

{

    namespace math

    {

#if        defined(__GNUC__) && (__GNUC__ < 3)

        // gcc 2.x ignores function scope using declarations,

        // put them in the scope of the enclosing namespace instead:

        using    ::std::abs;

        using    ::std::sqrt;

        using    ::std::sin;

        using    ::std::numeric_limits;

#endif    /* defined(__GNUC__) && (__GNUC__ < 3) */

        // This is the "Sinus Cardinal" of index Pi.

        template<typename T>

        inline T    sinc_pi(const T x)

        {

#ifdef    BOOST_NO_STDC_NAMESPACE

            using    ::abs;

            using    ::sin;

            using    ::sqrt;

#else    /* BOOST_NO_STDC_NAMESPACE */

            using    ::std::abs;

            using    ::std::sin;

            using    ::std::sqrt;

#endif    /* BOOST_NO_STDC_NAMESPACE */

            using    ::std::numeric_limits;

            static T const    taylor_0_bound = numeric_limits<T>::epsilon();

            static T const    taylor_2_bound = sqrt(taylor_0_bound);

            static T const    taylor_n_bound = sqrt(taylor_2_bound);

            if    (abs(x) >= taylor_n_bound)

            {

                return(sin(x)/x);

            }

            else

            {

                // approximation by taylor series in x at 0 up to order 0

                T    result = static_cast<T>(1);

                if    (abs(x) >= taylor_0_bound)

                {

                    T    x2 = x*x;

                    // approximation by taylor series in x at 0 up to order 2

                    result -= x2/static_cast<T>(6);

                    if    (abs(x) >= taylor_2_bound)

                    {

                        // approximation by taylor series in x at 0 up to order 4

                        result += (x2*x2)/static_cast<T>(120);

                    }

                }

                return(result);

            }

        }

#ifdef    BOOST_NO_TEMPLATE_TEMPLATES

#else    /* BOOST_NO_TEMPLATE_TEMPLATES */

        template<typename T, template<typename> class U>

        inline U<T>    sinc_pi(const U<T> x)

        {

#if defined(BOOST_FUNCTION_SCOPE_USING_DECLARATION_BREAKS_ADL) || defined(__GNUC__)

            using namespace std;

#elif    defined(BOOST_NO_STDC_NAMESPACE)

            using    ::abs;

            using    ::sin;

            using    ::sqrt;

#else    /* BOOST_NO_STDC_NAMESPACE */

            using    ::std::abs;

            using    ::std::sin;

            using    ::std::sqrt;

#endif    /* BOOST_NO_STDC_NAMESPACE */

            using    ::std::numeric_limits;

            static T const    taylor_0_bound = numeric_limits<T>::epsilon();

            static T const    taylor_2_bound = sqrt(taylor_0_bound);

            static T const    taylor_n_bound = sqrt(taylor_2_bound);

            if    (abs(x) >= taylor_n_bound)

            {

                return(sin(x)/x);

            }

            else

            {

                // approximation by taylor series in x at 0 up to order 0

#ifdef __MWERKS__

                U<T>    result = static_cast<U<T> >(1);

#else

                U<T>    result = U<T>(1);

#endif

                if    (abs(x) >= taylor_0_bound)

                {

                    U<T>    x2 = x*x;

                    // approximation by taylor series in x at 0 up to order 2

                    result -= x2/static_cast<T>(6);

                    if    (abs(x) >= taylor_2_bound)

                    {

                        // approximation by taylor series in x at 0 up to order 4

                        result += (x2*x2)/static_cast<T>(120);

                    }

                }

                return(result);

            }

        }

#endif    /* BOOST_NO_TEMPLATE_TEMPLATES */

    }

}

#endif /* BOOST_SINC_HPP */