Polynomial Evaluation

[JonTidswell]

Sleef currently uses a full Estrin method.

Estrin often has lower latency than Horner evaluation, but it sacrifices accuracy and throughput.

Unfortunately, the gains in latency are effectively limited by the ability to exploit the hardware pipeline parallelism, and expanding the polynomial evaluation tree completely can be slower than a partial expansion.
The pipeline is very processor specific, so it would be a large amount of work to cover all the possibilities correctly, so presuming that a full Estrin is justified is a reasonable default position.

However, when there are multiple polynomials that can be evaluated concurrently, then the default position is sub-optimal. Three obvious cases where this occurs are:

1. When unrolling loops, it is better to fill the pipeline with separate Horner evaluations, and not use Estrin.

2. When evaluating rational polynomials, the pipeline needs to be shared between the polynomials.
   Note: Division exacerbates numerical errors, so instead of using a formula like p1(x)/p2(x)
   its better to use a formula like p0(x) + p1(x)/p2(x) where p0(x) is an accurate enough approximation that the errors in the division are below rounding threshold (typically this means p0(x) is evaluated in extended precision).

3. When the input domain is partitioned it may be faster to evaluate multiple partitions in parallel and then to use conditional assignment, rather than use conditional branches.

I have had a first attempt (evalpoly branch) at extending estrin.h to support Horner in order to support a discussion, but I don't assume my approach is suitable.
I also suspect it could be extended to support extended precision in the first couple of terms to avoid the common pattern of a POLY(...) call followed by some instrincs to complete the polynomial in extended precision.

[nibrunie]

This is still an active research topic, see for example Polynomial Evaluation on Superscalar Architecture, Applied to the Elementary Function ex (https://dl.acm.org/doi/abs/10.1145/3408893).
Horner is indeed considered "optimal" for some througput cases because it exposes the least number of operations.

[JonTidswell]

Interesting read.
Although, I think, the variations in accuracy make the results less easy to interpret.

[nibrunie]

I agree. I think this work could be extended with a more thorough evaluation on the accuracy impact and more theoretical / experimental work on evaluation schemes.
I am working with G. revy and C. Lauter towards a generalize framework for polynomial scheme evaluation (including accuracy management)

G. Revy has published work on exploration of fixed point evaluation scheme http://cgpe.gforge.inria.fr/index.php?page=home
and we have been able to publish a draft of our current work here: https://ieeexplore.ieee.org/abstract/document/8825121

I intend to investigate this with my own code generator https://github/com/kalray/metalibm which already implements some tools for polynomial approximation and evaluation scheme generation.
(Performance are behind what sleef provides today for most functions).

[shibatch]

Hello,

Oh, are you one of the authors of this book?

https://www.amazon.com/dp/3319765256

[JonTidswell]

I was hoping for a generic approach, not focused on any particular function.
(Part of why I chose Giles erfinv approximation as an introduction is because it has clear fast/slow paths, and is a simple polynomial not rational; this is not uniformly true of functions in which I am interested.)

[shibatch]

As for point 2, I don't understand what formula with a division you are talking about. You seems not using a division in your erfinv implementation.

As for POLY call, are you suggesting to incorporate the double-double evaluation of the lowest-degree terms? I think that would just complicate things, but I would like to know if you have a good idea.

[shibatch]

If there are two independent polynomials to evaluate, Horner is almost as fast as Estrin with the current typical hardware. I agree with this point.

[JonTidswell]

I am pondering incorporating double-double (or equally double-float) for the lowest degree terms.
I agree this will require slightly more complicated macro definitions, but it should correspondingly simplify function definitions.

I like the idea of this trade-off because I want to focus on choices in polynomials required to keep the pipeline full, not about the intrinsics required to evaluate them, so I want to implement and benchmark multiple variations; which I think I can do better if I abstract this detail into a macro.

[ PS As an example of a formula with a division consider any (every?) Inverse CDF of Normal (a.k.a probit), although William Shaw's (https://arxiv.org/abs/0901.0638) is the only one targeted for SIMD. ]

[shibatch]

Sounds good.
But please note that there are differences in the ways the lowest-degree terms are evaluated in DD precision.
For example, if the coefficient is 1, the multiplication for that term is omitted.