# Chebyshev polynomial evaluation class using C++1z fold expressions

I've been messing around with templates in C++14 and C++1z recently, and since I have somewhat of an obsession with optimization (and since my field of work relates to it) I decided to try to implement an optimized class for evaluating Chebychev polynomials of the first kind.

I eventually want to work this class into a new class for generating look-up tables at compile-time, but want to make sure what I have so far is acceptable. I'm particularly new to templates so any suggestions there would be appreciated, as well as anything related to optimization, but anything at all would be helpful really.

1. I wasn't sure if I should be forwarding types around (in case someone was using a custom class for numerical computation), so I stuck for by-value for now.

2. The code currently improves significantly with -ffast-math (roughly 60% decrease in run time). I expect that it's due to how the cpow() function is set up, but I'm not sure if there's an easy way to make it so I don't have to rely on that.

A significant amount of code is taken up by helper functions/templates:

#include <cstdint>
#include <utility>

////////////////////////////////////////////////////////////////////////////////
// Utility/Helper templates and functions.                                    //
////////////////////////////////////////////////////////////////////////////////

/// simple template class to store parameter packs
template<class ...Ts>
struct pack {using type = pack<Ts...>;};

/// helper template for repeat_t, which repeats a given type T
/// N times into a pack<> type
template<std::size_t, class, class = pack<>>
struct repeat_helper;

template<std::size_t N, class T, class ...Ts>
struct repeat_helper<N, T, pack<Ts...>> :
std::conditional_t<
(N == 0),
// if N is zero, we're done
pack<Ts...>,
// otherwise, subtract one, add a copy of T to the parameter pack,
// and continue
repeat_helper<(N - 1), T, pack<T, Ts...>>
> {};

/// make a pack<> type with the type T repeated N times
template<std::size_t N, class T = void>
using repeat_t = typename repeat_helper<N, T>::type;

/// utility function that accepts template arguments and does nothing
/// (used with expanding fold expressions)
template<class...>
static inline constexpr void noop() {}

/// implementation for compile time integer powers
template<unsigned int N, class T, class ...Types>
static inline constexpr auto cpow_impl(T x, pack<Types...>)
{
// expands to x * x * x ...
return ((noop<Types>(), x) * ...);
}

/// compile time template integer power function
template<unsigned int N, class T>
static inline constexpr auto cpow(T x)
{
return cpow_impl<N, T>(x, repeat_t<N>{});
}

/// simple factorial
constexpr uintmax_t factorial(uintmax_t n)
{
return n == 0 ? 1 : n * factorial(n - 1);
}

////////////////////////////////////////////////////////////////////////////////
// Chebyshev polynomial evaluation class                                      //
////////////////////////////////////////////////////////////////////////////////

/// template class for evaluating Chebyshev polynomials of the first kind,
/// with order N. Example use: chebychev_poly<7>::eval(0.2)
template<intmax_t N, class = std::make_integer_sequence<intmax_t, N / 2 + 1>>
struct chebyshev_poly;

template<intmax_t N, intmax_t ...Ks>
struct chebyshev_poly<N, std::integer_sequence<intmax_t, Ks...>>
{
/// coefficient numerator is N * (-1)^K * (N - K - 1)! * 2^N
template<intmax_t K>
static constexpr intmax_t numerator = N * cpow<K>(-1)
* factorial(N - K - 1) * cpow<N - 2 * K>(2);

/// coefficient denominator is 2 * K! * (N - 2 * K)!
template<intmax_t K>
static constexpr intmax_t denominator = 2 * factorial(K)
* factorial(N - 2 * K);

/// the actual coefficient for a given type
template<class T, intmax_t K>
static constexpr T coeff = numerator<K> / denominator<K>;

template<class T>
static constexpr auto eval(T x)
{
// Evaluate the polynomial at a given point x by summing all terms
// using a parameter pack. Each term is of the form:
//     coeff<K> * x^(N - 2 * K)
// where K is the term index that goes from 0 to floor(N / 2) inclusive.
//
// For each term, the coefficient is obtained through the coeff<>
// variable template and the power of x via cpow(). All terms are
// summed using a fold expression.
return ((coeff<T, Ks> * cpow<(N - 2 * Ks)>(x)) + ...);
}
};


And here's some example output, complete with terrible gnuplot default color scheme:

# Fold expressions

While the fold expressions are properly used, they do not always contribute to make the code more readable; I dare say that they are not the tool you need to solve your problem (for example in the implementation of cpow, you need to resort to yet another template expansion trick to make the whole thing work). Moreover, be aware that N4358 proposes to remove the defaults for an empty pack for the operators +, *, & and |, so if a follow-up paper is accepted, then your implementation of cpow will become ill-formed when the exponent is $0$.

Anyway, if you choose to stick to fold expressions, it is safer to explicitly write the identity element of the multiplication to avoid surprises:

template<unsigned int N, class T, class ...Types>
static inline constexpr auto cpow_impl(T x, pack<Types...>)
{
// expands to x * x * x ...
return ((noop<Types>(), x) * ... * 1);
//                              ^^^^
}


# Exponentiation by squaring

Your algorithm for exponentiation may not be the most efficient in the world. While it is "free" at runtime since everything is computed at compile-time, C++ is already known to have long compilation times so using a more efficient algorithm at compile time may help reducing the compilation time. Therefore, instead of the naive "multiply $n$ times" algorithm, I would use the exponentiation by squaring algorithm instead.

I have an old C++11 constexpr implementation, so I simply pasted it below, but you could probably use the new capabilities of constexpr in C++14 to implement a better version of it. It also handles negative exponents:

template<typename T, typename Unsigned>
constexpr auto pow_impl(T x, Unsigned exponent)
-> std::common_type_t<T, Unsigned>
{
// Exponentiation by squaring
return (exponent == 0) ? 1 :
(exponent % 2 == 0) ? pow_impl(x*x, exponent/2) :
x * pow_impl(x*x, (exponent-1)/2);
}

template<typename T, typename Integer>
constexpr auto pow(T x, Integer exponent)
-> std::common_type_t<T, Integer>
{
return (exponent == 0) ? 1 :
(exponent > 0) ? pow_impl(x, exponent) :
1 / pow_impl(x, -exponent);
}


# Specialized algorithms

Calling a generic algorithm with a constant value in a formula can sometimes be considered to be a special algorithm. In your case, $2^n$ and $-1^n$ may be considered special. Their respective implementations can easily be made $O(1)$, which is always something you might want at some point:

template<typename Integer>
constexpr auto pow_of_two(Integer exponent)
-> Integer
{
return 1 << exponent;
}

template<typename Integer>
constexpr auto pow_of_minus_one(Integer exponent)
-> Integer
{
return (exponent % 2) ? 1 : -1;
}


A better, (faster, more accurate) way to evaluate polynomials, given the coefficients, is to use horner's rule, eg

evaluate c + b*x + a*x*x via c + x*(b + x*a)


This lends itself very well to being implemented as a fold.

In the particular case of chebychev polynomials (and other orthogonal families) it is even better to use their recurrence relation. Not only does this avoid the (sometimes unpleasant) calculation of the coefficients, but will also compute the first n polynomials in a single order n loop.

For chebychev polynomials T the recurrence is

T[0](x) = 1
T[1](x) = x
T[n+1](x) = 2*x*T[n](x) - T[n-1](x)


You didn't say what you wanted the polynomials for; a common case is to evaluate things like:

Sum{ 0<=j<n | cs[j]*T[j](x) }


This can be evaluated efficiently using Clenshaw's recurrence, like this:

double  cheb_eval( int n, const double* cs, double x)
{
double  r = 0.0;
double  s = 0.0;
double  t;
double  x2 = 2.0*x;
while( --n>=1)
{   t = r;
r = x2*r - s + cs[n];
s = t;
}
return x*r - s + cs[0];
}

• Accepted the other answer since it pointed out a bug, but would have marked yours otherwise since I totally didn't think of the recurrence relations (and they are the most relevant to my primary concern of runtime performance). – cartographer May 18 '15 at 18:31
• @cartographer Don't hesitate to accept this one instead if it helped you the most; algorithm improvements are always more interesting than bugfixes and such :) I don't even know whether N4358 got accepted or not (even though I wrote it) so it's not a bug, it's a potential future bug, which is vague enough. – Morwenn May 19 '15 at 7:08