Add Estrin's method for polynomial evaluation

N.B.: This is a slightly modified version of the code provided by Thomas Dybdahl Ahle in a github issue.

[CI SKIP] [ci skip]

include/boost/math/tools/estrin.hpp

@@ -0,0 +1,58 @@
+/*
+ * Copyright Thomas Dybdahl Ahle, Nick Thompson, Matt Borland, 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)
+ */
+#ifndef BOOST_MATH_TOOLS_ESTRIN_HPP
+#define BOOST_MATH_TOOLS_ESTRIN_HPP
+
+#include <array>
+#include <vector>
+
+namespace boost::math::tools {
+
+template<typename RandomAccessContainer1, typename RandomAccessContainer2, typename RealOrComplex>
+inline auto estrin(RandomAccessContainer1 const & coeffs, RandomAccessContainer2& scratch, RealOrComplex z) {
+  auto n = coeffs.size();
+  for (decltype(n) i = 0; i < n / 2; i++) {
+    scratch[i] = coeffs[2 * i] + coeffs[2 * i + 1] * z;
+  }
+  if (n & 1) {
+    scratch[n / 2] = coeffs[n - 1];
+  }
+  auto m = (n + 1) / 2;
+
+  while (m != 1) {
+    z = z * z;
+    for (decltype(n) i = 0; i < m / 2; i++) {
+      scratch[i] = scratch[2 * i] + scratch[2 * i + 1] * z;
+    }
+    if (m & 1) {
+      scratch[m / 2] = scratch[m - 1];
+    }
+    m = (m + 1) / 2;
+  }
+  return scratch[0];
+}
+
+
+// The std::array template specialization doesn't need to allocate:
+template <typename RealOrComplex1, size_t n, typename RealOrComplex2>
+inline RealOrComplex2 estrin(const std::array<RealOrComplex1, n> &coeffs, RealOrComplex2 z) {
+  std::array<RealOrComplex2, (n + 1) / 2> ds;
+  return estrin(coeffs, ds, z);
+}
+
+template <typename RandomAccessContainer, typename RealOrComplex>
+inline RealOrComplex estrin(const RandomAccessContainer &coeffs, RealOrComplex z) {
+  auto n = coeffs.size();
+  // Normally, I'd make `ds` a RandomAccessContainer, but its value type needs to be RealOrComplex,
+  // and the value_type of the passed RandomAccessContainer can just be Real.
+  // Allocation of the std::vector is not ideal, but I have no other ideas at the moment:
+  std::vector<RealOrComplex> ds((n+1)/2);
+  return estrin(coeffs, ds, z);
+}
+}
+
+#endif

test/test_estrin.cpp

@@ -0,0 +1,79 @@
+#include <random>
+#include <boost/math/tools/estrin.hpp>
+#include "math_unit_test.hpp"
+
+using boost::math::tools::estrin;
+
+template<typename Real, typename Complex>
+void test_polynomial_properties() {
+   std::random_device rd;
+   auto seed = rd();
+   std::mt19937_64 gen(seed);
+   std::uniform_real_distribution<Real> dis(-1, 1);
+   size_t n = seed % 128;
+   std::vector<Real> coeffs(n);
+   for (auto& c : coeffs) {
+        c = dis(gen);
+   }
+   // Test evaluation at zero:
+   auto p0 = estrin(coeffs, Complex(0));
+   if(!CHECK_ULP_CLOSE(p0.real(), coeffs[0], 0)) {
+       std::cerr << "  p(0) != c0 with seed " << seed << "\n";  
+   }
+   // Just given the structure of the calculation, I believe the complex part should be identically zero:
+   if (!CHECK_ULP_CLOSE(p0.imag(), Real(0), 0)) {
+       std::cerr << "  p(0) != c0 with seed " << seed << "\n";
+   }
+
+   auto p1_computed = estrin(coeffs, Complex(1));
+   // Might want to consider using Kahan summation so the expected value is super accurate:
+   auto p1_expected = std::reduce(coeffs.begin(), coeffs.end());
+   if (!CHECK_ULP_CLOSE(p1_expected, p1_computed.real(), coeffs.size())) {
+       std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
+   }
+   if (!CHECK_ULP_CLOSE(Real(0), p1_computed.imag(), coeffs.size())) {
+       std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
+   }
+}
+
+template<typename Real>
+void test_std_array_overload() {
+   std::random_device rd;
+   auto seed = rd();
+   std::mt19937_64 gen(seed);
+   std::uniform_real_distribution<Real> dis(-1, 1);
+   std::array<Real, 12> coeffs;
+   for (auto& c : coeffs) {
+        c = dis(gen);
+   }
+   // Test evaluation at zero:
+   auto p0 = estrin(coeffs, std::complex<Real>(0));
+   if(!CHECK_ULP_CLOSE(p0.real(), coeffs[0], 0)) {
+       std::cerr << "  p(0) != c0 with seed " << seed << "\n";  
+   }
+   if (!CHECK_ULP_CLOSE(p0.imag(), Real(0), 0)) {
+       std::cerr << "  p(0) != c0 with seed " << seed << "\n";
+   }
+   auto p0_real = estrin(coeffs, Real(0));
+   if(!CHECK_ULP_CLOSE(p0_real, coeffs[0], 0)) {
+       std::cerr << "  p(0) != c0 with seed " << seed << "\n";  
+   }
+
+   auto p1_computed = estrin(coeffs, std::complex<Real>(1));
+   auto p1_expected = std::reduce(coeffs.begin(), coeffs.end());
+   if (!CHECK_ULP_CLOSE(p1_expected, p1_computed.real(), coeffs.size())) {
+       std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
+   }
+   if (!CHECK_ULP_CLOSE(Real(0), p1_computed.imag(), coeffs.size())) {
+       std::cerr << " p(1) != sum(coeffs) with seed " << seed << "\n";
+   }
+}
+
+
+int main() {
+   test_polynomial_properties<float, std::complex<float>>();
+   test_polynomial_properties<double, std::complex<double>>();
+   test_std_array_overload<float>();
+   test_std_array_overload<double>();
+   return boost::math::test::report_errors();
+}