math/tools/estrin.hpp

0001 /*
0002  * Copyright Thomas Dybdahl Ahle, Nick Thompson, Matt Borland, John Maddock, 2023
0003  * Use, modification and distribution are subject to the
0004  * Boost Software License, Version 1.0. (See accompanying file
0005  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  */
0007 #ifndef BOOST_MATH_TOOLS_ESTRIN_HPP
0008 #define BOOST_MATH_TOOLS_ESTRIN_HPP
0009
0010 #include <array>
0011 #include <vector>
0012 #include <type_traits>
0013 #include <boost/math/tools/assert.hpp>
0014
0015 namespace boost {
0016 namespace math {
0017 namespace tools {
0018
0019 template <typename RandomAccessContainer1, typename RandomAccessContainer2, typename RealOrComplex>
0020 inline RealOrComplex evaluate_polynomial_estrin(RandomAccessContainer1 const &coeffs, RandomAccessContainer2 &scratch, RealOrComplex z) {
0021   // Does anyone care about the complex coefficients, real argument case?
0022   // I've never seen it used, and this static assert makes the error messages much better:
0023   static_assert(std::is_same<typename RandomAccessContainer2::value_type, RealOrComplex>::value,
0024                 "The value type of the scratch space must be the same as the abscissa.");
0025   auto n = coeffs.size();
0026   BOOST_MATH_ASSERT_MSG(scratch.size() >= (n + 1) / 2, "The scratch space must be at least N+1/2");
0027
0028   if (n == 0) {
0029     return static_cast<RealOrComplex>(0);
0030   }
0031   for (decltype(n) i = 0; i < n / 2; i++) {
0032     scratch[i] = coeffs[2 * i] + coeffs[2 * i + 1] * z;
0033   }
0034   if (n & 1) {
0035     scratch[n / 2] = coeffs[n - 1];
0036   }
0037   auto m = (n + 1) / 2;
0038
0039   while (m != 1) {
0040     z = z * z;
0041     for (decltype(n) i = 0; i < m / 2; i++) {
0042       scratch[i] = scratch[2 * i] + scratch[2 * i + 1] * z;
0043     }
0044     if (m & 1) {
0045       scratch[m / 2] = scratch[m - 1];
0046     }
0047     m = (m + 1) / 2;
0048   }
0049   return scratch[0];
0050 }
0051
0052 // The std::array template specialization doesn't need to allocate:
0053 template <typename RealOrComplex1, size_t n, typename RealOrComplex2>
0054 inline RealOrComplex2 evaluate_polynomial_estrin(const std::array<RealOrComplex1, n> &coeffs, RealOrComplex2 z) {
0055   std::array<RealOrComplex2, (n + 1) / 2> ds;
0056   return evaluate_polynomial_estrin(coeffs, ds, z);
0057 }
0058
0059 template <typename RandomAccessContainer, typename RealOrComplex>
0060 inline RealOrComplex evaluate_polynomial_estrin(const RandomAccessContainer &coeffs, RealOrComplex z) {
0061   auto n = coeffs.size();
0062   // Normally, I'd make `ds` a RandomAccessContainer, but its value type needs to be RealOrComplex,
0063   // and the value_type of the passed RandomAccessContainer can just be Real.
0064   // Allocation of the std::vector is not ideal, but I have no other ideas at the moment:
0065   std::vector<RealOrComplex> ds((n + 1) / 2);
0066   return evaluate_polynomial_estrin(coeffs, ds, z);
0067 }
0068
0069 } // namespace tools
0070 } // namespace math
0071 } // namespace boost
0072 #endif