1

I need to evaluate a polynomial (1 dimension) of some degree (<10) in multiple points. For instance, when calculating for polynomial p(x) = 4x^3 + 2*x^2 -2*x + 5, at x = 0, x= 2 result should be : [5, 41]. Couldn't find an equivalent function to Matlab's Polyval.

The result accuracy is not critical for me (can be rounded to integer) but the calculation time is.

My straight-forward solution:

#include<iostream>
#include<vector>
#include<math.h>

std::vector<double> Polyval(std::vector<double> coeffs, std::vector<double> values)
{
  std::vector<double> results;
    for (auto const &val:values)
        {
          double result = 0;
          for (int i = 0, deg = coeffs.size() - 1; deg >= 0; deg--, i++)
            {
                result += coeffs[i] * std::pow(val, deg);
            }
          results.push_back (result);
        }
        return results;
}

int main()
{
  std::vector<double> coeffs = { 4, 2, -2, 5};
  std::vector<double> valuesToEvaluate = { 0, 2 , -4};
  std::vector<double> results = Polyval (coeffs, valuesToEvaluate);

  for (auto const &res:results)
    {
      std::cout << res << std::endl;
    }
}

Not sure if there is a better way in terms of performance.

joepol
  • 744
  • 6
  • 22
  • 1
    You can drop `std::pow` for starters, and accumulate the powers of `x` as you go. See the next comment for the name attributed to this method: – Bathsheba Feb 25 '20 at 14:22
  • 1
    Use [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method). – molbdnilo Feb 25 '20 at 14:23
  • Can this help you get there? https://www.boost.org/doc/libs/1_45_0/libs/math/doc/sf_and_dist/html/math_toolkit/toolkit/internals1/rational.html – parktomatomi Feb 25 '20 at 14:30
  • Unless your application has an upper limit of runtime in the order of microseconds, I don't think optimizing a simple for loop will ever have any impact in your code. "Premature optimization is the root of all evil" – Ander Biguri Feb 25 '20 at 14:40
  • Start from 0 degree and accumulate powers of value. No need to use `std::pow()` at all. – vsh Feb 25 '20 at 15:02
  • @parktomatomi I can't use Boost library in my project but it implements Horner's method as suggested by @molbdnilo and @Bathsheba. I now base my solution on Boost's `evaluate_polynomial`, will post my full implementation – joepol Feb 25 '20 at 15:32
  • @AnderBiguri I usually tend to agree with that saying, but I tested the performance of the two methods (same polynomial, same values) 1K times - using `pow()` takes 3 times longer, for this magnitude (a few milliseconds difference). Also see : [Why is pow(int, int) so slow?](https://stackoverflow.com/questions/41072787/why-is-powint-int-so-slow/41072811) – joepol Feb 26 '20 at 09:45

1 Answers1

2

As suggested in the comments I now use Horner's method, based on Boost implementation of polynomial evaluation, The major differences are :

  1. Polynomial order - in this solution, same as Matlab, highest polynomial degree is first. e.g : p(x) = 4x^3 + 2*x^2 -2*x + 5 is represented as a vector like this { 4, 2, -2, 5}

  2. Evaluates multiple values.

#include<assert.h>
#include<vector>

std::vector<double> Polyval(std::vector<double> coeffs, std::vector<double> values)
{
  assert(coeffs.size() > 0);
  std::vector<double> results;
  for (auto const &val:values)
  {
      double result = coeffs[0];
      for (int i = 1; i < coeffs.size(); i++)
        {
            result *= val;
            result += coeffs[i];
        }
        results.push_back (result);
    }
    return results;
}

Edit: Adding performance metrics of the two methods (using pow() vs. Horner's method )

Metrics

Polynomial : p(x) = 4*x^5 + 2*x^4 -2*x^3 + 5*x + 15

Runs : 10,000

Points to evaluate : {0, 2, -4, 8, 15, 1.25, 512 ,-5.3 ,12.215, 153, 23, -11}

Build type : RELEASE

Duration times : pow - 46,576[microseconds] vs. Horner -6,500[microseconds]

Duration difference : ~ 7 times faster (in favor of Horner's method)

Measured duration like this, for both implementations:

#include<iostream>
#include<assert.h>
#include<vector>
#include<chrono>

    int iter = 10000;
    std::vector<double> coeffs = { 4, 2, -2, 5, 0, 15};
    std::vector<double> valuesToEvaluate = {0, 2, -4, 8, 15, 1.25, 512 ,-5.3 
    ,12.215, 153, 23, -11};
    auto start = std::chrono::high_resolution_clock::now();

    do{
    std::vector<double> results = Polyval(coeffs, valuesToEvaluate);
    }
    while(iter-->0);

    auto end = std::chrono::high_resolution_clock::now();
    long duration = std::chrono::duration_cast<std::chrono::microseconds>(end - 
    start).count();

    std::cout << duration << std::endl;
joepol
  • 744
  • 6
  • 22