HiveBrain v1.2.0
Get Started
← Back to all entries
patterncppMinor

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

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
chebyshevevaluationpolynomialfoldexpressionsusingclass

Problem

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.

Two comments about the code:

-
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.

-
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
#include

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

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

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

template
struct repeat_helper> :
std::conditional_t,
// otherwise, subtract one, add a copy of T to the parameter pack,
// and continue
repeat_helper>
> {};

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

/// utility function that accepts template arguments and does nothing
/// (used with expand

Solution

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
static inline constexpr auto cpow_impl(T x, pack)
{
    // expands to x * x * x ...
    return ((noop(), 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
constexpr auto pow_impl(T x, Unsigned exponent)
    -> std::common_type_t
{
    // 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
constexpr auto pow(T x, Integer exponent)
    -> std::common_type_t
{
    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
constexpr auto pow_of_two(Integer exponent)
    -> Integer
{
    return 1 
constexpr auto pow_of_minus_one(Integer exponent)
    -> Integer
{
    return (exponent % 2) ? 1 : -1;
}

Code Snippets

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

Context

StackExchange Code Review Q#86555, answer score: 5

Revisions (0)

No revisions yet.