More accurate evalpoly

As far as I can tell, you almost implemented a compensated Horner algorithm, only with the notable difference that the error term prod - (hi - pi) is only guaranteed to be exact if abs(prod) < abs(pi). I’m not sure whether this is the assumption you were thinking about but if so, then under this assumption you have a compensated Horner algorithm which should be as accurate as a standard Horner scheme using twice the native precision:
https://www.sciencedirect.com/science/article/pii/S0898122108001120

If this assumption does not hold in practice, maybe we should try and measure how much performance would be lost if we wanted to accurately compute that error term (i.e. the one related to the hi = pi+prod operation) in all cases. Or is that what you already tested and refer to as a “fully correctly double float accurate implementation”?

5 Likes