[857] | 1 | // boost sinc.hpp header file
|
---|
| 2 |
|
---|
| 3 | // (C) Copyright Hubert Holin 2001.
|
---|
| 4 | // Distributed under the Boost Software License, Version 1.0. (See
|
---|
| 5 | // accompanying file LICENSE_1_0.txt or copy at
|
---|
| 6 | // http://www.boost.org/LICENSE_1_0.txt)
|
---|
| 7 |
|
---|
| 8 | // See http://www.boost.org for updates, documentation, and revision history.
|
---|
| 9 |
|
---|
| 10 | #ifndef BOOST_SINC_HPP
|
---|
| 11 | #define BOOST_SINC_HPP
|
---|
| 12 |
|
---|
| 13 |
|
---|
| 14 | #include <cmath>
|
---|
| 15 | #include <boost/limits.hpp>
|
---|
| 16 | #include <string>
|
---|
| 17 | #include <stdexcept>
|
---|
| 18 |
|
---|
| 19 |
|
---|
| 20 | #include <boost/config.hpp>
|
---|
| 21 |
|
---|
| 22 |
|
---|
| 23 | // These are the the "Sinus Cardinal" functions.
|
---|
| 24 |
|
---|
| 25 | namespace boost
|
---|
| 26 | {
|
---|
| 27 | namespace math
|
---|
| 28 | {
|
---|
| 29 | #if defined(__GNUC__) && (__GNUC__ < 3)
|
---|
| 30 | // gcc 2.x ignores function scope using declarations,
|
---|
| 31 | // put them in the scope of the enclosing namespace instead:
|
---|
| 32 |
|
---|
| 33 | using ::std::abs;
|
---|
| 34 | using ::std::sqrt;
|
---|
| 35 | using ::std::sin;
|
---|
| 36 |
|
---|
| 37 | using ::std::numeric_limits;
|
---|
| 38 | #endif /* defined(__GNUC__) && (__GNUC__ < 3) */
|
---|
| 39 |
|
---|
| 40 | // This is the "Sinus Cardinal" of index Pi.
|
---|
| 41 |
|
---|
| 42 | template<typename T>
|
---|
| 43 | inline T sinc_pi(const T x)
|
---|
| 44 | {
|
---|
| 45 | #ifdef BOOST_NO_STDC_NAMESPACE
|
---|
| 46 | using ::abs;
|
---|
| 47 | using ::sin;
|
---|
| 48 | using ::sqrt;
|
---|
| 49 | #else /* BOOST_NO_STDC_NAMESPACE */
|
---|
| 50 | using ::std::abs;
|
---|
| 51 | using ::std::sin;
|
---|
| 52 | using ::std::sqrt;
|
---|
| 53 | #endif /* BOOST_NO_STDC_NAMESPACE */
|
---|
| 54 |
|
---|
| 55 | using ::std::numeric_limits;
|
---|
| 56 |
|
---|
| 57 | static T const taylor_0_bound = numeric_limits<T>::epsilon();
|
---|
| 58 | static T const taylor_2_bound = sqrt(taylor_0_bound);
|
---|
| 59 | static T const taylor_n_bound = sqrt(taylor_2_bound);
|
---|
| 60 |
|
---|
| 61 | if (abs(x) >= taylor_n_bound)
|
---|
| 62 | {
|
---|
| 63 | return(sin(x)/x);
|
---|
| 64 | }
|
---|
| 65 | else
|
---|
| 66 | {
|
---|
| 67 | // approximation by taylor series in x at 0 up to order 0
|
---|
| 68 | T result = static_cast<T>(1);
|
---|
| 69 |
|
---|
| 70 | if (abs(x) >= taylor_0_bound)
|
---|
| 71 | {
|
---|
| 72 | T x2 = x*x;
|
---|
| 73 |
|
---|
| 74 | // approximation by taylor series in x at 0 up to order 2
|
---|
| 75 | result -= x2/static_cast<T>(6);
|
---|
| 76 |
|
---|
| 77 | if (abs(x) >= taylor_2_bound)
|
---|
| 78 | {
|
---|
| 79 | // approximation by taylor series in x at 0 up to order 4
|
---|
| 80 | result += (x2*x2)/static_cast<T>(120);
|
---|
| 81 | }
|
---|
| 82 | }
|
---|
| 83 |
|
---|
| 84 | return(result);
|
---|
| 85 | }
|
---|
| 86 | }
|
---|
| 87 |
|
---|
| 88 |
|
---|
| 89 | #ifdef BOOST_NO_TEMPLATE_TEMPLATES
|
---|
| 90 | #else /* BOOST_NO_TEMPLATE_TEMPLATES */
|
---|
| 91 | template<typename T, template<typename> class U>
|
---|
| 92 | inline U<T> sinc_pi(const U<T> x)
|
---|
| 93 | {
|
---|
| 94 | #if defined(BOOST_FUNCTION_SCOPE_USING_DECLARATION_BREAKS_ADL) || defined(__GNUC__)
|
---|
| 95 | using namespace std;
|
---|
| 96 | #elif defined(BOOST_NO_STDC_NAMESPACE)
|
---|
| 97 | using ::abs;
|
---|
| 98 | using ::sin;
|
---|
| 99 | using ::sqrt;
|
---|
| 100 | #else /* BOOST_NO_STDC_NAMESPACE */
|
---|
| 101 | using ::std::abs;
|
---|
| 102 | using ::std::sin;
|
---|
| 103 | using ::std::sqrt;
|
---|
| 104 | #endif /* BOOST_NO_STDC_NAMESPACE */
|
---|
| 105 |
|
---|
| 106 | using ::std::numeric_limits;
|
---|
| 107 |
|
---|
| 108 | static T const taylor_0_bound = numeric_limits<T>::epsilon();
|
---|
| 109 | static T const taylor_2_bound = sqrt(taylor_0_bound);
|
---|
| 110 | static T const taylor_n_bound = sqrt(taylor_2_bound);
|
---|
| 111 |
|
---|
| 112 | if (abs(x) >= taylor_n_bound)
|
---|
| 113 | {
|
---|
| 114 | return(sin(x)/x);
|
---|
| 115 | }
|
---|
| 116 | else
|
---|
| 117 | {
|
---|
| 118 | // approximation by taylor series in x at 0 up to order 0
|
---|
| 119 | #ifdef __MWERKS__
|
---|
| 120 | U<T> result = static_cast<U<T> >(1);
|
---|
| 121 | #else
|
---|
| 122 | U<T> result = U<T>(1);
|
---|
| 123 | #endif
|
---|
| 124 |
|
---|
| 125 | if (abs(x) >= taylor_0_bound)
|
---|
| 126 | {
|
---|
| 127 | U<T> x2 = x*x;
|
---|
| 128 |
|
---|
| 129 | // approximation by taylor series in x at 0 up to order 2
|
---|
| 130 | result -= x2/static_cast<T>(6);
|
---|
| 131 |
|
---|
| 132 | if (abs(x) >= taylor_2_bound)
|
---|
| 133 | {
|
---|
| 134 | // approximation by taylor series in x at 0 up to order 4
|
---|
| 135 | result += (x2*x2)/static_cast<T>(120);
|
---|
| 136 | }
|
---|
| 137 | }
|
---|
| 138 |
|
---|
| 139 | return(result);
|
---|
| 140 | }
|
---|
| 141 | }
|
---|
| 142 | #endif /* BOOST_NO_TEMPLATE_TEMPLATES */
|
---|
| 143 | }
|
---|
| 144 | }
|
---|
| 145 |
|
---|
| 146 | #endif /* BOOST_SINC_HPP */
|
---|