15,118,323 members
Articles / Programming Languages / C++
Tip/Trick
Posted 17 Jun 2020

5.9K views
5 bookmarked

# Evaluation of Polynomials and Powers

Rate me:
How to efficiently evaluate a polynomial in C++
A short reminder of the basic tricks for polynomial evaluation

## The Problem

In a recent article on this site, the following fragment caught my eye:

C++
```for (j = 0; j <= Grado; j++)
zy[j] += Val[i][0] * pow(Val[i][1], j);
```

If we call Val[i][1] a and Val[i][0] b, we can see that the author evaluates b *(1+a+a2+a3+...) in a remarkably inefficient way.

This reminded me of a quote from "Numerical Recipes in C" where the authors, in their opinionated style, were saying:

Quote:

We assume that you know enough never to evaluate a polynomial this way:
p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x;
or (even worse!),
p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0);

Come the (computer) revolution, all persons found guilty of such criminal behavior will be summarily executed, and their programs won’t be!

(" William H. Press ... [and others]. (1992). Numerical recipes in C : the art of scientific computing. Cambridge [Cambridgeshire] ; New York :Cambridge University Press," page 173)

## The Code

We still don't know when the (computer) revolution envisaged by Dr. Press et al. will come but, in the meantime, you can save yourself by using the following template functions:

C++
```/*!
Evaluate a polynomial using Horner's scheme
\param x Evaluation point
\param coeff polynomial coefficients in order from lowest power (coeff[0]) to
highest power (coeff[N-1])
\param n size of coefficient's array.

\return Polynomial value in point x
coeff[n-1]*x^(n-1) + coeff[n-2]*x^(n-2) + ... + coeff[1]*x + coeff[0]
*/
template <typename T>
T poly (T x, const T* coeff, int n)
{
T val = coeff[n - 1];
for (int i = n - 2; i >= 0; i--)
{
val *= x;
val += coeff[i];
}
return val;
}

///  Evaluate a polynomial using Horner's scheme
template <typename T, size_t N>
T poly (T x, std::array<T, N>coeff)
{
return poly (x, coeff.data (), (int)N);
}

///  Evaluate a polynomial using Horner's scheme
template <typename T>
T poly (T x, std::vector<T>coeff)
{
return poly (x, coeff.data (), (int)coeff.size());
}```

There are three variants of the `poly` function: the first one accepts a "classical" C array while the second and the third accept a C++ array or, respectively, a vector.

Below are some examples of how to use them:

C++
```//coeffs for (X+1)^3
int cube[] { 1,3,3,1 };

int v = poly (2, cube, _countof(cube));
// v is 27
// ...

// f(x) = x^4 + 2x^3 + 3x^2 + 4x + 5
// f(2) = 57
v = poly (2, std::array<int, 5>{ 5, 4, 3, 2, 1 });

// Coefficients for sine approximation using Taylor series
// sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
std::vector<double> a(8);
double fact = -1;
for (int i = 1; i <= 7; i++)
{
if (i % 2)
{
fact *= -i;
a[i] = 1 / fact;
}
else
fact *= i;
}

double vv = poly (M_PI/2, a);
// vv is very close to 1.
```

## Integer Powers

Now that I hope I convinced you not to use the `pow` function to compute polynomials, let me show you a function that can be more efficient for computing integer powers:

C++
```///integer exponentiation function
template <typename T>
T ipow (T base, int exp)
{
T result = 1;
while (exp)
{
if (exp & 1)
result *= base;
exp >>= 1;
base *= base;
}
return result;
};```

Do not use this function for polynomial evaluation. You can use it in other cases where you have to compute some integer powers and you want to show off to friends and family.

In case they ask you if this function is more efficient than classical `pow` function, here are some numbers from my computer (no speed monster this one). The results are for 1000000 calls to `pow` and `ipow` functions with exponent 2 and 32 respectively.

```One second has 1000493 usec
Pow 2 results (usec):
ipow - integer base, integer power - 1219
ipow - double base, integer power - 1267
pow  - double base, integer power - 1010
pow  - double base, double power - 26184

One second has 1002864 usec
Pow 32 results (usec):
ipow - integer base, integer power - 9178
ipow - double base, integer power - 2614
pow  - double base, integer power - 24876
pow  - double base, double power - 24718```

In C++, the `pow` function is highly overloaded and, for small exponents, it is quite efficient. For larger exponents however, our `ipow` function wins the battle.

The code for the timing test is:

C++
```#define NMAX 1000000
volatile double dx;
volatile int ix;

void go (int exp)
{
UnitTest::Timer t;

t.Start ();
for (int i = 0; i < NMAX; i++)
ix = ipow (i, exp);
long long dt1 = t.GetTimeInUs ();
t.Start ();
for (int i = 0; i < NMAX; i++)
dx = ipow ((double)i, exp);
long long dt2 = t.GetTimeInUs ();
t.Start ();
for (int i = 0; i < NMAX; i++)
dx = pow ((double)i, exp);
long long dt3 = t.GetTimeInUs ();
t.Start ();
for (int i = 0; i < NMAX; i++)
dx = pow ((double)i, (double)exp);
long long dt4 = t.GetTimeInUs ();

// Check timer
t.Start ();
Sleep (1000);
long long dt5 = t.GetTimeInUs ();
cout << endl << "One second has " << dt5 << " usec"
<< endl;

cout << "Pow " << exp << " results (usec): " << endl <<
" ipow - integer base, integer power - " << dt1 << endl <<
" ipow - double base, integer power - " << dt2 << endl <<
" pow  - double base, integer power - " << dt3 << endl <<
" pow  - double base, double power - " << dt4 << endl;
}

go(2);
go(32);
```

## History

• 17th June, 2020: Initial version
• 19th June, 2020: Revised timing results; added timing code