/usr/include/boost/math/special_functions
Edit: /usr/include/boost/math/special_functions/chebyshev_transform.hpp (6209B)
// (C) Copyright Nick Thompson 2017.
// Use, modification and distribution are subject to 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)
#ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP
#define BOOST_MATH_SPECIAL_CHEBYSHEV_TRANSFORM_HPP
#include
#include
#include
#include
#include
#ifdef BOOST_HAS_FLOAT128
#include
#endif
namespace boost { namespace math {
namespace detail{
template
struct fftw_cos_transform;
template<>
struct fftw_cos_transform
{
fftw_cos_transform(int n, double* data1, double* data2)
{
plan = fftw_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
}
~fftw_cos_transform()
{
fftw_destroy_plan(plan);
}
void execute(double* data1, double* data2)
{
fftw_execute_r2r(plan, data1, data2);
}
static double cos(double x) { return std::cos(x); }
static double fabs(double x) { return std::fabs(x); }
private:
fftw_plan plan;
};
template<>
struct fftw_cos_transform
{
fftw_cos_transform(int n, float* data1, float* data2)
{
plan = fftwf_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
}
~fftw_cos_transform()
{
fftwf_destroy_plan(plan);
}
void execute(float* data1, float* data2)
{
fftwf_execute_r2r(plan, data1, data2);
}
static float cos(float x) { return std::cos(x); }
static float fabs(float x) { return std::fabs(x); }
private:
fftwf_plan plan;
};
template<>
struct fftw_cos_transform
{
fftw_cos_transform(int n, long double* data1, long double* data2)
{
plan = fftwl_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
}
~fftw_cos_transform()
{
fftwl_destroy_plan(plan);
}
void execute(long double* data1, long double* data2)
{
fftwl_execute_r2r(plan, data1, data2);
}
static long double cos(long double x) { return std::cos(x); }
static long double fabs(long double x) { return std::fabs(x); }
private:
fftwl_plan plan;
};
#ifdef BOOST_HAS_FLOAT128
template<>
struct fftw_cos_transform<__float128>
{
fftw_cos_transform(int n, __float128* data1, __float128* data2)
{
plan = fftwq_plan_r2r_1d(n, data1, data2, FFTW_REDFT10, FFTW_ESTIMATE);
}
~fftw_cos_transform()
{
fftwq_destroy_plan(plan);
}
void execute(__float128* data1, __float128* data2)
{
fftwq_execute_r2r(plan, data1, data2);
}
static __float128 cos(__float128 x) { return cosq(x); }
static __float128 fabs(__float128 x) { return fabsq(x); }
private:
fftwq_plan plan;
};
#endif
}
template
class chebyshev_transform
{
public:
template
chebyshev_transform(const F& f, Real a, Real b,
Real tol = 500 * std::numeric_limits::epsilon(),
size_t max_refinements = 15) : m_a(a), m_b(b)
{
if (a >= b)
{
throw std::domain_error("a < b is required.\n");
}
using boost::math::constants::half;
using boost::math::constants::pi;
using std::cos;
using std::abs;
Real bma = (b-a)*half();
Real bpa = (b+a)*half();
size_t n = 256;
std::vector vf;
size_t refinements = 0;
while(refinements < max_refinements)
{
vf.resize(n);
m_coeffs.resize(n);
detail::fftw_cos_transform plan(static_cast(n), vf.data(), m_coeffs.data());
Real inv_n = 1/static_cast(n);
for(size_t j = 0; j < n/2; ++j)
{
// Use symmetry cos((j+1/2)pi/n) = - cos((n-1-j+1/2)pi/n)
Real y = detail::fftw_cos_transform::cos(pi()*(j+half())*inv_n);
vf[j] = f(y*bma + bpa)*inv_n;
vf[n-1-j]= f(bpa-y*bma)*inv_n;
}
plan.execute(vf.data(), m_coeffs.data());
Real max_coeff = 0;
for (auto const & coeff : m_coeffs)
{
if (detail::fftw_cos_transform::fabs(coeff) > max_coeff)
{
max_coeff = detail::fftw_cos_transform::fabs(coeff);
}
}
size_t j = m_coeffs.size() - 1;
while (abs(m_coeffs[j])/max_coeff < tol)
{
--j;
}
// If ten coefficients are eliminated, the we say we've done all
// we need to do:
if (n - j > 10)
{
m_coeffs.resize(j+1);
return;
}
n *= 2;
++refinements;
}
}
Real operator()(Real x) const
{
using boost::math::constants::half;
if (x > m_b || x < m_a)
{
throw std::domain_error("x not in [a, b]\n");
}
Real z = (2*x - m_a - m_b)/(m_b - m_a);
return chebyshev_clenshaw_recurrence(m_coeffs.data(), m_coeffs.size(), z);
}
// Integral over entire domain [a, b]
Real integrate() const
{
Real Q = m_coeffs[0]/2;
for(size_t j = 2; j < m_coeffs.size(); j += 2)
{
Q += -m_coeffs[j]/((j+1)*(j-1));
}
return (m_b - m_a)*Q;
}
const std::vector& coefficients() const
{
return m_coeffs;
}
Real prime(Real x) const
{
Real z = (2*x - m_a - m_b)/(m_b - m_a);
Real dzdx = 2/(m_b - m_a);
if (m_coeffs.size() < 2)
{
return 0;
}
Real b2 = 0;
Real d2 = 0;
Real b1 = m_coeffs[m_coeffs.size() -1];
Real d1 = 0;
for(size_t j = m_coeffs.size() - 2; j >= 1; --j)
{
Real tmp1 = 2*z*b1 - b2 + m_coeffs[j];
Real tmp2 = 2*z*d1 - d2 + 2*b1;
b2 = b1;
b1 = tmp1;
d2 = d1;
d1 = tmp2;
}
return dzdx*(z*d1 - d2 + b1);
}
private:
std::vector m_coeffs;
Real m_a;
Real m_b;
};
}}
#endif