]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // (C) Copyright Nick Thompson 2017. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP | |
7 | #define BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP | |
8 | #include <cmath> | |
9 | #include <type_traits> | |
10 | #include <fftw3.h> | |
11 | #include <boost/math/constants/constants.hpp> | |
12 | #include <boost/math/special_functions/chebyshev.hpp> | |
13 | ||
14 | #ifdef BOOST_HAS_FLOAT128 | |
15 | #include <quadmath.h> | |
16 | #endif | |
17 | ||
18 | namespace boost { namespace math { | |
19 | ||
20 | namespace detail{ | |
21 | ||
22 | template <class T> | |
23 | struct fftw_cos_transform; | |
24 | ||
25 | template<> | |
26 | struct fftw_cos_transform<double> | |
27 | { | |
28 | fftw_cos_transform(int n, double* data1, double* data2) | |
29 | { | |
30 | plan = fftw_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); | |
31 | } | |
32 | ~fftw_cos_transform() | |
33 | { | |
34 | fftw_destroy_plan(plan); | |
35 | } | |
36 | void execute(double* data1, double* data2) | |
37 | { | |
38 | fftw_execute_r2r(plan, data1, data2); | |
39 | } | |
40 | static double cos(double x) { return std::cos(x); } | |
41 | static double fabs(double x) { return std::fabs(x); } | |
42 | private: | |
43 | fftw_plan plan; | |
44 | }; | |
45 | ||
46 | template<> | |
47 | struct fftw_cos_transform<float> | |
48 | { | |
49 | fftw_cos_transform(int n, float* data1, float* data2) | |
50 | { | |
51 | plan = fftwf_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); | |
52 | } | |
53 | ~fftw_cos_transform() | |
54 | { | |
55 | fftwf_destroy_plan(plan); | |
56 | } | |
57 | void execute(float* data1, float* data2) | |
58 | { | |
59 | fftwf_execute_r2r(plan, data1, data2); | |
60 | } | |
61 | static float cos(float x) { return std::cos(x); } | |
62 | static float fabs(float x) { return std::fabs(x); } | |
63 | private: | |
64 | fftwf_plan plan; | |
65 | }; | |
66 | ||
67 | template<> | |
68 | struct fftw_cos_transform<long double> | |
69 | { | |
70 | fftw_cos_transform(int n, long double* data1, long double* data2) | |
71 | { | |
72 | plan = fftwl_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); | |
73 | } | |
74 | ~fftw_cos_transform() | |
75 | { | |
76 | fftwl_destroy_plan(plan); | |
77 | } | |
78 | void execute(long double* data1, long double* data2) | |
79 | { | |
80 | fftwl_execute_r2r(plan, data1, data2); | |
81 | } | |
82 | static long double cos(long double x) { return std::cos(x); } | |
83 | static long double fabs(long double x) { return std::fabs(x); } | |
84 | private: | |
85 | fftwl_plan plan; | |
86 | }; | |
87 | #ifdef BOOST_HAS_FLOAT128 | |
88 | template<> | |
89 | struct fftw_cos_transform<__float128> | |
90 | { | |
91 | fftw_cos_transform(int n, __float128* data1, __float128* data2) | |
92 | { | |
93 | plan = fftwq_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE); | |
94 | } | |
95 | ~fftw_cos_transform() | |
96 | { | |
97 | fftwq_destroy_plan(plan); | |
98 | } | |
99 | void execute(__float128* data1, __float128* data2) | |
100 | { | |
101 | fftwq_execute_r2r(plan, data1, data2); | |
102 | } | |
103 | static __float128 cos(__float128 x) { return cosq(x); } | |
104 | static __float128 fabs(__float128 x) { return fabsq(x); } | |
105 | private: | |
106 | fftwq_plan plan; | |
107 | }; | |
108 | ||
109 | #endif | |
110 | } | |
111 | ||
112 | template<class Real> | |
113 | class chebyshev_transform | |
114 | { | |
115 | public: | |
116 | template<class F> | |
117 | chebyshev_transform(const F& f, Real a, Real b, | |
118 | Real tol = 500 * std::numeric_limits<Real>::epsilon(), | |
119 | size_t max_refinements = 15) : m_a(a), m_b(b) | |
120 | { | |
121 | if (a >= b) | |
122 | { | |
123 | throw std::domain_error("a < b is required.\n"); | |
124 | } | |
125 | using boost::math::constants::half; | |
126 | using boost::math::constants::pi; | |
127 | using std::cos; | |
128 | using std::abs; | |
129 | Real bma = (b-a)*half<Real>(); | |
130 | Real bpa = (b+a)*half<Real>(); | |
131 | size_t n = 256; | |
132 | std::vector<Real> vf; | |
133 | ||
134 | size_t refinements = 0; | |
135 | while(refinements < max_refinements) | |
136 | { | |
137 | vf.resize(n); | |
138 | m_coeffs.resize(n); | |
139 | ||
140 | detail::fftw_cos_transform<Real> plan(static_cast<int>(n), vf.data(), m_coeffs.data()); | |
141 | Real inv_n = 1/static_cast<Real>(n); | |
142 | for(size_t j = 0; j < n/2; ++j) | |
143 | { | |
144 | // Use symmetry cos((j+1/2)pi/n) = - cos((n-1-j+1/2)pi/n) | |
145 | Real y = detail::fftw_cos_transform<Real>::cos(pi<Real>()*(j+half<Real>())*inv_n); | |
146 | vf[j] = f(y*bma + bpa)*inv_n; | |
147 | vf[n-1-j]= f(bpa-y*bma)*inv_n; | |
148 | } | |
149 | ||
150 | plan.execute(vf.data(), m_coeffs.data()); | |
151 | Real max_coeff = 0; | |
152 | for (auto const & coeff : m_coeffs) | |
153 | { | |
154 | if (detail::fftw_cos_transform<Real>::fabs(coeff) > max_coeff) | |
155 | { | |
156 | max_coeff = detail::fftw_cos_transform<Real>::fabs(coeff); | |
157 | } | |
158 | } | |
159 | size_t j = m_coeffs.size() - 1; | |
160 | while (abs(m_coeffs[j])/max_coeff < tol) | |
161 | { | |
162 | --j; | |
163 | } | |
164 | // If ten coefficients are eliminated, the we say we've done all | |
165 | // we need to do: | |
166 | if (n - j > 10) | |
167 | { | |
168 | m_coeffs.resize(j+1); | |
169 | return; | |
170 | } | |
171 | ||
172 | n *= 2; | |
173 | ++refinements; | |
174 | } | |
175 | } | |
176 | ||
177 | Real operator()(Real x) const | |
178 | { | |
179 | using boost::math::constants::half; | |
180 | if (x > m_b || x < m_a) | |
181 | { | |
182 | throw std::domain_error("x not in [a, b]\n"); | |
183 | } | |
184 | ||
185 | Real z = (2*x - m_a - m_b)/(m_b - m_a); | |
186 | return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), z); | |
187 | } | |
188 | ||
189 | // Integral over entire domain [a, b] | |
190 | Real integrate() const | |
191 | { | |
192 | Real Q = m_coeffs[0]/2; | |
193 | for(size_t j = 2; j < m_coeffs.size(); j += 2) | |
194 | { | |
195 | Q += -m_coeffs[j]/((j+1)*(j-1)); | |
196 | } | |
197 | return (m_b - m_a)*Q; | |
198 | } | |
199 | ||
200 | const std::vector<Real>& coefficients() const | |
201 | { | |
202 | return m_coeffs; | |
203 | } | |
204 | ||
205 | Real prime(Real x) const | |
206 | { | |
207 | Real z = (2*x - m_a - m_b)/(m_b - m_a); | |
208 | Real dzdx = 2/(m_b - m_a); | |
209 | if (m_coeffs.size() < 2) | |
210 | { | |
211 | return 0; | |
212 | } | |
213 | Real b2 = 0; | |
214 | Real d2 = 0; | |
215 | Real b1 = m_coeffs[m_coeffs.size() -1]; | |
216 | Real d1 = 0; | |
217 | for(size_t j = m_coeffs.size() - 2; j >= 1; --j) | |
218 | { | |
219 | Real tmp1 = 2*z*b1 - b2 + m_coeffs[j]; | |
220 | Real tmp2 = 2*z*d1 - d2 + 2*b1; | |
221 | b2 = b1; | |
222 | b1 = tmp1; | |
223 | ||
224 | d2 = d1; | |
225 | d1 = tmp2; | |
226 | } | |
227 | return dzdx*(z*d1 - d2 + b1); | |
228 | } | |
229 | ||
b32b8144 FG |
230 | private: |
231 | std::vector<Real> m_coeffs; | |
232 | Real m_a; | |
233 | Real m_b; | |
234 | }; | |
235 | ||
236 | }} | |
237 | #endif |