C++ class to create and evaluate Chebyshev approximations of arbitrary functions

Question

I recently needed to create a function to approximate a complex trigonometric function on an embedded system without a floating point unit and without a fast trigonometric library. So I pulled out my ancient copy of "Numerical Recipes in C, Second Edition" by Press, Teukolsky, Vetterling and Flannery, and updated the Chebyshev polynomial approximation functions to use modern C++.

For an approximation of order \$N\$ within the range \$[-1,1]\$, the code calculates the Chebyshev approximation coefficients of the passed function \$f\$ as: $$ c_j = \frac{2}{N}\sum_{k=1}^{N} f\left[ \cos\left(\frac{\pi(k - \frac{1}{2})}{N}\right)\right] \cos\left(\frac{\pi j (k - \frac{1}{2})}{N}\right) $$

To make the code more general, we map any arbitrary input range \$[a, b]\$ to \$[-1, 1]\$ using $$ x' = \frac{x - \frac{1}{2}(b + a)}{\frac{1}{2}(b - a)} $$

This code then uses Clenshaw's algorithm to evaluate the approximation.

The demo function takes a function to be approximated, low and high limits and the order of the approximation as inputs and calculates the coefficients and prints 100 values and the error value of the approximation at each of those points, and the maximum of the absolute error. It then exercises the checked evaluation of the approximation (which does range checking) to provoke a thrown std::range_error. The checked evaluation uses an unchecked evaluation to actually perform the Clenshaw algorithm and evaluate the approximation.

A plot of the error values for a 7th order approximation of the cosine function in the range \$[0, \frac{\pi}{2}]\$ is shown below. It has the characteristic undulating error value of a Chebyshev approximation.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <limits>
#include <numbers>
#include <stdexcept>
#include <vector>

template <typename T>
class ChebApprox {
public:
    // given a function implementation, create a Chebyshev polynomial approximation
    ChebApprox(T lowerbound, T upperbound, unsigned n, T (*func)(T));
    // range-check evaluation of Chebyshev approximation
    T operator()(T x) const;
    // unchecked evaluation of Chebyshev approximation
    T chebeval(T x) const;
    // returns beginning iterator of coefficients list
    auto begin() const { return coefficients.cbegin(); }
    // returns ending iterator of coefficients list
    auto end() const { return coefficients.cend(); }
private:
    T lowerbound; 
    T upperbound; 
    std::vector<T> coefficients;
};

template <typename T>
ChebApprox<T>::ChebApprox(T lowerbound, T upperbound, unsigned n, T (*func)(T)) 
    : lowerbound{lowerbound}
    , upperbound{upperbound}
{
    const auto bma{(upperbound - lowerbound)/2};
    const auto bpa{(upperbound + lowerbound)/2};
    std::vector<T> f;
    f.reserve(n);
    coefficients.reserve(n);
    for (unsigned k = 0; k < n; ++k) {
        T y{std::cos(std::numbers::pi * (k + 0.5) / n)};
        f.push_back((*func)(y * bma + bpa));
    }
    const T fac = 2.0/n;
    for (unsigned j=0; j < n; ++j) {
        T sum{0};
        for (unsigned k=0; k < n; ++k) {
            sum += f[k] * std::cos(std::numbers::pi * j * (k+0.5) / n);
        }
        coefficients.push_back(fac*sum);
    }
}

template <typename T>
T ChebApprox<T>::operator()(T x) const {
    if (x < lowerbound || x > upperbound) {
        throw std::range_error("Approximation function input out of allowed range");
    }
    return chebeval(x);
}

template <typename T>
T ChebApprox<T>::chebeval(T x) const {
    T y{(2 * x - lowerbound - upperbound) / (upperbound - lowerbound)};
    T y2{2 * y};
    T d{0};
    T dd{0};
    for (auto j{coefficients.size()-1}; j; --j) {
        auto sv{d};
        d = y2 * d - dd + coefficients[j];
        dd = sv;
    }
    return y * d - dd + coefficients[0] / 2;
}

void demo(double (*func)(double), double lo, double hi, unsigned n) 
{
    ChebApprox<double> ch{lo, hi, n, func};
    std::cout << "Coefficients:\n" << std::setprecision(std::numeric_limits<double>::digits10 + 1);
    std::copy(ch.begin(), ch.end(), std::ostream_iterator<double>(std::cout, ", "));
    std::cout << "\nEvaluation\n";
    auto delta{(hi-lo)/100};
    double maxerror{0};
    for (auto x{lo}; x < hi; x += delta) {
        auto calc{ch(x)};
        auto error{calc - func(x)};
        if (std::fabs(error) > std::fabs(maxerror)) { maxerror = error; }
        std::cout << x << '\t' << calc << '\t' << error << '\n';
    }
    std::cout << "Maximum absolute error was " << maxerror << '\n';
    try {
        auto bad{ch(hi + 1)};  // outside defined range
    } catch (std::range_error &err) {
        std::cout << err.what() << '\n';
    }
}

int main() {
    demo(std::cos, 0, std::numbers::pi/2, 7);
    demo(std::sqrt, 0, 100, 5);
}

Answer 1

Naming things

You have a bounds checking evaluation function, operator(), and one that doesn't do bounds checking named chebeval(). If the idea is to have this mimic an STL containers operator[] vs. at(), then perhaps make it more like that: have operator() do the unchecked evaluation, and have a member function named at() that does the bounds checks. Alternatively, keep the safe operator() as it is, but perhaps rename chebeval() to something else that makes it clear that this is not doing bounds checks (eval_unchecked()?).

You use a lot of one-letter variables, making the code hard to read. Consider making them more descriptive. Alternatively, add a comment with a link to a page/paper with the formulas and algorithms you are using, and stay as close to the notation as used in that page/paper as possible.

Consider replacing begin()/end() with coefficients()

What is the begin and end of a Chebyshev approximation? Not everyone would assume that begin() and end() would give you access to the coefficients. So I would make it more explicit, and just have a member function returning a const reference:

    auto &coefficients() const { return m_coefficients; }
private:
    std::vector<T> m_coefficients;

And then you can use the ranges library to write for example:

std::ranges::copy(ch.coefficients(), std::ostream_iterator<double>(std::cout, ", "));

Consider replacing std::vector with a std::array

You know the number of coefficients up front at construction time, and it will never change. And since class ChebApprox is already a template, you can just add a template parameter for the number of coefficients instead, and then use std::array<T, n> coefficients for a small performance improvement.

Allow function objects

You are using a raw function pointer for func, but this precludes you from passing function objects and lambdas to the constructor. In particular, you cannot create a new ChebApprox using a previously created ChebApprox as the input function. One way to do this is to use std::function, like so:

template <typename T>
class ChebApprox {
public:
    ChebApprox(T lowerbound, T upperbound, unsigned n, std::function<T(T)> func);

However, as pointed out by Ray Hamel, std::function might allocate memory, although a decent implementation will not do so if it only needs to store a raw function pointer. But we can avoid this issue by not using std::function, but instead adding a template parameter for the function type:

template <typename T, typename Func>
class ChebApprox {
    using 
public:
    ChebApprox(T lowerbound, T upperbound, unsigned n, Func &&func);

The drawback of the latter method is that Func is now no longer constrained, and passing it an overloaded function such as std::sqrt will be problematic, and would require you to explicitly cast it to something unambiguous first.

Consider using std::midpoint()

When calculating bma and bpa, use std::midpoint(). While unlikely to give problems, it is nonetheless safer to use than what you wrote, especially if the bounds are close to being denormal or infinity.

It might also be possible to use std::lerp() instead:

for (unsigned k = 0; k < n; ++k) {
    T y{1 + 0.5 * std::cos(std::numbers::pi * (k + 0.5) / n)};
    f.push_back(func(std::lerp(lowerbound, upperbound, y));
}

Avoid an expensive division in chebeval()

Divisions are one of the most expensive CPU instructions. If you have a small number of coefficients, the calculation of y might take a significant fraction of the CPU time in chebeval(). Instead of storing lowerbound and upperbound as member variables, consider storing lowerbound + upperbound and 1 / (upperbound - lowerbound) instead, so the calculation of y gets sped up.

The range check can be done by first calculating y, then checking if it is in the range \$-1\dots 1\$.