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

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);
}

• What were the memory restrictions like? Apr 19 at 15:14
• Memory restrictions were fairly tight. The target processor was a Renesas μPD78P078Y with 60K total PROM and 1K RAM. Computations on the target were actually performed using 48-bit fixed-point. Apr 19 at 15:24

# 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\$$.

• Regarding chebeval - it's not intended to be used outside the range; it just doesn't check for that. It is to operator() as std::vector::operator[]() is to std::vector::at(). Apr 19 at 16:49
• Ah, shouldn't it be made private then? Apr 19 at 16:50
• @alephzero It is indeed absolutely standard in mathematics, but if you've ever read a paper, you know that they contain an explanation for every variable name used in a formula. It is missing in this code. And while you and I might know how to find a description of Clenshaw's algorithm, not everyone does, and adding a proper reference doesn't cost you anything but might greatly help someone else. Apr 20 at 6:27
• I disagree with your suggestion to use std::function in place of function pointers. In particular, it seems incongruous with your suggestion to replace std::vector with a std::array (std::function performs a dynamic allocation). Instead I would make a second template parameter F and replace the parameter T (*func)(T) with F &&func and the line f.push_back(func(y*bma + bpa)) with f.push_back(std::invoke(func, y*bma + bpa)). This allows any appropriate invocable object to be passed to the function without any overhead. Apr 20 at 14:00
• And, speaking of, I also disagree with replacing std::vector with a std::array. What about the cases where n is not known at compile-time, or is large enough that the coefficients array would overflow the stack? If avoiding an unnecessary dynamic allocation is important, then in C++20 perhaps coefficients could be a std::span<T> with the function that constructs ChebApprox passing it as a parameter and owning the associated memory. Apr 20 at 14:31