ceph/src/boost/boost/math/special_functions/chebyshev_transform.hpp

   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
 230     void print_coefficients() const
 231     {
 232         std::cout << "{";
 233         for(auto const & coeff : m_coeffs) {
 234           std::cout << coeff << ", ";
 235         }
 236         std::cout << "}\n";
 237     }
 238
 239
 240 private:
 241     std::vector<Real> m_coeffs;
 242     Real m_a;
 243     Real m_b;
 244 };
 245
 246 }}
 247 #endif