blob: 84fbf0e324ace1e6602daccb674ca9adf0399e14 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
// 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
#ifdef _MSC_VER
#pragma once
#endif
#include <boost/math/tools/config.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/policies/policy.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/limits.hpp>
#include <string>
#include <stdexcept>
#include <boost/config.hpp>
// These are the the "Sinus Cardinal" functions.
namespace boost
{
namespace math
{
namespace detail
{
// This is the "Sinus Cardinal" of index Pi.
template<typename T>
inline T sinc_pi_imp(const T x)
{
BOOST_MATH_STD_USING
T const taylor_0_bound = tools::epsilon<T>();
T const taylor_2_bound = tools::root_epsilon<T>();
T const taylor_n_bound = tools::forth_root_epsilon<T>();
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);
}
}
} // namespace detail
template <class T>
inline typename tools::promote_args<T>::type sinc_pi(T x)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::sinc_pi_imp(static_cast<result_type>(x));
}
template <class T, class Policy>
inline typename tools::promote_args<T>::type sinc_pi(T x, const Policy&)
{
typedef typename tools::promote_args<T>::type result_type;
return detail::sinc_pi_imp(static_cast<result_type>(x));
}
#ifndef BOOST_NO_TEMPLATE_TEMPLATES
template<typename T, template<typename> class U>
inline U<T> sinc_pi(const U<T> x)
{
BOOST_MATH_STD_USING
using ::std::numeric_limits;
T const taylor_0_bound = tools::epsilon<T>();
T const taylor_2_bound = tools::root_epsilon<T>();
T const taylor_n_bound = tools::forth_root_epsilon<T>();
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);
}
}
template<typename T, template<typename> class U, class Policy>
inline U<T> sinc_pi(const U<T> x, const Policy&)
{
return sinc_pi(x);
}
#endif /* BOOST_NO_TEMPLATE_TEMPLATES */
}
}
#endif /* BOOST_SINC_HPP */
|
