boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-25 18:32:08 +00:00
Code Issues Projects Releases Wiki Activity

Math.Chebyshev: tidy up and simplify fftw dispatch code with a traits class to abstract the differences.

Browse Source
This commit is contained in:
John Maddock
2017-09-25 19:37:36 +01:00
parent 2eabd87de1
commit f7acd35e51
1 changed files with 97 additions and 200 deletions
Download Patch File Download Diff File Expand all files Collapse all files

297
include/boost/math/special_functions/chebyshev_transform.hpp
View File

@@ -17,6 +17,98 @@
namespace boost { namespace math {
namespace detail{
template <class T>
struct fftw_cos_transform;
template<>
struct fftw_cos_transform<double>
{
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<float>
{
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<long double>
{
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 Real>
class chebyshev_transform
{
@@ -25,20 +117,6 @@ public:
chebyshev_transform(const F& f, Real a, Real b,
Real tol = 500 * std::numeric_limits<Real>::epsilon(),
size_t max_refinements = 15) : m_a(a), m_b(b)
{
typedef mpl::int_ <
boost::is_same<float, Real>::value ? 0 :
boost::is_same<double, Real>::value ? 1 :
boost::is_same<long double, Real>::value ? 2 :
#ifdef BOOST_HAS_FLOAT128
boost::is_same<__float128, Real>::value ? 3 :
#endif
-1
> tag_type;
init(f, a, b, tol, max_refinements, tag_type());
}
template<class F>
void init(const F& f, Real a, Real b, Real tol, size_t max_refinements, const mpl::int_<1>&)
{
if (a >= b)
{
@@ -59,24 +137,23 @@ public:
vf.resize(n);
m_coeffs.resize(n);
fftw_plan plan = fftw_plan_r2r_1d(static_cast<int>(n), vf.data(), m_coeffs.data(), FFTW_REDFT10, FFTW_ESTIMATE);
detail::fftw_cos_transform<Real> plan(static_cast<int>(n), vf.data(), m_coeffs.data());
Real inv_n = 1/static_cast<Real>(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 = cos(pi<Real>()*(j+half<Real>())*inv_n);
Real y = detail::fftw_cos_transform<Real>::cos(pi<Real>()*(j+half<Real>())*inv_n);
vf[j] = f(y*bma + bpa)*inv_n;
vf[n-1-j]= f(bpa-y*bma)*inv_n;
}
fftw_execute_r2r(plan, vf.data(), m_coeffs.data());
fftw_destroy_plan(plan);
plan.execute(vf.data(), m_coeffs.data());
Real max_coeff = 0;
for (auto const & coeff : m_coeffs)
{
if (abs(coeff) > max_coeff)
if (detail::fftw_cos_transform<Real>::fabs(coeff) > max_coeff)
{
max_coeff = abs(coeff);
max_coeff = detail::fftw_cos_transform<Real>::fabs(coeff);
}
}
size_t j = m_coeffs.size() - 1;
@@ -97,186 +174,6 @@ public:
}
}
template<class F>
void init(const F& f, Real a, Real b, Real tol, size_t max_refinements, const mpl::int_<0>&)
{
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>();
Real bpa = (b+a)*half<Real>();
size_t n = 256;
std::vector<Real> vf;
size_t refinements = 0;
while(refinements < max_refinements)
{
vf.resize(n);
m_coeffs.resize(n);
fftwf_plan plan = fftwf_plan_r2r_1d(static_cast<int>(n), vf.data(), m_coeffs.data(), FFTW_REDFT10, FFTW_ESTIMATE);
Real inv_n = 1/static_cast<Real>(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 = cos(pi<Real>()*(j+half<Real>())*inv_n);
vf[j] = f(y*bma + bpa)*inv_n;
vf[n-1-j]= f(bpa-y*bma)*inv_n;
}
fftwf_execute_r2r(plan, vf.data(), m_coeffs.data());
fftwf_destroy_plan(plan);
Real max_coeff = 0;
for (auto const & coeff : m_coeffs)
{
if (abs(coeff) > max_coeff)
{
max_coeff = abs(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;
}
}
template<class F>
void init(const F& f, Real a, Real b, Real tol, size_t max_refinements, const mpl::int_<2>&)
{
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>();
Real bpa = (b+a)*half<Real>();
size_t n = 256;
std::vector<Real> vf;
size_t refinements = 0;
while(refinements < max_refinements)
{
vf.resize(n);
m_coeffs.resize(n);
fftwl_plan plan = fftwl_plan_r2r_1d(static_cast<int>(n), vf.data(), m_coeffs.data(), FFTW_REDFT10, FFTW_ESTIMATE);
Real inv_n = 1/static_cast<Real>(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 = cos(pi<Real>()*(j+half<Real>())*inv_n);
vf[j] = f(y*bma + bpa)*inv_n;
vf[n-1-j]= f(bpa-y*bma)*inv_n;
}
fftwl_execute_r2r(plan, vf.data(), m_coeffs.data());
fftwl_destroy_plan(plan);
Real max_coeff = 0;
for (auto const & coeff : m_coeffs)
{
if (abs(coeff) > max_coeff)
{
max_coeff = abs(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;
}
}
#ifdef BOOST_HAS_FLOAT128
template<class F>
void init(const F& f, __float128 a, __float128 b, __float128 tol, size_t max_refinements, const mpl::int_<3>&)
{
if (a >= b)
{
throw std::domain_error("a < b is required.\n");
}
using boost::math::constants::half;
using boost::math::constants::pi;
Real bma = (b-a)*half<Real>();
Real bpa = (b+a)*half<Real>();
size_t n = 256;
std::vector<Real> vf;
size_t refinements = 0;
while(refinements < max_refinements)
{
vf.resize(n);
m_coeffs.resize(n);
fftwq_plan plan = fftwq_plan_r2r_1d(static_cast<int>(n), vf.data(), m_coeffs.data(), FFTW_REDFT10, FFTW_ESTIMATE);
Real inv_n = 1/static_cast<Real>(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 = cosq(pi<Real>()*(j+half<Real>())*inv_n);
vf[j] = f(y*bma + bpa)*inv_n;
vf[n-1-j]= f(bpa-y*bma)*inv_n;
}
fftwq_execute_r2r(plan, vf.data(), m_coeffs.data());
fftwq_destroy_plan(plan);
Real max_coeff = 0;
for (auto const & coeff : m_coeffs)
{
if (fabsq(coeff) > max_coeff)
{
max_coeff = fabsq(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;
}
}
#endif
Real operator()(Real x) const
{
using boost::math::constants::half;