Evaluation of Polynomials and Powers

The Problem

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

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+a²+a³+...) 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:

/*!
  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:

//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:

///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:

#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

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