0

I want to write a simple polynomial class that can take an array of coefficients and expand it into a function a compile time so I don't need to loop over the coefficients at run time. I want to do something like this:

template <PARAM_TYPE, PARAMS>
class P {
public:
    PARAM_TYPE eval(PARAM_TYPE p){
        //Does PARAMS[0] * pow(p, PARAMS.length() -1) + ... + PARAMS[N-1] 
    }
}

Sample call

P<double,{2,4,3}> quadratic;
quadratic.eval(5); //returns 73

I don't want to be doing the loop since that will take time. Ideally I want to be able to form the expression above at compile time. Is this possible? Thanks

6
  • 1
    isn't that a variadic template of non-type values? Commented Mar 17, 2016 at 23:17
  • Can you elaborate? Commented Mar 17, 2016 at 23:29
  • something like this maybe? stackoverflow.com/questions/7877293/… not an exact answer, but gives the idea.. basically you get a parameter pack of integers. You can operate over them using the ... notation but you might need an en.cppreference.com/w/cpp/utility/integer_sequence as well. but I don't actually know how to do what you want... Commented Mar 17, 2016 at 23:35
  • Re "I don't want to be doing the loop since that will take time.", there will be a loop in the underlying machine code no matter what. Commented Mar 17, 2016 at 23:49
  • Do you want quadratic.eval(5) to be calculated at runtime or compile time? Commented Mar 18, 2016 at 1:05

2 Answers 2

1

Here is an example of doing what you want. The compiler is finicky about whether or not it optimizes away all the code into constants depending on the usage I noticed and which compiler you use.

test here

#include <type_traits>

template<class T, unsigned Exponent>
inline constexpr typename std::enable_if<Exponent == 0, T>::type
pow2(const T base)
{
    return 1;
}

template<class T, unsigned Exponent>
inline constexpr typename std::enable_if<Exponent % 2 != 0, T>::type
pow2(const T base)
{
      return base * pow2<T, (Exponent-1)/2>(base) * pow2<T, (Exponent-1)/2>(base);
}

template<class T, unsigned Exponent>
inline constexpr typename std::enable_if<Exponent != 0 && Exponent % 2 == 0, T>::type
pow2(const T base)
{
    return pow2<T, Exponent / 2>(base) * pow2<T, Exponent / 2>(base);
}

template<typename ParamType>
inline constexpr ParamType polynomial(const ParamType&, const ParamType& c0)
{
  return c0;
}

template<typename ParamType, typename Coeff0, typename ...Coeffs>
inline constexpr ParamType polynomial(const ParamType& x, const Coeff0& c0, const Coeffs& ...cs)
{
    return (static_cast<ParamType>(c0) * pow2<ParamType, sizeof...(cs)>(x)) + polynomial(x, static_cast<ParamType>(cs)...);
}

unsigned run(unsigned x)
{
   return polynomial(x, 2, 4, 3);
}

double run(double x)
{
  return polynomial(x, 2, 4, 3);
}

unsigned const_unsigned()
{
    static const unsigned value = polynomial(5, 2, 4, 3);
    return value;
}

double const_double()
{
    static const double value = polynomial(5, 2, 4, 3);
    return value;
}

EDIT: I have updated the code to a use a tweaked version of pow2<>() that aggressively performs calculations at compile time. This version optimizes so well at -O2 that it actually surprised me. You can see the generated assembly for the full program using the button above the code. If all arguments are constant, the compiler will generate the entire constant value at compile time. If the first argument is runtime-dependent, it generates very tight code for it still.

(Thanks to @dyp on this question for the inspiration to pow)

Sign up to request clarification or add additional context in comments.

4 Comments

So if I have polynomial(x, 0.5, 4, 3, 2, 7, 9) the compiler will generate the unrolled version of 0.5*x^5 + ... + 9?
May I ask why you needed to implement your own pow?
I have updated the code in the answer to more aggressively inline in the context of x not being a compile-time constant. This code produces surprisingly optimal code and removes all the function calls the previous version would generate when x was not known at compile time.
For your question about pow, the standard versions will convert both arguments to double and use a floating point power function when it's not necessary. Also this version allows for better inlining.
0

To evaluate a polynom, a good algorithm is Horner (See https://en.wikipedia.org/wiki/Horner%27s_method). The main idea is to compute the polynom recursively. Let's a polynom P of order n with coefficient ai. It is easy to see that the sequence Pk = Pk-1*x0 + an-k with P0 = an, that P(x0) = Pn.

So let's implement this algorithm using constexpr function:

template<class T>
constexpr double horner(double x, T an) { return an; }

template<class... T, class U = T>
constexpr double horner(double x, U an, T... a) { return horner(x, a...) * x + an; }

std::cout << horner(5., 1, 2, 1) << std::endl;

//test if the invocation of the constexpr function takes the constant expression branch
std::cout << noexcept(horner(5., 1, 2, 1)) << std::endl;

As you see, it is really easy to implement the evaluation of a polynom with constexpr functions using the recursive formula.

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.