Use Estrin's Scheme for polynomial evaluation

[thomasahle]

According to https://www.boost.org/doc/libs/1_81_0/libs/math/doc/html/math_toolkit/rational.html Boost currently evaluates Polynomials (and rational functions) using Horner's rule.
However, a simple change of the algorithm, known as Estrin's scheme is known to both improve performance and reduce numerical instability.

Maybe there's a reason this method is not currently used?
Otherwise I can provide a pull request.

There are also quite a few other algorithms for polynomial evaluation: https://en.wikipedia.org/wiki/Polynomial_evaluation though most of them requires some preprocessing of the polynomial, so they are only useful if you are going to evaluate it many times. I don't know if something like that would be useful for Boost?

[mborland]

There is no compelling reason we use Horner's rule over any other algorithm to my knowledge. If you provide a pull request I will certainly review it. For performance comparison we typically use google benchmark. You'll find numerous examples in reporting/performance folder.

[NAThompson]

@thomasahle : This is a good idea, and would be a welcome addition by me!

Only thing I worry about is that it's sometimes very hard to get the AVX instructions required by the scheme to generate; probably gotta spend some time with godbolt.

In any case, maybe I should try Estrin's scheme to help with the struggles I'm having in my other MR!

[jeremy-murphy]

I remember reading somewhere that depending on the ratio of some two values, maybe the variable and something else, Horner's method was more accurate when the ratio was above x but another method, perhaps the naive one, was more accurate below x. Now I'm going to have to go back through the books to find it!

[thomasahle]

@NAThompson I can try coding the AVX manually, but from my experience Estrin can often give a factor 2 speedup over Horner, even without AVX. I honestly don't get how it happens simply from rearranging the operations but it works for me at least 😅. Hopefully the Google benchmarks will agree.
For numerical precision I'm not sure how best to test that.

[jzmaddock]

Generally speaking, for simple code like this, trying to second guess the optimizer results in worse code in most cases.

We also have a few different methods for polynomial evaluation already, see https://www.boost.org/doc/libs/1_72_0/libs/math/doc/html/math_toolkit/tuning.html The second order Horner which is the default, is already 1.5-2x faster than naive Horner. BTW Most of the code is machine generated via: https://github.com/boostorg/math/blob/develop/tools/generate_rational_code.cpp

That said, this is all old code now, so may well be worth revisiting and testing for modern processors.

There's also some performance testing in https://github.com/boostorg/math/blob/develop/reporting/performance/test_poly_method.cpp

[NAThompson]

For numerical precision I'm not sure how best to test that.

@thomasahle : We generally use my garbage code ulps_plot.hpp. All the plots I made here are made with that tool. For an introduction to ulps plots, see here.

[thomasahle]

@jzmaddock There's a cpp program for generating the cpp programs? I guess you don't fully trust the compiler after all :D

Where do I see the generated code measured for the tables?

[NAThompson]

Where do I see the generated code measured for the tables?

For instance:

include/boost/math/tools/detail/polynomial_horner3_9.hpp

[NAThompson]

@thomasahle : Just a quick google/benchmark I hacked up for comparison:

//  (C) Copyright Nick Thompson 2023.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <complex>
#include <random>
#include <array>
#include <vector>
#include <iostream>
#include <benchmark/benchmark.h>

template<typename Real, typename RealOrComplex>
inline auto horners_method(const Real* coeffs, long n, RealOrComplex z) {
    RealOrComplex result = coeffs[n - 1];
    for (long i = n - 2; i >= 0; --i) {
        if constexpr (std::is_same_v<Real, RealOrComplex>) {
           result = std::fma(result, z, coeffs[i]);
        } else {
           result = result * z + coeffs[i];
        }
    }
    return result;
}


template<class Real>
void HornersMethodRealCoeffsRealArg(benchmark::State& state)
{
    long n = state.range(0);
    std::random_device rd;
    auto seed = rd();
    std::mt19937_64 mt(seed);
    std::uniform_real_distribution<Real> unif(-10, 10);

    std::vector<Real> c(n);
    std::fill(c.begin(), c.end(), unif(mt));
    Real x = unif(mt);
    for (auto _ : state)
    {
        Real output = horners_method(c.data(), n, x);
        benchmark::DoNotOptimize(output);
        x += std::sqrt(std::numeric_limits<Real>::epsilon());
    }
    state.SetComplexityN(state.range(0));    
}

BENCHMARK_TEMPLATE(HornersMethodRealCoeffsRealArg, float)->Range(1, 1<<15)->Complexity();
BENCHMARK_TEMPLATE(HornersMethodRealCoeffsRealArg, double)->Range(1, 1 << 15)->Complexity();

template<class Real>
void HornersMethodRealCoeffsComplexArg(benchmark::State& state)
{
    long n = state.range(0);
    std::random_device rd;
    auto seed = rd();
    std::mt19937_64 mt(seed);
    std::uniform_real_distribution<Real> unif(-10, 10);

    std::vector<Real> c(n);
    std::fill(c.begin(), c.end(), unif(mt));
    std::complex<Real> x{unif(mt), unif(mt)};
    for (auto _ : state)
    {
        std::complex<Real> output = horners_method(c.data(), n, x);
        benchmark::DoNotOptimize(output);
        x += std::sqrt(std::numeric_limits<Real>::epsilon());
    }
    state.SetComplexityN(state.range(0));    
}

BENCHMARK_TEMPLATE(HornersMethodRealCoeffsComplexArg, float)->Range(1, 1<<15)->Complexity();
BENCHMARK_TEMPLATE(HornersMethodRealCoeffsComplexArg, double)->Range(1, 1 << 15)->Complexity();


template<typename Real, typename RealOrComplex>
inline auto second_order_horner(const Real* coeffs, long n, RealOrComplex z) {
   RealOrComplex p1 = coeffs[n - 1];
   RealOrComplex p2 = coeffs[n - 2];
   RealOrComplex zsq = z*z;
   for (long i = n - 3; i >= 0; i -= 2) {
       p1 = p1 * zsq + coeffs[i];
       p2 = p2 * zsq + coeffs[i-1];
   }
   return p1 + z*p2;
}


template<class Real>
void SecondOrderHornersMethodRealCoeffsRealArg(benchmark::State& state)
{
    long n = state.range(0);
    std::random_device rd;
    auto seed = rd();
    std::mt19937_64 mt(seed);
    std::uniform_real_distribution<Real> unif(-10, 10);

    std::vector<Real> c(n);
    std::fill(c.begin(), c.end(), unif(mt));
    Real x = unif(mt);
    for (auto _ : state)
    {
        Real output = second_order_horner(c.data(), n, x);
        benchmark::DoNotOptimize(output);
        x += std::sqrt(std::numeric_limits<Real>::epsilon());
    }
    state.SetComplexityN(state.range(0));    
}

BENCHMARK_TEMPLATE(SecondOrderHornersMethodRealCoeffsRealArg, float)->Range(1, 1<<15)->Complexity();
BENCHMARK_TEMPLATE(SecondOrderHornersMethodRealCoeffsRealArg, double)->Range(1, 1 << 15)->Complexity();


BENCHMARK_MAIN();

If you'd like to add your Estrin to this it'd be pretty useful I think . . .

[thomasahle]

Ok, I ran the benchmarks now. I got pretty good results, unless I did something wrong, which is very possible. Looks like a factor 2-3 improvement over Horner and Second Order Horner for larger n.

Run on (16 X 2300 MHz CPU s)
CPU Caches:
  L1 Data 32 KiB
  L1 Instruction 32 KiB
  L2 Unified 256 KiB (x8)
  L3 Unified 16384 KiB
Load Average: 2.31, 2.39, 2.54
--------------------------------------------------------------------------------------------------
Benchmark                                                        Time             CPU   Iterations
--------------------------------------------------------------------------------------------------
HornersMethodRealCoeffsRealArg<float>/1                      0.481 ns        0.481 ns   1000000000
HornersMethodRealCoeffsRealArg<float>/8                       3.10 ns         3.10 ns    223132461
HornersMethodRealCoeffsRealArg<float>/64                      31.1 ns         31.1 ns     22366933
HornersMethodRealCoeffsRealArg<float>/512                      448 ns          448 ns      1475651
HornersMethodRealCoeffsRealArg<float>/4096                    3933 ns         3931 ns       174970
HornersMethodRealCoeffsRealArg<float>/32768                  31729 ns        31726 ns        21835
HornersMethodRealCoeffsRealArg<float>_BigO                    0.97 N          0.97 N
HornersMethodRealCoeffsRealArg<float>_RMS                        0 %             0 %
HornersMethodRealCoeffsRealArg<double>/1                     0.258 ns        0.257 ns   1000000000
HornersMethodRealCoeffsRealArg<double>/8                      2.44 ns         2.44 ns    281804677
HornersMethodRealCoeffsRealArg<double>/64                     30.4 ns         30.4 ns     23361289
HornersMethodRealCoeffsRealArg<double>/512                     450 ns          450 ns      1571946
HornersMethodRealCoeffsRealArg<double>/4096                   3919 ns         3919 ns       179403
HornersMethodRealCoeffsRealArg<double>/32768                 32309 ns        32297 ns        20706
HornersMethodRealCoeffsRealArg<double>_BigO                   0.99 N          0.99 N
HornersMethodRealCoeffsRealArg<double>_RMS                       1 %             1 %
SecondOrderHornersMethodRealCoeffsRealArg<float>/1            1.00 ns        1.000 ns    645607563
SecondOrderHornersMethodRealCoeffsRealArg<float>/8            2.15 ns         2.15 ns    310807211
SecondOrderHornersMethodRealCoeffsRealArg<float>/64           26.7 ns         26.7 ns     26769564
SecondOrderHornersMethodRealCoeffsRealArg<float>/512           438 ns          438 ns      1636317
SecondOrderHornersMethodRealCoeffsRealArg<float>/4096         3977 ns         3976 ns       175106
SecondOrderHornersMethodRealCoeffsRealArg<float>/32768       30364 ns        30353 ns        23074
SecondOrderHornersMethodRealCoeffsRealArg<float>_BigO         0.93 N          0.93 N
SecondOrderHornersMethodRealCoeffsRealArg<float>_RMS             1 %             1 %
SecondOrderHornersMethodRealCoeffsRealArg<double>/1           35.4 ns         35.2 ns     20775582
SecondOrderHornersMethodRealCoeffsRealArg<double>/8           2.34 ns         2.30 ns    296232348
SecondOrderHornersMethodRealCoeffsRealArg<double>/64          32.0 ns         30.9 ns     24793506
SecondOrderHornersMethodRealCoeffsRealArg<double>/512          441 ns          441 ns      1388255
SecondOrderHornersMethodRealCoeffsRealArg<double>/4096        3962 ns         3960 ns       180663
SecondOrderHornersMethodRealCoeffsRealArg<double>/32768      31181 ns        31144 ns        22893
SecondOrderHornersMethodRealCoeffsRealArg<double>_BigO        0.95 N          0.95 N
SecondOrderHornersMethodRealCoeffsRealArg<double>_RMS            1 %             1 %
EstrinsMethodRealCoeffsRealArg<float>/1                       2.82 ns         2.82 ns    241334372
EstrinsMethodRealCoeffsRealArg<float>/8                       11.1 ns         11.1 ns     61738195
EstrinsMethodRealCoeffsRealArg<float>/64                      26.1 ns         26.1 ns     27562311
EstrinsMethodRealCoeffsRealArg<float>/512                      139 ns          139 ns      5021485
EstrinsMethodRealCoeffsRealArg<float>/4096                    1068 ns         1067 ns       645542
EstrinsMethodRealCoeffsRealArg<float>/32768                   8637 ns         8632 ns        72112
EstrinsMethodRealCoeffsRealArg<float>_BigO                    0.26 N          0.26 N
EstrinsMethodRealCoeffsRealArg<float>_RMS                        0 %             0 %
EstrinsMethodRealCoeffsRealArg<double>/1                      2.87 ns         2.87 ns    229098074
EstrinsMethodRealCoeffsRealArg<double>/8                      12.5 ns         12.5 ns     54883881
EstrinsMethodRealCoeffsRealArg<double>/64                     27.2 ns         27.2 ns     26123010
EstrinsMethodRealCoeffsRealArg<double>/512                     163 ns          162 ns      4516595
EstrinsMethodRealCoeffsRealArg<double>/4096                   1386 ns         1382 ns       518822
EstrinsMethodRealCoeffsRealArg<double>/32768                 14601 ns        14590 ns        46269
EstrinsMethodRealCoeffsRealArg<double>_BigO                   0.03 NlgN       0.03 NlgN
EstrinsMethodRealCoeffsRealArg<double>_RMS                       1 %             1 %

Where I used this version of Estrin:

template<typename Real, typename RealOrComplex>
inline auto estrin(const Real* coeffs, long n, RealOrComplex z) {
    if (n == 1)
        return coeffs[0];

    RealOrComplex ds[n/2];
    for (int i = 0; i < n / 2; i++) {
        ds[i] = coeffs[2*i] + coeffs[2*i+1] * z;
    }
    n /= 2;

    while (n != 1) {
        z = z * z;
        for (int i = 0; i < n / 2; i++) {
            ds[i] = ds[2*i] + ds[2*i+1] * z;
        }
        n /= 2;
    }

    return ds[0];
}

Note this only supports n which are powers of two, which is what the test uses.
Actually I don't really understand what the SecondOrder method does in the case n=1. Looks like it reads memory out of bounds?

I normally have n as a template argument, to make sure the compiler will unroll all the loops, and not actually allocate an array. But this test seems to use n values much larger than I'm used to (up to 32K).

[NAThompson]

Actually I don't really understand what the SecondOrder method does in the case n=1. Looks like it reads memory out of bounds?

@thomasahle : Yeah, I just hacked this up really quick as an example; I should've actually made a bit of effort . . .

[thomasahle]

Here is a version that doesn't assume n is a power of two. It runs just as fast.

template<typename Real, typename RealOrComplex>
inline auto estrin2(const Real* coeffs, long n, RealOrComplex z) {
    if (n == 1)
        return coeffs[0];

    RealOrComplex ds[(n+1)/2];
    for (int i = 0; i < n / 2; i++) {
        ds[i] = coeffs[2*i] + coeffs[2*i+1] * z;
    }
    if ((n & 1) == 1) ds[n / 2] = coeffs[n - 1];
    n = (n + 1) / 2;

    while (n != 1) {
        z = z * z;
        for (int i = 0; i < n / 2; i++) {
            ds[i] = ds[2*i] + ds[2*i+1] * z;
        }
        if ((n & 1) == 1) ds[n / 2] = ds[n - 1];
        n = (n + 1) / 2;
    }

    return ds[0];
}

I ran the benchmarks on a "Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz" with "-O3 and -match=native".

[NAThompson]

I normally have n as a template argument, to make sure the compiler will unroll all the loops, and not actually allocate an array.

At least in my wavelet MR, I do know the polynomial coefficients at compile time and ostensibly that's where this method would really shine. It just would be a bit awkward in the googlebenchmark file but who cares.

[thomasahle]

If I run this on godbolt with n known in advance, it gives me the same code as a manual Estrin like

  double mult5(double x) {
    // P4(x) = (C0 + C1x) + (C2 + C3x) x2 + C4x4
    // We basically do bottom-up horner
    double c0 = ms[0] + ms[1] * x;  // C0 = M0 + M1 x
    double c1 = ms[2] + ms[3] * x;  // C1 = M2 + M3 x
    double c2 = ms[4];              // C2 = M4
    // Combine first two terms using horner with x^2
    x = x * x;
    double d0 = c0 + c1 * x;        // D0 = C0 + C1 x^2
    double d1 = c2;                 // D1 = C2
    // Combine last term using horner with x^4
    x = x * x;
    double h = d0 + d1 * x;         // H = D0 + D1 x^4
    return h;
  }

I wonder if the lack of compile time unfolding is why the n=8 case above suffers...

[thomasahle]

There are also some methods that does preprocessing on the coefficients, and are able to use just n/2 multiplications. In my experiments these are even faster than Estrin. But of course, it requires a different API.