238 lines · plain
1// (C) Copyright Nick Thompson 2017.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP7#define BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP8#include <cmath>9#include <type_traits>10#include <boost/math/constants/constants.hpp>11#include <boost/math/special_functions/chebyshev.hpp>12 13#ifdef BOOST_HAS_FLOAT12814#include <quadmath.h>15#endif16 17#ifdef __has_include18# if __has_include(<fftw3.h>)19# include <fftw3.h>20# else21# error "This feature is unavailable without fftw3 installed"22#endif23#endif24 25namespace boost { namespace math {26 27namespace detail{28 29template <class T>30struct fftw_cos_transform;31 32template<>33struct fftw_cos_transform<double>34{35 fftw_cos_transform(int n, double* data1, double* data2)36 {37 plan = fftw_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);38 }39 ~fftw_cos_transform()40 {41 fftw_destroy_plan(plan);42 }43 void execute(double* data1, double* data2)44 {45 fftw_execute_r2r(plan, data1, data2);46 }47 static double cos(double x) { return std::cos(x); }48 static double fabs(double x) { return std::fabs(x); }49private:50 fftw_plan plan;51};52 53template<>54struct fftw_cos_transform<float>55{56 fftw_cos_transform(int n, float* data1, float* data2)57 {58 plan = fftwf_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);59 }60 ~fftw_cos_transform()61 {62 fftwf_destroy_plan(plan);63 }64 void execute(float* data1, float* data2)65 {66 fftwf_execute_r2r(plan, data1, data2);67 }68 static float cos(float x) { return std::cos(x); }69 static float fabs(float x) { return std::fabs(x); }70private:71 fftwf_plan plan;72};73 74template<>75struct fftw_cos_transform<long double>76{77 fftw_cos_transform(int n, long double* data1, long double* data2)78 {79 plan = fftwl_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);80 }81 ~fftw_cos_transform()82 {83 fftwl_destroy_plan(plan);84 }85 void execute(long double* data1, long double* data2)86 {87 fftwl_execute_r2r(plan, data1, data2);88 }89 static long double cos(long double x) { return std::cos(x); }90 static long double fabs(long double x) { return std::fabs(x); }91private:92 fftwl_plan plan;93};94#ifdef BOOST_HAS_FLOAT12895template<>96struct fftw_cos_transform<__float128>97{98 fftw_cos_transform(int n, __float128* data1, __float128* data2)99 {100 plan = fftwq_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);101 }102 ~fftw_cos_transform()103 {104 fftwq_destroy_plan(plan);105 }106 void execute(__float128* data1, __float128* data2)107 {108 fftwq_execute_r2r(plan, data1, data2);109 }110 static __float128 cos(__float128 x) { return cosq(x); }111 static __float128 fabs(__float128 x) { return fabsq(x); }112private:113 fftwq_plan plan;114};115 116#endif117}118 119template<class Real>120class chebyshev_transform121{122public:123 template<class F>124 chebyshev_transform(const F& f, Real a, Real b,125 Real tol = 500 * std::numeric_limits<Real>::epsilon(),126 size_t max_refinements = 16) : m_a(a), m_b(b)127 {128 if (a >= b)129 {130 throw std::domain_error("a < b is required.\n");131 }132 using boost::math::constants::half;133 using boost::math::constants::pi;134 using std::cos;135 using std::abs;136 Real bma = (b-a)*half<Real>();137 Real bpa = (b+a)*half<Real>();138 size_t n = 256;139 std::vector<Real> vf;140 141 size_t refinements = 0;142 while(refinements < max_refinements)143 {144 vf.resize(n);145 m_coeffs.resize(n);146 147 detail::fftw_cos_transform<Real> plan(static_cast<int>(n), vf.data(), m_coeffs.data());148 Real inv_n = 1/static_cast<Real>(n);149 for(size_t j = 0; j < n/2; ++j)150 {151 // Use symmetry cos((j+1/2)pi/n) = - cos((n-1-j+1/2)pi/n)152 Real y = detail::fftw_cos_transform<Real>::cos(pi<Real>()*(j+half<Real>())*inv_n);153 vf[j] = f(y*bma + bpa)*inv_n;154 vf[n-1-j]= f(bpa-y*bma)*inv_n;155 }156 157 plan.execute(vf.data(), m_coeffs.data());158 Real max_coeff = 0;159 for (auto const & coeff : m_coeffs)160 {161 if (detail::fftw_cos_transform<Real>::fabs(coeff) > max_coeff)162 {163 max_coeff = detail::fftw_cos_transform<Real>::fabs(coeff);164 }165 }166 size_t j = m_coeffs.size() - 1;167 while (abs(m_coeffs[j])/max_coeff < tol)168 {169 --j;170 }171 // If ten coefficients are eliminated, the we say we've done all172 // we need to do:173 if (n - j > 10)174 {175 m_coeffs.resize(j+1);176 return;177 }178 179 n *= 2;180 ++refinements;181 }182 }183 184 inline Real operator()(Real x) const185 {186 return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), m_a, m_b, x);187 }188 189 // Integral over entire domain [a, b]190 Real integrate() const191 {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() const201 {202 return m_coeffs;203 }204 205 Real prime(Real x) const206 {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 230private:231 std::vector<Real> m_coeffs;232 Real m_a;233 Real m_b;234};235 236}}237#endif238