math/special_functions/chebyshev_transform.hpp

0001 //  (C) Copyright Nick Thompson 2017.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP
0007 #define BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP
0008 #include <cmath>
0009 #include <type_traits>
0010 #include <boost/math/constants/constants.hpp>
0011 #include <boost/math/special_functions/chebyshev.hpp>
0012
0013 #ifdef BOOST_HAS_FLOAT128
0014 #include <quadmath.h>
0015 #endif
0016
0017 #ifdef __has_include
0018 #  if __has_include(<fftw3.h>)
0019 #    include <fftw3.h>
0020 #  else
0021 #    error "This feature is unavailable without fftw3 installed"
0022 #endif
0023 #endif
0024
0025 namespace boost { namespace math {
0026
0027 namespace detail{
0028
0029 template <class T>
0030 struct fftw_cos_transform;
0031
0032 template<>
0033 struct fftw_cos_transform<double>
0034 {
0035    fftw_cos_transform(int n, double* data1, double* data2)
0036    {
0037       plan = fftw_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
0038    }
0039    ~fftw_cos_transform()
0040    {
0041       fftw_destroy_plan(plan);
0042    }
0043    void execute(double* data1, double* data2)
0044    {
0045       fftw_execute_r2r(plan, data1, data2);
0046    }
0047    static double cos(double x) { return std::cos(x); }
0048    static double fabs(double x) { return std::fabs(x); }
0049 private:
0050    fftw_plan plan;
0051 };
0052
0053 template<>
0054 struct fftw_cos_transform<float>
0055 {
0056    fftw_cos_transform(int n, float* data1, float* data2)
0057    {
0058       plan = fftwf_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
0059    }
0060    ~fftw_cos_transform()
0061    {
0062       fftwf_destroy_plan(plan);
0063    }
0064    void execute(float* data1, float* data2)
0065    {
0066       fftwf_execute_r2r(plan, data1, data2);
0067    }
0068    static float cos(float x) { return std::cos(x); }
0069    static float fabs(float x) { return std::fabs(x); }
0070 private:
0071    fftwf_plan plan;
0072 };
0073
0074 template<>
0075 struct fftw_cos_transform<long double>
0076 {
0077    fftw_cos_transform(int n, long double* data1, long double* data2)
0078    {
0079       plan = fftwl_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
0080    }
0081    ~fftw_cos_transform()
0082    {
0083       fftwl_destroy_plan(plan);
0084    }
0085    void execute(long double* data1, long double* data2)
0086    {
0087       fftwl_execute_r2r(plan, data1, data2);
0088    }
0089    static long double cos(long double x) { return std::cos(x); }
0090    static long double fabs(long double x) { return std::fabs(x); }
0091 private:
0092    fftwl_plan plan;
0093 };
0094 #ifdef BOOST_HAS_FLOAT128
0095 template<>
0096 struct fftw_cos_transform<__float128>
0097 {
0098    fftw_cos_transform(int n, __float128* data1, __float128* data2)
0099    {
0100       plan = fftwq_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
0101    }
0102    ~fftw_cos_transform()
0103    {
0104       fftwq_destroy_plan(plan);
0105    }
0106    void execute(__float128* data1, __float128* data2)
0107    {
0108       fftwq_execute_r2r(plan, data1, data2);
0109    }
0110    static __float128 cos(__float128 x) { return cosq(x); }
0111    static __float128 fabs(__float128 x) { return fabsq(x); }
0112 private:
0113    fftwq_plan plan;
0114 };
0115
0116 #endif
0117 }
0118
0119 template<class Real>
0120 class chebyshev_transform
0121 {
0122 public:
0123     template<class F>
0124     chebyshev_transform(const F& f, Real a, Real b,
0125        Real tol = 500 * std::numeric_limits<Real>::epsilon(),
0126        size_t max_refinements = 16) : m_a(a), m_b(b)
0127     {
0128         if (a >= b)
0129         {
0130             throw std::domain_error("a < b is required.\n");
0131         }
0132         using boost::math::constants::half;
0133         using boost::math::constants::pi;
0134         using std::cos;
0135         using std::abs;
0136         Real bma = (b-a)*half<Real>();
0137         Real bpa = (b+a)*half<Real>();
0138         size_t n = 256;
0139         std::vector<Real> vf;
0140
0141         size_t refinements = 0;
0142         while(refinements < max_refinements)
0143         {
0144             vf.resize(n);
0145             m_coeffs.resize(n);
0146
0147             detail::fftw_cos_transform<Real> plan(static_cast<int>(n), vf.data(), m_coeffs.data());
0148             Real inv_n = 1/static_cast<Real>(n);
0149             for(size_t j = 0; j < n/2; ++j)
0150             {
0151                 // Use symmetry cos((j+1/2)pi/n) = - cos((n-1-j+1/2)pi/n)
0152                 Real y = detail::fftw_cos_transform<Real>::cos(pi<Real>()*(j+half<Real>())*inv_n);
0153                 vf[j] = f(y*bma + bpa)*inv_n;
0154                 vf[n-1-j]= f(bpa-y*bma)*inv_n;
0155             }
0156
0157             plan.execute(vf.data(), m_coeffs.data());
0158             Real max_coeff = 0;
0159             for (auto const & coeff : m_coeffs)
0160             {
0161                 if (detail::fftw_cos_transform<Real>::fabs(coeff) > max_coeff)
0162                 {
0163                     max_coeff = detail::fftw_cos_transform<Real>::fabs(coeff);
0164                 }
0165             }
0166             size_t j = m_coeffs.size() - 1;
0167             while (abs(m_coeffs[j])/max_coeff < tol)
0168             {
0169                 --j;
0170             }
0171             // If ten coefficients are eliminated, the we say we've done all
0172             // we need to do:
0173             if (n - j > 10)
0174             {
0175                 m_coeffs.resize(j+1);
0176                 return;
0177             }
0178
0179             n *= 2;
0180             ++refinements;
0181         }
0182     }
0183
0184     inline Real operator()(Real x) const
0185     {
0186         return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), m_a, m_b, x);
0187     }
0188
0189     // Integral over entire domain [a, b]
0190     Real integrate() const
0191     {
0192           Real Q = m_coeffs[0]/2;
0193           for(size_t j = 2; j < m_coeffs.size(); j += 2)
0194           {
0195               Q += -m_coeffs[j]/((j+1)*(j-1));
0196           }
0197           return (m_b - m_a)*Q;
0198     }
0199
0200     const std::vector<Real>& coefficients() const
0201     {
0202         return m_coeffs;
0203     }
0204
0205     Real prime(Real x) const
0206     {
0207         Real z = (2*x - m_a - m_b)/(m_b - m_a);
0208         Real dzdx = 2/(m_b - m_a);
0209         if (m_coeffs.size() < 2)
0210         {
0211             return 0;
0212         }
0213         Real b2 = 0;
0214         Real d2 = 0;
0215         Real b1 = m_coeffs[m_coeffs.size() -1];
0216         Real d1 = 0;
0217         for(size_t j = m_coeffs.size() - 2; j >= 1; --j)
0218         {
0219             Real tmp1 = 2*z*b1 - b2 + m_coeffs[j];
0220             Real tmp2 = 2*z*d1 - d2 + 2*b1;
0221             b2 = b1;
0222             b1 = tmp1;
0223
0224             d2 = d1;
0225             d1 = tmp2;
0226         }
0227         return dzdx*(z*d1 - d2 + b1);
0228     }
0229
0230 private:
0231     std::vector<Real> m_coeffs;
0232     Real m_a;
0233     Real m_b;
0234 };
0235
0236 }}
0237 #endif